初识Prophet模型（二）-- 应用篇

相关学习：初识Prophet模型（一）--理论篇

7、Prophet 模型应用

7.0 背景描述

该案例使用的是wiki网站日访问量（数值经过log处理）的csv数据文件
描述的是美国著名橄榄球四分卫的维基页面浏览量，他是美国球员，一年里的周期规律会起很大作用，而一周里的周期规律也很明显。

7.1 导入数据

import pandas as pd
import numpy as np
from fbprophet import Prophet
import matplotlib.pyplot as plt
%matplotlib inline

df = pd.read_csv('data.csv')
df.head()

df.dtypes #检查下df的数据类型
df['ds'] = df['ds'].apply(pd.to_datetime)# ds列必须是pandas的datetime数据类型，使用pandas自带的pd.to_datetime将日期转为datetime类型

plt.rcParams['figure.figsize']=(20,10)
plt.style.use('ggplot')
df.set_index('ds').y.plot()

7.2 拟合模型

model = Prophet(daily_seasonality=True)
model.fit(df)

<fbprophet.forecaster.Prophet at 0x10715a0b8>

7.3 预测（使用默认参数）

生成一个未来的日期的dataframe，然后用训练好的模型prophet来predict。

future = model.make_future_dataframe(periods=730)
future.tail()

有了未来的日期,就可以使用学习到的趋势来预测未来日期的走势。
预测的结果包括如下变量

'ds', 'trend', 'yhat_lower', 'yhat_upper', 'trend_lower', 'trend_upper',
       'additive_terms', 'additive_terms_lower', 'additive_terms_upper',
       'weekly', 'weekly_lower', 'weekly_upper', 'yearly', 'yearly_lower',
       'yearly_upper', 'multiplicative_terms', 'multiplicative_terms_lower',
       'multiplicative_terms_upper', 'yhat'

我们只用 'ds', 'yhat', 'yhat_lower', 'yhat_upper'

forecast=model.predict(future)

forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].head()

print(fig1)

成分分析**

趋势是由不同的成分组成，比如总趋势、年、季节、月、周等等，我们要将这些成分从趋势中抽取出来看看不同成分的趋势情况

预测的结果包括如下变量

'ds', 'trend', 'yhat_lower', 'yhat_upper', 'trend_lower', 'trend_upper',
       'additive_terms', 'additive_terms_lower', 'additive_terms_upper',
       'weekly', 'weekly_lower', 'weekly_upper', 'yearly', 'yearly_lower',
       'yearly_upper', 'multiplicative_terms', 'multiplicative_terms_lower',
       'multiplicative_terms_upper', 'yhat'

下面图1是根据trend画出来的，图2是根据weekly画出来的，图3是根据yearly画出来的，图4是根据daily画出来的
因为是加法模型，所以：
- forecast['additive_terms'] = forecast['daily']+forecast['weekly'] + forecast['yearly']；
- forecast['yhat'] = forecast['trend'] + forecast['additive_terms']
- forecast['yhat'] = forecast['trend'] +forecast['daily']+ forecast['weekly'] + forecast['yearly']。
如果有节假日因素，那么就会有：
- forecast['yhat'] = forecast['trend'] +forecast['daily']+forecast['weekly'] + forecast['yearly'] + forecast['holidays']。
对于那些是节假日的天数， forecast['holidays']才会有值
不是节假日的天数，forecast['holidays']为0
因为是加法模型，'multiplicative_terms', 'multiplicative_terms_lower', 'multiplicative_terms_upper'这3列为空。

因此，在下面的拆解图中，weekly中的Monday为0.3的意思就是，在trend的基础上，加0.3；Saturday为-0.3的意思就是，在trend的基础上，减0.3。因此，这条线的高低也在一定程度上反应了“销量的趋势“。

fig2 = model.plot_components(forecast)
print(fig2)

第一幅趋势图里，可以看到按页面浏览量的总趋势。是因为作者最近退休，所以浏览量逐渐下降。
第二幅按周的周期规律图里能看出，在比赛当天和赛后（周日和周一）访问量明显较高。
美国NFL橄榄球比赛主要集中在九月到次年1月初，和二月初的超级碗。这也反映在第三幅图中，按年的周期规律图。

