一文看懂 Facebook 开源的时序预测框架算法

分享于2021年05月13日 算法
网易严选有很多业务决策依赖预测模型的输入,其中时间序列预测是一类比较基础的模型,服务于采购、营销、仓配、客服等业务。在客服话务量预测和单量预测场景中我们用到了开源时序预测框架Prophet。本篇介绍Prophet的基本原理和使用方法。1. 简介2. 安装3. 基本用法4. 算法原理4.1 趋势模型4.2 季节模型4.3 节假日模型4.4 Prophet模型5. 源码5.1 模型初始化5.2 模型训练5.3 生成测试集5.4 预测6. 总结7. 参考1. 简介Prophet[1]是 Facebook 在2017年开源的时间序列预测算法,提供了 R 和 Python 语言的支持[2]。Prophet容易上手,短短几行代码就能建立时序预测模型,算法的基本思想类似于时间序列分解,将时间序列分成为趋势(Trend),季节性(seasonality)和节日(holiday)。Prophet能自动识别一些简单周期(年/周/日),并内置了节假日模块,甚至支持了中国节日(包括农历节日,如端午节、中秋节)。此外,它还能自动剔除异常值和处理缺失值。Prophet提供了大量可配置的参数,使用者可根据具体需求调整模型,比如调整季节性的拟合度,添加自定义节假日,添加自定义变量等。不同的参数配置可能得到完全不一样的模型,所以使用Prophet时需要对自定义参数有一些了解。2. 安装通过Anaconda安装Prophet非常简单(用pip安装可能会踩坑):# 1. 创建虚拟环境conda create -n prophet python=3.7# 2. 激活虚拟环境source activate prophet# 3. 通过conda-forge安装prophetconda install -c conda-forge fbprophet3. 基本用法Prophet的基本用法分以下几步:读取训练数据初始化prophet模型训练prophet模型生成测试集进行预测# 导入包import pandas as pdfrom fbprophet import Prophet# 读数据# df(dataframe)需要包含两列,一列 date,一列 ydf = pd.DataFrame({    'ds': pd.date_range('20071210',periods=3700),    'y': [np.random.rand() + np.cos(0.05 * x) for x in range(3700)]})# 初始化 prophet 模型m = Prophet()# 训练 prophet 模型m.fit(df)# 生成预测数据(默认天粒度)future = m.make_future_dataframe(periods=365)# 使用训练好的模型进行预测forecast = m.predict(future)# 预测绘图fig1 = m.plot(forecast)# 分量绘图fig2 = m.plot_components(forecast)数据示例输入和输出数据格式如下:输入数据 df: pd.DataFrame-dsy02007-12-109.59076112007-12-118.51959022007-12-128.18367732007-12-138.07246742007-12-147.893572预测结果  forecast: pd.DataFramedsyhatyhat_loweryhat_upper32652017-01-158.2129427.4635608.93721532662017-01-168.5379937.7902599.26749232672017-01-178.3254287.5256759.05939132682017-01-188.1580597.4336348.88362732692017-01-198.1700467.4318018.840703更多的参数配置请参考官方示例[3]。4. 算法原理Prophet基于时间序列的分解,把时间序列分解成四个部分,分别是趋势项,季节项,节日项和剩余项,将这几项之和作为预测值:其中表示趋势项,它表示时间序列非周期性的变化趋势; 表示周期项,或者称为季节项,一般来说是以日或周或者年为单位; 表示节假日项,表示在当天是否存在节假日; 表示误差项或者称为剩余项。Prophet 算法就是通过拟合这几项,然后最后把它们累加起来就得到了时间序列的预测值。4.1 趋势模型在 Prophet 中,趋势项可以有两种选择,一个是饱和增长模型(saturating model),也称逻辑回归模型(logistic function);另一个是分段线性模型(piecewise linear model)。下面分别介绍这两个模型:饱和增长模型在自然界中,很多趋势是非线性增长的,比如人口,经济增速等。它们通常都有一定的承载上界,比如facebook的用户数量,这种增长常用 Logistic 增长模型进行建模。Logistic 函数的基本形式是:函数图像如下图所示:在实际情况中,拟合趋势需要在 Logistic 模型上增加一些参数:其中称为曲线的最大渐近值(也称承载上界),表示曲线的增长率, 表示曲线的中点。下图为发生变动时的函数变化。当时,恰好就是 sigmoid 函数的形式。渐进值可能会随时间发生变化,所以将替换为随时间变化的函数。增长率也可能会随时间变化,例如特定的事件会改变增长率。为了适应增长率的变化,Prophet引入变点(changepoint)的概念,变点前后时序的趋势发生了变化。Prophet默认会设置25个变点(n_changepoints=25),设置范围是前80%的时序数据(changepoint_range=0.8),也就是在时间序列的前 80% 的区间内会设置变点,而且是均匀分布,如下图所示。假设在时间序列中设置个变点,每个变点的时间戳为。定义在变点处增长率的变化率为(可以理解为曲线的二阶导数,增长率的一阶导数),表示变点在时间上增长率的变化,若定义基准变化率为,那么时间上的增长率可以定义为:为时间之前变点的增长变化率累加和,再加上基准变化率。使用向量的形式重写上述公式,简化为:其中,增长率变化的时候,偏移量也必须跟着变化,否则函数会不连续。如图所示,纵轴表示时序趋势值,横轴表示时间。假设一开始的趋势变化如紫线所示,增长率,偏移量,时出现变点,增长率变换为,如红线所示。但是在变点处,也就是处红紫线不相交,造成趋势不连续。可通过修正偏移量使变点前后趋势连续。如下图所示,当偏移量修正为0.8时,变点前后趋势项连续。在变点偏移量的更新可以表示为:最后logistic model表示为:注意:在选择逻辑回归模型作为趋势项时,需要人为指定。分段线性模型对于那些不显式为饱和增长的预测问题,Prophet还提供了另外一种简洁且有用的趋势模型 — 分段线性模型。变点之间的趋势是线性函数,如下图所示。线性函数指的是, 而分段线性函数指的是在每一个子区间上,函数都是线性函数,但是在整段区间上,函数并不完全是线性的。和逻辑回归模型相同,分段线性模型同样加入了变点,分段线性模型可表示为:其中,表示初始增长率,表示时间的增长率,为修正后的偏移量,其中,的设定和逻辑回归模型有所区别。4.2 季节模型通常时间序列包含多周期季节性(multi-period seasonality),比如寒暑假的影响以年为周期,工作日的影响以周为周期等。Prophet采用傅里叶级数来拟合周期效应。傅里叶级数能将任何周期函数或周期信号分解成一个(可能由无穷个元素组成的)简单振荡函数的集合,即正弦函数和余弦函数。举一个例子,以为周期的时间序列可以用下式傅里叶级数表示:可视化效果如下:越大,也就是级数越多,越趋向于对周期的过拟合,在Prophet中,可以通过参数调整傅里叶级数。假设一个以为周期的时间序列,它的傅里叶级数形式为:Prophet内置了三种周期:年、周和日。当时,表示以年为周期,当时,表示以周为周期,当时,表示以天为周期。Prophet 作者按照经验,设置了这些周期的默认级数,对于以年为周期的序列, ;对于以周为周期的序列, ,对于以天为周期的序列, 。使用傅里叶级数拟合周期性需要估计个参数,用向量表示:对于年周期()来说,它的季节性傅里叶特征可以表示为:时间的季节性分量就可以表示为:4.3 节假日模型节假日和事件往往对时间序列预测有较大的影响,且通常不遵循周期模式,所以无法用季节性模型去建模节日。另外,节日和时间也有可能是非稳定周期的,比如中国的农历节日,NBA的决赛时间,手机的新品发布时间等等。不同的节假日可以看成相互独立的模型,且不同的节假日有不同的前后窗口值,表示该节假日会影响前后一段时间的时间序列。比如春节有15天的窗口,而清明节只有3天窗口期。对于第个节假日来说,表示该节假日过去和未来时间的集合,比如下表列举的是感恩节的集合。用参数表示节日的分量,向量表示个节假日的分量。用指示函数表示在时间所处的节日,类似于one-hot形式:时间的节日分量可以表示为其中,服从正态分布,其中是可调参数,默认值是10。当值越大时,表示节假日对模型的影响越大;当值越小时,表示节假日对模型的效果越小。4.4 Prophet模型综上,Prophet的预测模型[4]可以表示为:若选择逻辑回归模型作为趋势模型,上式可演变为:若趋势项为分段线性模型:训练时,采用L-BFGS做最大似然估计(maximum likelihood estimate)训练出模型中的参数。5. 源码若想完全掌握Prophet模型,还是需要去撸一遍源码,一方面是在使用时心里有个底,另一方面是为了更好的调参。通过 debug 以下代码走一遍 Prophet 的流程。# 初始化 prophet 模型m = Prophet()# 训练 prophet 模型m.fit(df)# 生成预测数据(默认天粒度)future = m.make_future_dataframe(periods=365)# 使用训练好的模型进行预测forecast = m.predict(future)5.1 模型初始化首先,初始化prophet模型:m = Prophet()初始化的时候可以指定一些参数,若不指定,Prophet会指定参数的默认值。具体的参数默认初始化如下:# func __init__()# 初始化参数(可设置)growth='linear',              # 趋势选择 'linear'或'logistic'(分段线性 或 逻辑回归)changepoints=None,            # 变点时间(Series)n_changepoints=25,            # changepoint个数:默认25changepoint_range=0.8,        # changepoint范围:默认前80%yearly_seasonality='auto',    # 年周期性,可设 True or Falseweekly_seasonality='auto',    # 周周期性,可设 True or Falsedaily_seasonality='auto',     # 日周期性,可设 True or Falseholidays=None,                # 节假日,默认为Noneseasonality_mode='additive',  # 周期叠加方式,'additive'或'multiplicative'seasonality_prior_scale=10.0, # 周期先验值holidays_prior_scale=10.0,    # 节假日先验值changepoint_prior_scale=0.05, # 变点先验值mcmc_samples=0,               # interval_width=0.80,          #uncertainty_samples=1000,     # # 其他参数start = None                   # 训练集开始时间y_scale = None                 # 训练集 y 最大值logistic_floor = Falset_scale = None                 # 训练集时间跨度changepoints_t = None          # 变点在时间轴的比例seasonalities = OrderedDict({})extra_regressors = OrderedDict({})country_holidays = None       # 指定国家节日,中国为'CN'stan_fit = None               # pystan拟合时的参数params = {}                   # 训练得到的参数history = None                # 输入数据history_dates = None          # 历史数据中的时间train_component_cols = None   # 每列特征对应的类别component_modes = None        # {'additive':[],'multiplicative':[]}train_holiday_names = None    # 配置的 holiday 名称其他参数都比较好理解,重点说一下后缀带有prior_scale的参数。它们是控制变量的拟合程度的参数,可以理解为正则项系数,值越高,越趋向过拟合,值越低越趋向平滑。初始化时,Prophet 会验证以下参数:func validate_inputs() # 验证输入1. growth not in ('linear', 'logistic') # 验证'growth'2. (changepoint_range < 0) or (changepoint_range > 1) # 验证changepoint_range区间3. holidays lower_window/upper_window # 验证holiday名称和窗口值4. seasonality_mode not in ['additive', 'multiplicative'] # 验证'seasonality_mode'以上,初始化的部分就结束了。5.2 模型训练接下来便可以喂入数据训练模型了:m.fit(train_data)其中,train_data是 pd.DataFrame 类型,必须包含两列:dsy日期值如果growth设为'logistic',必须还有一列cap:dsycap日期值承载上界这里的cap其实就是趋势模型中的,需要人为输入。喂入数据后,开始训练Prophet模型,fit() 函数中的主要步骤如下:扩充dataframe设置季节性参数生成季节性特征设置变点初始化变点参数初始化stan模型(prophet采用 pystan[5] 作为线性回归模型)训练stan模型简单来说,这些步骤可以合并成最基本的机器学习范式:1-5步生成特征矩阵(如下图所示),6-7步训练线性回归模型(采用pystan)。1. 扩充dataframe首先通过 setup_dataframe() 函数为 train_data 增加三列,分别对时间和训练数据进行归一化。df['floor']=0 # 最小值(默认)df['t'] = (df['ds'] - self.start) / self.t_scale # 时间归一化df['y_scaled'] = (df['y'] - df['floor']) / self.y_scale # y 值归一化其中t_scale = df['ds'].max() - df['ds'].min() # 最大时间和最小时间的差值y_scale = (df['y'] - floor).abs() # y最大值的绝对值train_data:2. 设置季节性参数若初始化模型时没有设置季节性参数,Prophet 会通过 set_auto_seasonalities() 自动检测训练数据并设置季节性参数,主要有以下几个参数:min_dt: 检测最小时间间隔yearly_disable:若训练数据时间小于两年,禁用年周期weekly_disable:若训练数据时间小于两周或min_dt>=一周,禁用周周期daily_disable:若训练数据时间小于两天或min_dt>=一天,禁用天周期fourier_order:parse_seasonality_args(): 获取内置季节性的傅里叶级数(见下图)若季节性没有disable,设置默认季节性参数:# prophet 内置的 [年/月/日] 周期季节性参数:# period:周期# fourier_order:傅里叶级数# prior_scale:先验值,用于调整拟合度# mode:加法或乘法模型self.seasonalities['yearly'] = {        'period': 365.25,        'fourier_order': fourier_order, # 默认为 10        'prior_scale': self.seasonality_prior_scale, # 10        'mode': self.seasonality_mode, # additive        'condition_name': None    }self.seasonalities['weekly'] = {        'period': 7,        'fourier_order': fourier_order, # 默认为 3        'prior_scale': self.seasonality_prior_scale, # 10        'mode': self.seasonality_mode, # additive        'condition_name': None    }self.seasonalities['daily'] = {        'period': 1,        'fourier_order': fourier_order, # 默认为 4        'prior_scale': self.seasonality_prior_scale, # 10        'mode': self.seasonality_mode, # additive        'condition_name': None    }3. 生成季节性特征设定季节性参数后,通过 make_all_seasonality_features() 函数生成包含季节性特征的dataframe,包括季节特征、节日特征和自定义特征。make_all_seasonality_features的代码结构如下所示:make_seasonality_features(): 生成季节性特征fourier_series():生成傅里叶级数特征make_holiday_features():生成节日特征make_holidays_df(): 根据国家和年份构建节假日dataframeconstruct_holiday_dataframe(): 构建节假日dataframemake_holiday_features():构建节假日特征additional regressors():生成自定义特征regressor_column_matrix()add_group_component返回:seasonal_features, prior_scales, component_cols, modes首先,通过 make_seasonality_features 构建季节性特征,这里主要采用傅里叶级数构建特征(原理见第4节)。代码如下:# fourier_series: 生成傅里叶级数特征np.column_stack([    fun((2.0 * (i + 1) * np.pi * t / period)) # t 为 1970-01-01 到现在的天数 array    for i in range(series_order)    for fun in (np.sin, np.cos)])下面两张图展示了年和周两种周期的傅里叶级数特征:year_seasonality: 维的dataframe,每一行是该时间的20维的 year_fourier 特征。week_seasonality: 维的dataframe,每一行是该时间的6维的 week_fourier 特征。其次, make_holiday_features 函数根据节日配置参数构建节假日dataframe:holiday_df,如下表所示:并根据节假日表生成对应的节假日特征。holiday_feature: 维的特征,类似于one-hot表示。另外,Prophet还支持加入自定义特征,源码中 additional_regressors 函数绑定自定义特征。最后,合并seasonality、holiday、addtional特征,生成总的季节性特征seasonal_feature:regressor_column_matrix 函数对 seasonal_feature 进行加工,生成每一维特征对应的 one-hot 向量,方便之后的模型解释和可视化:modes:函数make_all_seasonality_features的返回结果,它是一个记录加法变量和乘法变量的字典 (dict) 。: {'additive': [  'yearly',  'weekly',  'spring_festival',  'double_11',  'double_12',  'yx_411',  'yx_618',  \"New Year's Day\",  'Chinese New Year',  'Tomb-Sweeping Day',  'Labor Day',  'Dragon Boat Festival',  'Mid-Autumn Festival',  'National Day',  'additive_terms',  'extra_regressors_additive',  'holidays'],'multiplicative': [ 'multiplicative_terms', 'extra_regressors_multiplicative']}4. 设置变点set_changepoint: 变点的设置有以下几种情况。参数显式传递变点自动生成变点不使用变点如果自动设置变点,Prophet会默认在前80%的数据中选择25个变点。changepoints 示例如下:5. 初始化变点参数linear_growth_init:分段线性增长模型,初始化斜率和截距logistic_growth_init:logistic增长模型,初始化增长率和偏移量下图为实际计算过程中和的初始化值。6. 初始化stan模型Prophet预置了一份模型,初始化时会加载进来:model = prophet_stan_model以下是 stan 模型所用到的参数:dat = {    'T': history.shape[0],    'K': seasonal_features.shape[1],    'S': len(self.changepoints_t),    'y': history['y_scaled'],    't': history['t'],    't_change': self.changepoints_t,    'X': seasonal_features,    'sigmas': prior_scales,    'tau': self.changepoint_prior_scale,    'trend_indicator': int(self.growth == 'logistic'),    's_a': component_cols['additive_terms'],    's_m': component_cols['multiplicative_terms'],}7. 训练stan模型这里简单介绍一下 Stan。Stan其实是一种用于统计推断的概率编程语言。以线性回归为例,下面是最简单的线性回归模型,具有单个预测变量,斜率和截距系数以及正态分布的噪声。该模型可以编写为:上式等效于对于残差的采样:并进一步简化为:Prophet的 stan 模型可直接写为:# stan分段线性模型# y ~ 趋势项 * (1 + 乘法项)+ 加法项(季节、节假日) + 残差y ~ normal(  linear_trend(k, m, delta, t, A, t_change)  .* (1 + X * (beta .* s_m))  + X * (beta .* s_a),  sigma_obs)vector linear_trend(  real k, # 斜率  real m, # 截距  vector delta, # 增长率变化率向量  vector t, # 时间  matrix A, # 标志位  vector t_change ) {  return (k + A * delta) .* t + (m + A * (-t_change .* delta));}5.3 生成测试集通过 make_future_dataframe 函数在训练集基础上往后延 periods 的时间作为测试集。future = m.make_future_dataframe(periods=90)5.4 预测forecast = m.predict(future)第一步,通过setup_dataframe()函数初始化测试集dataframe。第二步,对测试集进行趋势预测,这里展示分段线法(piecewise_linear)的预测过程。再回顾一下分段线性模型的公式:源码中的趋势计算如下:# 梯度和截距累加# Intercept changesgammas = -changepoint_ts * deltas# Get cumulative slope and intercept at each tk_t = k * np.ones_like(t) # k 为初始斜率m_t = m * np.ones_like(t) # m 为初始截距for s, t_s in enumerate(changepoint_ts):    indx = t >= t_s # 大于t_s的部分进行累加    k_t[indx] += deltas[s] # 累加斜率    m_t[indx] += gammas[s] # 累加截距trend = k_t * t + m_ttrend = trend * self.y_scale + df['floor'] # scale还原通过累加计算每一个变点之间的和,最后乘以缩放因子 y_scale,预测出趋势。第三步,周期性预测:通过 predict_seasonal_components 函数预测 seasonality, holidays 和 added regressors。与训练过程类似,通过 make_all_seasonality_features 函数构建预测日期的季节性特征,生成类似于训练集的特征矩阵。计算lower_p, upper_p,也就是预测结果的上下界。for component in component_cols.columns:    beta_c = self.params['beta'] * component_cols[component].values    comp = np.matmul(X, beta_c.transpose())    if component in self.component_modes['additive']:        comp *= self.y_scale    data[component] = np.nanmean(comp, axis=1)    data[component + '_lower'] = np.nanpercentile(        comp, lower_p, axis=1,    )    data[component + '_upper'] = np.nanpercentile(        comp, upper_p, axis=1,    )计算每个周期性特征分量 seasonal_components第四步,计算趋势和yhat的上下界(uncertainty)。sample_posterior_predictive:第五步,计算y_hat。df2 = pd.concat((df[cols], intervals, seasonal_components), axis=1)df2['yhat'] = (        df2['trend'] * (1 + df2['multiplicative_terms'])        + df2['additive_terms'])下图为预测结果的可视化示例,分别展示了trend、holidays、weekly、yearly分量,提供了一定的可解释性。6. 总结Prophet 针对一维时间序列进行建模,将时间序列分为趋势、季节、节日三部分特征,并输出预测结果的上下界,同时保证了预测结果的可解释性。从原理上来说,Prophet是一个自动特征提取的线性回归模型,适用于少特征、规律性高、单维度的时间序列场景。若业务场景有多维度特征,多时间序列,复杂度高,则不适用Prophet。7. 参考[1]: Prophet官网:https://facebook.github.io/prophet/docs/quick_start.html#python-api[2]: Prophet代码地址:https://github.com/facebook/prophet[3]: Prophet示例:https://github.com/facebook/prophet/tree/master/notebooks[4]: Prophet预测模型论文:https://peerj.com/preprints/3190/[5]: stan模型:https://pystan.readthedocs.io/en/latest/[6]: Facebook 时间序列预测算法 Prophet 的研究. https://zhuanlan.zhihu.com/p/52330017- EOF -推荐阅读  点击标题可跳转1、Facebook 工程师总结的 14 种算法面试模式2、Facebook 首次开源图片/视频匹配算法3、经典算法面试题:有序矩阵的快速查找觉得本文有帮助?请分享给更多人推荐关注「算法爱好者」,修炼编程内功点赞和在看就是最大的支持❤️