1. Prophet简介
Facebook 在 2017 年开源了一个叫 fbprophet 的时间序列预测的算法,Facebook 所提供的 prophet 算法Prophet 是 Facebook 发布的基于可分解(趋势+季节+节假日)模型的开源库。该算法支持自定义季节和节假日,解决了像春节、618 和双十一这种周期性节假日的指标预测难题。prophet 不仅可以处理时间序列存在一些异常值的情况,也可以处理部分缺失值的情形,还能够几乎全自动地预测时间序列未来的走势。而且 Prophet 包提供了直观易调的参数,即使是对缺乏模型知识的人来说,也可以据此对各种商业问题做出有意义的预测。它的官方网址与基本介绍来自于:
Github:github.com/facebook/prophet
官方网址:facebook.github.io/prophet
2. Prophet算法原理
数据的输入和输出
Prophet 的输入包含两列数据:ds 和 y 。ds 列为日期(YYYY-MM-DD)或者是具体的时间点(YYYY-MM-DD HH:MM:SS)。y 列是数值变量,即预测量。
算法实现
一般来说,在实际生活和生产环节中,除了季节项,趋势项,剩余项之外,通常还有节假日的效应。所以,在 prophet 算法里面同时考虑了以上四项:
$$y(t)=g(t)+s(t)+h(t)+\epsilon(t)$$
$g(t)$ 表示趋势项,它表示时间序列在非周期上面的变化趋势;
$s(t)$ 表示周期项,或者称为季节项,一般来说是以周或者年为单位;
$h(t)$ 表示节假日项,表示在当天是否存在节假日;
$\epsilon(t)$ 表示误差项或者称为剩余项。Prophet 算法就是通过拟合这几项,然后最后把它们累加起来就得到了时间序列的预测值。
趋势项模型 $g(t)$
趋势是对时间序列中非周期性部分的趋势进行拟合。在 Prophet 算法里面,趋势项有两个重要的函数,一个是基于逻辑回归函数,另一个是基于分段线性函数。这里来介绍一下基于逻辑回归的趋势项是怎么做的。
逻辑回归函数形式为:$\sigma(x)=\frac{1}{1+e^{-x}}$, 如果增加一些参数的话,那么逻辑回归就可以改写成:$f(x)=\frac{C}{1+e^{-k(x-m)}}$,这里的 $C,k,m$ 分别为曲线的最大渐近值,曲线的增长率,曲线的中点。当 $C=1,k=1,m=0$ 时,恰好就是大家常见的 sigmoid 函数的形式。那么这里增加了参数的一般函数形式就为:
$$f(x)=\frac{C}{1+e^{-k(x-m)}}$$
在现实环境中,参数 $C,k,m$ 不可能都是常数,而很有可能是随着时间的迁移而变化的,因此,在 Prophet 里面,作者考虑把这三个参数全部换成了随着时间而变化的函数,也就是 $C=C(t),k=k(t),m=m(t)$
通常在生产生活中,我们希望某指标在整个预测区间内持续增长或下降。但也会有最大最小容量的限制,比如某 app 预测在某地未来 12 个月的 UV,最大 UV 只能接近该地区人数。因此分析师在调参时可以定义时间序列预测的容量限制为 $C(t)$。

图1
另外,在现实的时间序列中,指标经常会因为一些事件的发生,曲线的走势速率不会是恒定的,在某些特定的时候会发生变化,那么需要在一些变点改变其速率。例如在下图中,点线代表给定时间序列中的突变点。