forecast_df = forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']]
df = pd.merge(df, forecast_df, on='ds', how='right')
df.set_index('ds').plot(figsize=(16,8), color=['royalblue', "green", "pink", "yellow"], grid=True);

x1 = forecast['ds']
y1 = forecast['yhat']
y2 = forecast['yhat_lower']
y3 = forecast['yhat_upper']
plt.plot(x1,y1)
plt.plot(x1,y2)
plt.plot(x1,y3)
plt.show()

print(forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].tail())

7.4 趋势突变点

自动检测变化点

步骤：
- 1、通过大量的速率改变的点检测变化点
- 2、对这些点做稀疏先验

默认情况下会检测出25个变化点，这些点均匀的分布在前80%的时间序列中

fig = model.plot(forecast)
for cp in model.changepoints:
    plt.axvline(cp, c='pink', ls='--', lw=2)

因为稀疏先验，大部分的变化点并没有用到，看看每个变化点的速率变化图：

deltas = model.params['delta'].mean(0)
fig = plt.figure(facecolor='w', figsize=(10, 6))
ax = fig.add_subplot(111)
ax.bar(range(len(deltas)), deltas, facecolor='#0072B2', edgecolor='#0072B2')
ax.grid(True, which='major', c='gray', ls='-', lw=1, alpha=0.2)
ax.set_ylabel('Rate change')
ax.set_xlabel('Potential changepoint')
fig.tight_layout()

变化点的数量可以通过参数n_changepoints指定，但最好还是通过调整正则化来修改
看看比较明显的变化点

from fbprophet.plot import add_changepoints_to_plot
fig = model.plot(forecast)
a = add_changepoints_to_plot(fig.gca(), model, forecast) #虚线处为给定时间序列中的变点

调整趋势灵活性

当趋势出现过拟合或者欠拟合的情况下，可以通过参数changepoint_prior_scale调整稀疏先验的程度，默认为0.05
该参数值越大，则趋势越灵活

增大灵活性

m = Prophet(changepoint_prior_scale=0.5)
forecast = m.fit(df).predict(future)
fig = m.plot(forecast)

减少灵活性

forecast = m.fit(df).predict(future)
fig = m.plot(forecast)

指定变化点的位置

通过参数changepoints手动指定位置，只有指定的这些点可以有速率变化

m = Prophet(changepoints=['2014-01-01'])
forecast = m.fit(df).predict(future)
fig = m.plot(forecast)

7.5季节性、假期效应和回归因子

假期和特殊事件建模

要给假期或者其他重复事件建模，就需要创建一个包含holiday和ds列的DataFrame
需要包含过去和将来所有的特殊日子，如果这些特殊日子没有出现在将来(要预测的日期)，那预测就不会用到
可以在这个数据框基础上再新建两列 lower_window 和 upper_window ，从而将节假日的时间扩展成一个区间 [ lower_window , upper_window ] ，例如：
- 如果想将平安夜也加入到 “圣诞节” 里，就设置 lower_window = -1 , upper_window = 0 ;
- 如果想将黑色星期五加入到 “感恩节” 里，就设置 lower_window = 0 , upper_window = 1

playoffs = pd.DataFrame({
  'holiday': 'playoff',
  'ds': pd.to_datetime(['2008-01-13', '2009-01-03', '2010-01-16',
                        '2010-01-24', '2010-02-07', '2011-01-08',
                        '2013-01-12', '2014-01-12', '2014-01-19',
                        '2014-02-02', '2015-01-11', '2016-01-17',
                        '2016-01-24', '2016-02-07']),
  'lower_window': 0,
  'upper_window': 1,
})
superbowls = pd.DataFrame({
  'holiday': 'superbowl',
  'ds': pd.to_datetime(['2010-02-07', '2014-02-02', '2016-02-07']),
  'lower_window': 0,
  'upper_window': 1,
})
holidays = pd.concat((playoffs, superbowls))

上面superbowl的日期也包含在playoff的日期中，也就是superbowl日期的影响会有个叠加效应

m = Prophet(holidays=holidays)
forecast = m.fit(df).predict(future)

可以通过forecast看看假期效应

forecast[(forecast['playoff'] + forecast['superbowl']).abs() > 0][['ds', 'playoff', 'superbowl']][-15:]

