初识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 = 1 n ∑ i = 1 n ∣ y ^ i − y i ∣ M A E=\frac{1}{n} \sum_{i=1}^{n}\left|\hat{y}_{i}-y_{i}\right| MAE=n1∑i=1n∣y^i−yi∣
范围[0,+∞),当预测值与真实值完全吻合时等于0,即完美模型;误差越大,该值越大。
-
MAPE:平均绝对百分比误差, M A P E = 100 % n ∑ i = 1 n ∣ y ^ i − y i y i ∣ M A P E=\frac{100 \%}{n} \sum_{i=1}^{n}\left|\frac{\hat{y}_{i}-y_{i}}{y_{i}}\right| MAPE=n100%∑i=1n∣∣∣yiy^i−yi∣∣∣
范围[0,+∞),MAPE 为0%表示完美模型,MAPE 大于 100 %则表示劣质模型。
- MAE:平均绝对误差,
-
定义计算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()
标签:plot,upper,--,Prophet,forecast,df,初识,yhat,ds 来源: https://blog.csdn.net/JanetHULAHA/article/details/112062700