Prophet 的算法实现¶
在时间序列分析领域，有一种常见的分析方法叫做时间序列的分解 (Decomposition of Time Series），它把时间序列 $y_{t}$ 分成几个部分，分别是季节项 $S_{t}$ ，趋势项 $T_{t}$ ，剩余项 $R_{t}$ 。也就是 说对所有的 $t \geq 0$ ，都有 $$ y_{t}=S_{t}+T_{t}+R_{t} $$ 除了加法的形式，还有乘法的形式，也就是: $$ y_{t}=S_{t} \times T_{t} \times R_{t} $$ 以上式子等价于 $\ln y_{t}=\ln S_{t}+\ln T_{t}+\ln R_{t}$ 。所以，有的时候在预测模型的时候，会先取 对数，然后再进行时间序列的分解，就能得到乘法的形式。在 fbprophet 算法中，作者们基于这 种方法进行了必要的改进和优化。 一般来说，在实际生活和生产环节中，除了季节项，趋势项，剩余项之外，通常还有节假日的效 应。所以，在 prophet 算法里面，作者同时考虑了以上四项，也就是: $$ y(t)=g(t)+s(t)+h(t)+\epsilon_{t} $$
其中 $g(t)$ 表示趋势项，它表示时间序列在非周期上面的变化趋势； $s(t)$ 表示周期项，或者称为 季节项，一般来说是以周或者年为单位； $h(t)$ 表示节假日项，表示在当天是否存在节假日； $\epsilon_{t}$ 表 示误差项或者称为剩余项。Prophet 算法就是通过拟合这几项，然后最后把它们累加起来就得到了 时间序列的预测值。
趋势项模型 $g(t)$¶
在 Prophet 算法里面，趋势项有两个重要的函数，一个是基于逻辑回归函数 (logistic function) 的，另一个是基于分段线性函数（piecewise linear function）的。 首先，我们来介绍一下基于逻辑回归的趋势项是怎么做的。 如果回顾逻辑回归函数的话，一般都会想起这样的形式： $\sigma(x)=1 /\left(1+e^{-x}\right)$, 它的导数是 $\sigma^{\prime}(x)=\sigma(x) \cdot(1-\sigma(x))$, 并且 $\lim _{x \rightarrow+\infty} \sigma(x)=1, \lim _{x \rightarrow-\infty} \sigma(x)=0$. 如果增加一些 参数的话，那么逻辑回归就可以改写成: $$ f(x)=C /\left(1+e^{-k(x-m)}\right) $$
这里的 $C$ 称为曲线的最大渐近值， $k$ 表示曲线的增长率， $m$ 表示曲线的中点。当 $C=1, k=1, m=0$ 时，恰好就是大家常见的 sigmoid 函数的形式。从 sigmoid 的函数表达 式来看，它满足以下的微分方程: $y^{\prime}=y(1-y)$ 。 那么，如果使用分离变量法来求解微分方程 $y^{\prime}=y(1-y)$ 就可以得到: $$ \frac{y^{\prime}}{y}+\frac{y^{\prime}}{1-y}=1 \Rightarrow \ln \frac{y}{1-y}=1 \Rightarrow y=1 /\left(1+K e^{-x}\right) $$ 但是在现实环境中，函数 $f(x)=C /\left(1+e^{-k(x-m)}\right)$ 的三个参数 $C, k, m$ 不可能都是常数， 而很有可能是随着时间的迁移而变化的，因此，在 Prophet 里面，作者考虑把这三个参数全部换 成了随着时间而变化的函数，也就是 $C=C(t), k=k(t), m=m(t)$. 除此之外，在现实的时间序列中，曲线的走势肯定不会一直保持不变，在某些特定的时候或者有着 某种潜在的周期曲线会发生变化，这种时候，就有学者会去研究变点检测，也就是所谓 change point detection。例如下面的这幅图的 $t_{1}^{*}, t_{2}^{*}$ 就是时间序列的两个变点。
在 Prophet 里面，是需要设置变点的位置的，而每一段的趋势和走势也是会根据变点的情况而改 变的。在程序里面有两种方法，一种是通过人工指定的方式指定变点的位置；另外一种是通过算法 来自动选择。在默认的函数里面，Prophet 会选择 n_changepoints $=25$ 个变点，然后设置变点 的范围是前 80\%（changepoint_range），也就是在时间序列的前 80\% 的区间内会设置变点。 通过 forecaster.py 里面的 set_changepoints 函数可以知道，首先要看一些边界条件是否合理， 例如时间序列的点数是否少于 nchangepoints 等内容；其次如果边界条件符合，那变点的位置 就是均匀分布的，这一点可以通过 np.linspace 这个函数看出来。 下面假设已经放置了 $S$ 个变点了，并且变点的位置是在时间戳 $s{j}, 1 \leq j \leq S$ 上，那么在这些时 间戳上，我们就需要给出增长率的变化，也就是在时间戳 $s{j}$ 上发生的 change in rate。可以假设 有这样一个向量: $\boldsymbol{\delta} \in \mathbb{R}^{S}$, 其中 $\delta{j}$ 表示在时间翟 $s{j}$ 上的增长率的变化量。如果一开始的增长率 我们使用 $k$ 来代替的话，那么在时间翟 $t$ 上的增长率就是 $k+\sum{j: t>s{j}} \delta{j}$ ，通过一个指示函数 $\mathbf{a}(t) \in{0,1}^{S}$ 就是 $$ a_{j}(t)=\left\{\begin{array}{l} 1, \text { if } t \geq s_{j} \\ 0, \text { otherwise } \end{array}\right. $$ 那么在时间翟 $t$ 上面的增长率就是 $k+\mathbf{a}^{T} \boldsymbol{\delta}$. 一旦变化量 $k$ 确定了，另外一个参数 $m$ 也要随之 确定。在这里需要把线段的边界处理好，因此通过数学计算可以得到: $$ \gamma_{j}=\left(s_{j}-m-\sum_{\ell<j} \gamma_{\ell}\right) \cdot\left(1-\frac{k+\sum_{\ell<j} \delta_{\ell}}{k+\sum_{a<j} \delta_{\ell}}\right) $$
所以，分段的逻辑回归增长模型就是: $$ g(t)=\frac{C(t)}{1+\exp \left(-\left(k+\boldsymbol{a}(t)^{t} \boldsymbol{\delta}\right) \cdot\left(t-\left(m+\boldsymbol{a}(t)^{T} \boldsymbol{\gamma}\right)\right.\right.}, $$ 其中， $$ \boldsymbol{a}(t)=\left(a_{1}(t), \cdots, a_{S}(t)\right)^{T}, \boldsymbol{\delta}=\left(\delta_{1}, \cdots, \delta_{S}\right)^{T}, \gamma=\left(\gamma_{1}, \cdots, \gamma_{S}\right)^{T} $$ 在逻辑回归函数里面，有一个参数是需要提前设置的，那就是 Capacity，也就是所谓的 $C(t)$ ， 在使用 Prophet 的 growth = 'logistic' 的时候，需要提前设置好 $C(t)$ 的取值才行。 再次，我们来介绍一下基于分段线性函数的趋势项是怎么做的。众所周知，线性函数指的是 $y=k x+b$, 而分段线性函数指的是在每一个子区间上，函数都是线性函数，但是在整段区间 上，函数并不完全是线性的。正如下图所示，分段线性函数就是一个折线的形状。
因此，基于分段线性函数的模型形如： $$ g(t)=(k+\boldsymbol{a}(t) \boldsymbol{\delta}) \cdot t+\left(m+\boldsymbol{a}(t)^{T} \gamma\right) $$ 其中 $k$ 表示增长率（growth rate）， $\boldsymbol{\delta}$ 表示增长率的变化量， $m$ 表示 offset parameter。而这 两种方法（分段线性函数与逻辑回归函数）最大的区别就是 $\gamma$ 的设置不一样，在分段线性函数 中， $\boldsymbol{\gamma}=\left(\gamma_{1}, \cdots, \gamma_{S}\right)^{T}, \gamma_{j}=-s_{j} \delta_{j}$. 注意：这与之前逻辑回归函数中的设置是不一样的。 在 Prophet 的源代码中， forecast.py 这个函数里面包含了最关键的步骤，其中 piecewise_logistic 函数表示了前面所说的基于逻辑回归的增长函数，它的输入包含了 cap 这个指 标，因此需要用户事先指定 capacity。而在 piecewise_linear 这个函数中，是不需要 capacity 这 个指标的，因此 $m=\operatorname{Prophet}()$ 这个函数默认的使用 growth = 'linear' 这个增长函数，也可以写 作 $m=$ Prophet(growth = 'linear')；如果想用 growth = 'logistic'，就要这样写:
from prophet import Prophet
model1 = Prophet(growth='logistic',yearly_seasonality=True, weekly_seasonality=True, daily_seasonality=True)
model1.add_country_holidays(country_name="CN")
data1['cap'] = 1000
变点的选择 (Changepoint Selection)¶
在介绍变点之前，先要介绍一下 Laplace 分布，它的概率密度函数为： $$ f(x \mid \mu, b)=\exp (-|x-\mu| / b) / 2 b $$ 其中 $\mu$ 表示位置参数， $b>0$ 表示尺度参数。Laplace 分布与正态分布有一定的差异。 在 Prophet 算法中，是需要给出变点的位置，个数，以及增长的变化率的。因此，有三个比较重 要的指标，那就是
- changepoint_range,
- n_changepoint,
- changepoint_prior_scale。 changepoint_range 指的是百分比，需要在前 changepoint_range 那么长的时间序列中设置变 点，在默认的函数中是 changepoint_range $=0.8$ 。 n_changepoint 表示变点的个数，在默认的 函数中是 n_changepoint $=25$ 。changepoint_priorscale 表示变点增长率的分布情况，在论文 中， $\delta{j} \sim$ Laplace $(0, \tau)$ ，这里的 $\tau$ 就是 change_point_scale。
在整个开源框架里面，在默认的场景下，变点的选择是基于时间序列的前 $80 \%$ 的历史数据，然后 通过等分的方法找到 25 个变点 (change points)，而变点的增长率是满足 Laplace 分布 $\delta_{j} \sim \operatorname{Laplace}(0,0.05)$ 的。因此，当 $\tau$ 趋近于零的时候， $\delta_{j}$ 也是趋向于零的，此时的增长函 数将变成全段的逻辑回归函数或者线性函数。这一点从 $g(t)$ 的定义可以轻易地看出。 对末来的预估 (Trend Forecast Uncertainty) 从历史上长度为 $T$ 的数据中，我们可以选择出 $S$ 个变点，它们所对应的增长率的变化量是 $\delta_{j} \sim \operatorname{Laplace}(0, \tau)$ 。此时我们需要预测末来，因此也需要设置相应的变点的位置，从代码中 看，在 forecaster.py 的 sample_predictive_trend 函数中，通过 Poisson 分布等概率分布方法找 到新增的 changepoint_tsnew 的位置，然后与 changepoint $t$ 拼接在一起就得到了整段序列的 changepoint_ts。 changepoint_ts_new $=1+\mathrm{np} \cdot$ random. rand(n_changes) $*(T-1)$ changepoint_ts $=n p$.concatenate ((self. changepoints_t, changepoint_ts_new $))$ 第一行代码的 1 保证了 changepoint_ts_new 里面的元素都大于 changets 里面的元素。除了变 点的位置之外，也需要考虑 $\delta$ 的情况。这里令 $\lambda=\sum{j=1}^{S}\left|\delta_{j}\right| / S$ ，于是新的增长率的变化量就 是按照下面的规则来选择的: 当 $j>T$ 时，
$$ \delta_{j}=\left\{\begin{array}{l} 0, \text { with probability }(T-S) / T \\ \sim \operatorname{Laplace}(0, \lambda), \text { with probability } S / T \end{array}\right. $$季节性趋势¶
几乎所有的时间序列预测模型都会考虑这个因素，因为时间序列通常会随着天，周，月，年等季节 性的变化而呈现季节性的变化，也称为周期性的变化。对于周期函数而言，大家能够马上联想到的 就是正弦余弦函数。而在数学分析中，区间内的周期性函数是可以通过正弦和余弦的函数来表示 的：假设 $f(x)$ 是以 $2 \pi$ 为周期的函数，那么它的傅立叶级数就是 $a_{0}+\sum_{n=1}^{\infty}\left(a_{n} \cos (n x)+b_{n} \sin (n x)\right)$ 。 在论文中，作者使用傅立叶级数来模拟时间序列的周期性。假设 $P$ 表示时间序列的周期， $P=365.25$ 表示以年为周期， $P=7$ 表示以周为周期。它的傅立叶级数的形式都是: $$ s(t)=\sum_{n=1}^{N}\left(a_{n} \cos \left(\frac{2 \pi n t}{P}\right)+b_{n} \sin \left(\frac{2 \pi n t}{P}\right)\right) $$ 就作者的经验而言，对于以年为周期的序列（ $P=365.25$ ) 而言， $N=10$ ；对于以周为周 期的序列 $(P=7)$ 而言， $N=3$ 。这里的参数可以形成列向量: $$ \boldsymbol{\beta}=\left(a_{1}, b_{1}, \cdots, a_{N}, b_{N}\right)^{T} . $$
当 $N=10$ 时， $$ X(t)=\left[\cos \left(\frac{2 \pi(1) t}{365.25}\right), \cdots, \sin \left(\frac{2 \pi(10) t}{365.25}\right)\right] $$ 当 $N=3$ 时， $$ X(t)=\left[\cos \left(\frac{2 \pi(1) t}{7}\right), \cdots, \sin \left(\frac{2 \pi(3) t}{7}\right)\right] $$ 因此，时间序列的季节项就是: $s(t)=X(t) \boldsymbol{\beta}$, 而 $\boldsymbol{\beta}$ 的初始化是 $\boldsymbol{\beta} \sim \operatorname{Normal}\left(0, \sigma^{2}\right)$ 。这里 的 $\sigma$ 是通过 seasonality_prior_scale 来控制的，也就是说 $\sigma=$ seasonality_prior_scale。这个 值越大，表示季节的效应越明显；这个值越小，表示季节的效应越不明显。同时，在代码里面， seasonality_mode 也对应着两种模式，分别是加法和乘法，默认是加法的形式。在开源代码中， $X(t)$ 函数是通过 fourier_series 来构建的。
节假日效应 (holidays and events) 在现实环境中，除了周末，同样有很多节假日，而且不同的国家有着不同的假期。在 Prophet 里 面，通过维基百科里面对各个国家的节假日的描述，hdays.py 收集了各个国家的特殊节假日。除 了节假日之外，用户还可以根据自身的情况来设置必要的假期，例如 The Super Bowl，双十一 等。
| Holiday | Country | Year | Date |
|---|---|---|---|
| Thanksgiving | US | 2015 | 26 Nov 2015 |
| Thanksgiving | US | 2016 | 24 Nov 2016 |
| Thanksgiving | US | 2017 | 23 Nov 2017 |
| Thanksgiving | US | 2018 | 22 Nov 2018 |
| Christmas | ** | 2015 | 25 Dec 2015 |
| Christmas | ** | 2016 | 25 Dec 2016 |
| Christmas | ** | 2017 | 25 Dec 2017 |
| Christmas | ** | 2018 | 25 Dec 2018 |
由于每个节假日对时间序列的影响程度不一样，例如春节，国庆节则是七天的假期，对于劳动节等 假期来说则假日较短。因此，不同的节假日可以看成相互独立的模型，并且可以为不同的节假日设 置不同的前后窗口值，表示该节假日会影响前后一段时间的时间序列。用数学语言来说，对与第 $i$ 个节假日来说， $D_{i}$ 表示该节假日的前后一段时间。为了表示节假日效应，我们需要一个相应的指 示函数（indicator function），同时需要一个参数 $\kappa_{i}$ 来表示节假日的影响范围。假设我们有 $L$ 个节假日，那么 $h(t)=Z(t) \kappa=\sum_{i=1}^{L} \kappa_{i} \cdot 1_{\left\{t \in D_{i}\right\}}$ 其中 $Z(t)=\left(1_{\left\{t \in D_{1}\right\}}, \cdots, 1_{\left\{t \in D_{L}\right\}}\right), \kappa=\left(\kappa_{1}, \cdots, \kappa_{L}\right)^{T}$. 其中 $\boldsymbol{\kappa} \sim \operatorname{Normal}\left(0, v^{2}\right)$ 并且该正态分布是受到 $v=$ holidays_prior_scale 这个指标影响 的。默认值是 10，当值越大时，表示节假日对模型的影响越大；当值越小时，表示节假日对模型 的效果越小。用户可以根据自己的情况自行调整。
模型拟合 (Model Fitting) 按照以上的解释，我们的时间序列已经可以通过增长项，季节项，节假日项来构建了，i.e. $$ y(t)=g(t)+s(t)+h(t)+\epsilon $$ 下一步我们只需要拟合函数就可以了，在 Prophet 里面，作者使用了 pyStan 这个开源工具中的 L-BFGS 方法来进行函数的拟合。具体可以参考 forecast.py 里面的 stan_init 函数。 Prophet 中可以设置的参数 在 Prophet 中，用户一般可以设置以下四种参数：
- Capacity：在增量函数是逻辑回归函数的时候，需要设置的容量值。
- Change Points：可以通过 n_changepoints 和 changepoint_range 来进行等距的变点设置， 也可以通过人工设置的方式来指定时间序列的变点。
- 季节性和节假日：可以根据实际的业务需求来指定相应的节假日。
- 光滑参数： $\tau=$ changepoint_prior_scale 可以用来控制趋势的灵活度， $\sigma=$ seasonality_prior_scale 用来控制季节项的灵活度， $v=$ holidays prior scale 用来控制节假 日的灵活度。
Prophet 的实际使用¶
Prophet 的简单使用¶
因为 Prophet 所需要的两列名称是 'ds' 和 ' $y$ '，其中， 'ds' 表示时间戳， ' $y$ ' 表示时间序列的值， 因此通常来说都需要修改 pd.dataframe 的列名字。如果原来的两列名字是 'timestamp' 和 'value' 的话，只需要这样写:
df = df.rename(columns={'timestamp':'ds', 'value':'y'})
如果 'timestamp' 是使用 unixtime 来记录的，需要修改成 YYYY-MM-DD hh:mm:ss 的形式:
df['ds'] = pd.to_datetime(df['ds'],unit='s')
在一般情况下，时间序列需要进行归一化的操作，而 pd.dataframe 的归一化操作也十分简单:
df['y'] = (df['y'] - df['y'].mean()) / (df['y'].std())
然后就可以初始化模型，然后拟合模型，并且进行时间序列的预测了。
#初始化模型：
m = Prophet()
# 拟合模型：
m.fit(df)
# 计算预测值：periods 表示需要预测的点数，freq 表示时间序列的频率。
future = m.make_future_dataframe(periods=30, freq='min')
future.tail()
forecast = m.predict(future)
在进行了预测操作之后，通常都希望把时间序列的预测趋势画出来：
#画出预测图：
m.plot(forecast)
#画出时间序列的分量：
m.plot_components(forecast)
<img width="60%" src="images/pre5905.png">
Prophet 的默认参数
```python
def __init__(
self,
growth='linear',
changepoints=None,
n_changepoints=25,
changepoint_range=0.8,
yearly_seasonality='auto',
weekly_seasonality='auto',
daily_seasonality='auto',
holidays=None,
seasonality_mode='additive',
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,
):
结束语¶
对于商业分析等领域的时间序列，Prophet 可以进行很好的拟合和预测，但是对于一些周期性或者趋势性不是很强的时间序列，用 Prophet 可能就不合适了。但是，Prophet 提供了一种时序预测的方法，在用户不是很懂时间序列的前提下都可以使用这个工具得到一个能接受的结果。具体是否用 Prophet 则需要根据具体的时间序列来确定
资料来源