看看假期效应在图上的显示，playoff日期有高峰，superbowl日期有更明显的高峰

fig = m.plot_components(forecast)

对假期单独画图

from fbprophet.plot import plot_forecast_component
plot_forecast_component(m, forecast, 'superbowl')

内置假期

可以通过add_country_holidays使用内置假期
通过模型的train_holiday_names方法查看哪些假期

m = Prophet(holidays=holidays)
m.add_country_holidays(country_name='CN')
m.fit(df)

m.train_holiday_names

INFO:fbprophet:Disabling daily seasonality. Run prophet with daily_seasonality=True to override this.

0                 playoff
1               superbowl
2          New Year's Day
3        Chinese New Year
4       Tomb-Sweeping Day
5               Labor Day
6    Dragon Boat Festival
7     Mid-Autumn Festival
8            National Day
dtype: object

m = Prophet(holidays=holidays)
m.add_country_holidays(country_name='US')
m.fit(df)

forecast = m.predict(future)
fig = m.plot_components(forecast)

季节性的傅里叶级数

年度季节性的傅里叶级数默认是10

from fbprophet.plot import plot_yearly
m = Prophet().fit(df)
a = plot_yearly(m)

默认值大多数是没问题的，但是季节性可能有更高频率的变化，并且通常没有这么平滑，这时可以增加这个值
增加这个值也可能导致过拟合,这里增加到20

from fbprophet.plot import plot_yearly
m = Prophet(yearly_seasonality=20).fit(df)
a = plot_yearly(m)

自定义季节性

时间序列超过两个周期时，Prophet默认训练星期和年的季节性
在sub-daily的时间序列时，会训练每天的季节性
可以使用函数add_seasonality添加小时/月/季度等其他季节性
函数add_seasonality的参数：
- name 哪种周期
- period 季节性的周期
- fourier_order 季节性的傅里叶级数
- prior_scale 可选参数
- 默认情况下，周的季节性傅里叶级数为3，年的季节性傅里叶级数为10
这里将每周的季节性替换为每月的季节性(period=30.5)

m = Prophet(weekly_seasonality=False)
m.add_seasonality(name='monthly', period=30.5, fourier_order=5)
forecast = m.fit(df).predict(future)
fig = m.plot_components(forecast)

依赖于其他因素的季节性

有时候季节性依赖于其他一些因素，比如每周的季节性在夏天和其他季节表现不一致，每天的季节性在周末和周内表现不一致，这种季节性可以使用带条件的季节性训练模型
在一年中默认每周的季节性表现是一致的，但是可能希望每周的季节性在淡季和旺季(每周末有比赛)表现不一致
可以使用带条件的季节性为淡季和旺季单独构建每周的季节性

这里先增加一列布尔类型的数据，来表示日期在淡季还是旺季：

def is_nfl_season(ds):
    date = pd.to_datetime(ds)
    return (date.month > 8 or date.month < 2)

df['on_season'] = df['ds'].apply(is_nfl_season)
df['off_season'] = ~df['ds'].apply(is_nfl_season)

接着禁用内置的每周季节性，使用淡季的周季节性和旺季的周季节性代替
因此，只有condition_name列为True的时候季节性才有日期
在预测的DataFrame上，也要做同样的操作

m = Prophet(weekly_seasonality=False)
m.add_seasonality(name='weekly_on_season', period=7, fourier_order=3, condition_name='on_season')
m.add_seasonality(name='weekly_off_season', period=7, fourier_order=3, condition_name='off_season')

future['on_season'] = future['ds'].apply(is_nfl_season)
future['off_season'] = ~future['ds'].apply(is_nfl_season)
forecast = m.fit(df).predict(future)
fig = m.plot_components(forecast)

从图中可以看到，在旺季的时候每周末都会打球，周日和周一都有大幅度增长，但在淡季则完全没有。

假期和季节性的prior scale

如果发现假期过拟合，可以设置参数holidays_prior_scale调整假期的prior scale使之平滑
这个参数默认是10，减少可以限制假期效果