图 2
在 Prophet 里面,是需要设置变点的位置的,而每一段的趋势和走势也是会根据变点的情况而改变的。在程序里面有两种方法,一种是通过人工指定的方式指定变点的位置;另外一种是通过算法来自动选择。因此 Prophet 算法需要给出变点的位置,个数,以及增长的变化率的。即changepoint_range,n_changepoint,changepoint_prior_scale 三个参数需要手动设置或自动默认算法的设置。changepoint_range 指的是百分比,需要在前 changepoint_range 那么长的时间序列中设置变点,默认是 changepoint_range = 0.8。n_changepoint 表示变点的个数,默认是 n_changepoint = 25。changepoint_prior_scale 表示变点增长率的分布情况,$\delta_j \sim Laplace(0,\tau)$ ,这里的就是 change_point_scale。默认为 0.05。因此总结起来就是变点的选择是基于时间序列的前 80% 的历史数据,然后通过等分的方法找到 25 个变点,而变点的增长率是满足 Laplace 分布 $\delta_j \sim Laplace(0,0.05)$ 的。即这三个参数决定了 $k(t)$ 和 $m(t)$。两个函数的推导过程如下:
假设已经放置了 S 个变点了,即 n_changepoint=S,并且变点的位置是在时间戳 $s_j,1\leq j \leq S$ 上。在这些时间戳上,增长率的变化 changepoint_prior_scale,为向量 $\delta \in R^S$ 其中 $\delta_j$ 表示在时间戳 $s_j$ 上的增长率的变化量,服从拉普拉斯分布。一开始的增长率使用 $k$ 来代替,那么在时间戳 $t$ 上的增长率就是 $k+\sum_{j:t>s_j} \delta_j$ ,通过 indicator 函数 $a(t) \in \{0,1\}^S$ 表示为:
$$a_j(t)=\begin{cases} 1,\,\, if \,\,t \geq s_j \\ 0,\,\, otheerwise \end{cases}$$
那么在时间戳 $t$ 上面的增长率就是 $k+a^T\delta$. 一旦变化量 $k$ 确定了,另外一个参数 $m$ 也要随之确定。在这里需要把线段的边界处理好,因此通过数学计算可以得到:
$$\gamma_j=(s_j-m-\sum_{l < j}\gamma l) * 1-\frac{k+\sum_{l < j}\delta_l}{k+\sum_{l \leq j}\delta_l}$$
所以,原函数 $f(x)=\frac{C}{1+e^{-k(x-m)}}$ 中的常数换成关于时间的函数形式时 $C=C(t),k=k(t),m=m(t)$ 就变成了分段的逻辑回归增长模型:
$$f(x)=\frac{C(t)}{1+e^{-(k+a(t)^t\delta)(t-(m+a(t)^T\gamma}}$$
因此,当 $\tau$ 趋近于零的时候,$\delta_j$ 也是趋向于零的,此时的增长函数将变成全段的逻辑回归函数或者线性函数。
季节项模型 $s(t)$
周期性的变化因子是时间序列预测模型都会考虑的因素,为了拟合并预测季节的效果,Prophet 基于傅里叶级数提出了一个灵活的模型。季节效应 $S(t)$ 根据以下方程进行估算: $P$ 表示时间序列的周期,$P=365.25$ 表示以年为周期,$P=7$ 表示以周为周期。季节效应 $S(t)$ 傅立叶级数形式是:
$$s(t)=\sum^N_{n=1}(a_ncos(\frac{2\pi nt}{P})+b_nsin(\frac{2\pi nt}{P}))$$
对季节性建模时,需要在给定 N 的情况下,估计参数 $\beta=(a_1,b_1,...,a_N,b_N)^T$ 傅里叶阶数 N 是一个重要的参数,它用来定义模型中是否考虑高频变化。对时间序列来说,如果分析师认为高频变化的成分只是噪声,没必要在模型中考虑,可以把 N 设为较低的值。如果不是,N 可以被设置为较高的值并用于提升预测精度。就作者的经验而言,对于以年为周期的序列 $P=365.25$ 而言,$N=10$;对于以周为周期的序列 $P=7$ 而言,$N=3$ 。
当 $N=10$ 时,
$$X(t)=[cos(\frac{2\pi(1)t}{365.25}),...,sin(\frac{2\pi(1)t}{365.25}))$$
当 $N=3$ 时,
$$X(t)=[cos(\frac{2\pi(1)t}{7}),...,sin(\frac{2\pi(1)t}{7}))$$
因此,时间序列的季节项就是:$s(t)=X(t)\beta$ 而 $\beta \sim N(0,\sigma^2)$ 。这里的 \sigma$ 值越大,表示季节的效应越明显;这个值越小,表示季节的效应越不明显。同时,在代码里面,seasonality_mode 也对应着两种模式,分别是加法和乘法,默认是加法的形式。
节假日项模型 $h(t)$
现实生活中的预测场景中有很多节假日,而且不同的国家有着不同的假期。还有类似于 618、双十一等这样不列入官方节日,但是对于指标预测影响非常重要的日期。在 Prophet 里面,通过维基百科里面对各个国家的节假日的描述,hdays.py 收集了各个国家的特殊节假日。Prophet 还允许分析师使用过去和未来事件的自定义列表,例如印度的 The Super Bowl,国内的双十一等。
由于每个节假日对时间序列的影响程度不一样,例如春节,国庆节则是七天的假期,对于劳动节等假期来说则假日较短。因此,不同的节假日可以看成相互独立的模型,并且可以为不同的节假日设置不同的前后窗口值,表示该节假日会影响前后一段时间的时间序列。用数学语言来说,对与第 $i$ 个节假日来说,$D_i$ 表示该节假日的前后一段时间。为了表示节假日效应,我们需要一个相应的 indicator 函数,同时需要一个参数 $k_i$ 来表示节假日的影响范围。假设我们有 $L$ 个节假日,那么
$$h(t)=Z(t)k=\sum_{i=1}^L k_i 1_{\{ t \in D_i\}} $$
其中 $Z(t)=(1_{\{ t \in D_1\}},...,1_{\{ t \in D_L\}}),k=(k_1,...,k_L)^T$
其中 $k \sim N(0,\nu^2)$ 该正态分布的标准差默认值是 10,当值越大时,表示节假日对模型的影响越大;当值越小时,表示节假日对模型的效果越小。用户可以根据自己的情况自行调整。
模型拟合
结合上面对增长项,季节项,节假日项三方面的详细讲解,现在可以用线性将三块结合一起来拟合时间序列:
$$y(t)=g(t)+s(t)+h(t)+\epsilon_t$$
在 Prophet 中,可以使用 Prophet 默认的参数,也可以自己设置以下四种参数:
Capacity:在增量函数是逻辑回归函数的时候,需要设置的容量值。
Change Points:可以通过 n_changepoints 和 changepoint_range 来进行等距的变点设置,也可以通过人工设置的方式来指定时间序列的变点。
季节性和节假日:可以根据实际的业务需求来指定相应的节假日。
光滑参数: $\tau=changepoint_prior_scale$ 可以用来控制趋势的灵活度, $\sigma=seasonality_prior_scale$ 用来控制季节项的灵活度, $\nu=holidays prior scale$ 用来控制节假日的灵活度。
1. 趋势参数

图 3 趋势参数
2. 季节&假日参数

图 4 季节&假日参数
3. Prophet实战(附Python代码)
该数据集是DATAHACK平台上的一个单变量的时间序列,即某新型公共交通服务的每小时客运量。基于给定的过去25个月的历史交通流量数据,我们可以尝试预测未来七个月的交通情况。
#import data import pandas as pd import numpy as np from fbprophet import Prophet #Read train and test train = pd.read_csv('Train_SU63ISt.csv') test = pd.read_csv('Test_0qrQsBZ.csv') #Convert to datetime format train['Datetime'] = pd.to_datetime(train.Datetime,format='%d-%m-%Y %H:%M') test['Datetime'] = pd.to_datetime(test.Datetime,format='%d-%m-%Y %H:%M') train['hour'] = train.Datetime.dt.hour

图 5 原始数据时间序列
我们可以看到时间序列中有很多噪声。我们可以对其进行重采样并汇总,得到一个噪声更少的新序列,进而更易建模。
# Calculate average hourly fraction hourly_frac = train.groupby(['hour']).mean()/np.sum(train.groupby(['hour']).mean()) hourly_frac.drop(['ID'], axis = 1, inplace = True) hourly_frac.columns = ['fraction'] # convert to time series from dataframe train.index = train.Datetime train.drop(['ID','hour','Datetime'], axis = 1, inplace = True) daily_train = train.resample('D').sum()
Prophet要求时间序列中的变量名为:
y -> 目标(Target)
ds -> 时间(Datetime)
因此,下一步是基于上述规范来转换数据文件:
daily_train['ds'] = daily_train.index daily_train['y'] = daily_train.Count daily_train.drop(['Count'],axis = 1, inplace = True)
拟合Prophet模型(未设置的参数表明默认原设置):
#参数设置 m = Prophet(yearly_seasonality = True, seasonality_prior_scale=0.1) #模型拟合 m.fit(daily_train) #预测窗口 future = m.make_future_dataframe(periods=213) #模型预测 forecast = m.predict(future) #预测结果可视化 fig = m.plot(forecast)
我们可以通过以下命令来查看各个成分:
m.plot_components(forecast)

图 6 时间序列拆解
基于每日数据的预测如下。

图 7 时间序列预测
参考链接:https://www.analyticsvidhya.com/blog/ 2018/05/generate-accurate-forecasts-facebook-prophet-python-r/
案例:https://www.kaggle.com/arindamgot/eda-prophet-mlp-neural-network-forecasting
github: https://github.com/facebook/prophet/blob/master/python/fbprophet/forecaster.py

评论