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
评论