m = Prophet(holidays=holidays, holidays_prior_scale=0.05).fit(df)
forecast = m.predict(future)
forecast[(forecast['playoff'] + forecast['superbowl']).abs() > 0][['ds', 'playoff', 'superbowl']][-10:]

可以看到，比起之前假期效应被减弱了，特别是在观看最少的superbowls上

可以用下面的方式设置每周季节性的prior_scale

额外的回归特征

可以使用函数add_regressor将其他回归特征添加到模型的线性部分
训练和预测的数据集上都需要包含这些回归特征的值

下面，为NFL赛季的每周日添加这样一个回归特征，再画图看看这个特征的效果

def nfl_sunday(ds):
    date = pd.to_datetime(ds)
    if date.weekday() == 6 and (date.month > 8 or date.month < 2):
        return 1
    else:
        return 0
df['nfl_sunday'] = df['ds'].apply(nfl_sunday)

m = Prophet()
m.add_regressor('nfl_sunday')
m.fit(df)

future['nfl_sunday'] = future['ds'].apply(nfl_sunday)

forecast = m.predict(future)
fig = m.plot_components(forecast)

也可以使用前面说过的holidays的接口，通过创建一个过去和未来的这些周日的list来处理NFL赛季周日的这种情况
函数add_regressor为定义额外的线性回归提供了一个更加通用的接口

7.6 模型诊断（内置方法）

Prophet有交叉验证功能，具体做法是通过在历史数据中选择一些截断点，对于这些截断点，只使用这些点之前的数据来拟合模型，然后比较真实值和预测值

下面模型使用前五年的数据训练，预测后一年的数据

m = Prophet()
m.fit(df)
future = m.make_future_dataframe(periods=366)

from fbprophet.diagnostics import cross_validation

df_cv = cross_validation(m, '365 days', initial='1825 days', period='365 days')
cutoff = df_cv['cutoff'].unique()[0]
df_cv = df_cv[df_cv['cutoff'].values == cutoff]

fig = plt.figure(facecolor='w', figsize=(10, 6))
ax = fig.add_subplot(111)
ax.plot(m.history['ds'].values, m.history['y'], 'k.')
ax.plot(df_cv['ds'].values, df_cv['yhat'], ls='-', c='#0072B2')
ax.fill_between(df_cv['ds'].values, df_cv['yhat_lower'],
                df_cv['yhat_upper'], color='#0072B2',
                alpha=0.2)
ax.axvline(x=pd.to_datetime(cutoff), c='gray', lw=4, alpha=0.5)
ax.set_ylabel('y')
ax.set_xlabel('ds')
ax.text(x=pd.to_datetime('2010-01-01'),y=12, s='Initial', color='black',
       fontsize=16, fontweight='bold', alpha=0.8)
ax.text(x=pd.to_datetime('2012-08-01'),y=12, s='Cutoff', color='black',
       fontsize=16, fontweight='bold', alpha=0.8)
ax.axvline(x=pd.to_datetime(cutoff) + pd.Timedelta('365 days'), c='gray', lw=4,
           alpha=0.5, ls='--')
ax.text(x=pd.to_datetime('2013-01-01'),y=6, s='Horizon', color='black',
       fontsize=16, fontweight='bold', alpha=0.8);

可以使用函数cross_validation给这些历史截断点自动完成交叉验证，参数如下：
- horizon：代表每次从cutoff往后预测多少天
- initial ：一开始的时间是多少
- period ：代表每隔多长时间设置一个cutoff
- 默认情况下，period是horizon的三倍，并且每隔半个horizon设置一个截断点
交叉验证的输出是一个DataFrame，包含真实的y和预测的yhat，可以用来评判效果

下面的交叉验证，horizon=365天，initial=730天，period=180天，在八年的时间序列中，等于有11((365*8-730-365)/180)个总的预测

from fbprophet.diagnostics import cross_validation
df_cv = cross_validation(m, initial='730 days', period='180 days', horizon = '365 days')
df_cv.head()

函数performance_metrics可以用来评判模型效果，计算得到：
- 统计信息包括均方误差(mean squared error, MSE)
- 均方根误差(root mean squared error, RMSE)
- 平均绝对误差(mean absolute error, MAE)
- 平均绝对误差(mean absolute percent error, MAPE)
- yhat_lower和yhat_upper估计的覆盖率

from fbprophet.diagnostics import performance_metrics
df_p = performance_metrics(df_cv)
df_p.head()

可以使用plot_cross_validation_metric来可视化这些指标，查看mape指标下的可视化，可以看到对未来一个月的预测有5%的误差，一年后的预测误差增加到了11%

from fbprophet.plot import plot_cross_validation_metric
fig = plot_cross_validation_metric(df_cv, metric='mape')

7.7 模型评估

通过历史数据对已知的最后一年数据进行预测并评估模型

prediction_size = 365
train_df = df[:-prediction_size]
train_df.tail()

model2 = Prophet(daily_seasonality=True)
model2.fit(train_df)
future2 = model2.make_future_dataframe(periods=365)
forecast2 = model2.predict(future2)

model2 .plot(forecast2);

定义一个辅助函数，用于从原数据集df中获取实际值y，然后和forecast对象中的预测值比较

def make_comparison_dataframe(historical, forecast):
    return forecast.set_index('ds')[['yhat', 'yhat_lower', 'yhat_upper']].join(historical.set_index('ds'))

cmp_df = make_comparison_dataframe(df, forecast2)
cmp_df.tail()

预测评价指标：
- MAE:平均绝对误差,
  
  $M A E=\frac{1}{n} \sum_{i=1}^{n}\left|\hat{y}_{i}-y_{i}\right|$
范围[0,+∞)，当预测值与真实值完全吻合时等于0，即完美模型；误差越大，该值越大。
- MAPE：平均绝对百分比误差,
  $M A P E=\frac{100 \%}{n} \sum_{i=1}^{n}\left|\frac{\hat{y}_{i}-y_{i}}{y_{i}}\right|$
  
  范围[0,+∞)，MAPE 为0%表示完美模型，MAPE 大于 100 %则表示劣质模型。
定义计算MAPE和MAE的辅助函数

def calculate_forecast_errors(df, prediction_size):
    df = df.copy()
    df['e'] = df['y'] - df['yhat']
    df['p'] = 1* df['e'] / df['y']
    predicted_part = df[-prediction_size:]
    error_mean = lambda error_name: np.mean(np.abs(predicted_part[error_name]))
    return {'MAPE': error_mean('p'), 'MAE': error_mean('e')}

计算MAPE和MAE

for err_name, err_value in calculate_forecast_errors(cmp_df, prediction_size).items():
    print(err_name, err_value)

MAPE 0.053184142465032766
MAE 0.4132303661978998

Box-Cox变换

定义计算函数

def inverse_boxcox(y, lambda_):
    return np.exp(y) if lambda_ == 0 else np.exp(np.log(lambda_ * y + 1) / lambda_)

准备数据，设置索引

train_df2 = train_df.copy().set_index('ds')

应用Box-Cox变换。这里它将返回两个值，第一个值是转换后的序列，第二个值是找到的最优λ值（最大似然）

from scipy import stats
import statsmodels.api as sm
train_df2['y'], lambda_prophet = stats.boxcox(train_df2['y'])
train_df2.reset_index(inplace=True)

创建一个新Prophet模型，并重复之前的拟合-预测流程

model3  = Prophet(daily_seasonality=True)
model3 .fit(train_df2)
future3 = model3.make_future_dataframe(periods=prediction_size)
forecast3 = model3.predict(future3)

通过逆函数和已知的λ值反转Box-Cox变换

for column in ['yhat']:
    forecast3[column] = inverse_boxcox(forecast3[column],lambda_prophet)

计算MAPE和MAE

cmp_df2 = make_comparison_dataframe(df, forecast3)

for err_name, err_value in calculate_forecast_errors(cmp_df2, prediction_size).items():
     print(err_name, err_value)

MAPE 0.04373071028220759
MAE 0.34260353853143777

对最后一年的真实值与预测值进行可视化对比

test_df = df[-prediction_size:]

test_df = test_df.set_index('ds')

forecast2 = forecast2[['ds','yhat']].set_index('ds')

df_all = forecast2.join(test_df).dropna()

df_all.head()

df_all.plot()
plt.rcParams['figure.figsize']=(30,20)
plt.style.use('ggplot')
plt.legend(['true', 'yhat'])
plt.show()