時序預測框架 — Prophet介紹

網易嚴選有很多業務決策依賴預測模型的輸入，其中時間序列預測是一類比較基礎的模型，服務於採購、營銷、倉配、客服等業務。

在客服話務量預測和單量預測場景中我們用到了開源時序預測框架Prophet。本篇介紹Prophet的基本原理和使用方法。

1. 簡介

2. 安裝

3. 基本用法

4. 算法原理

4.1 趨勢模型

4.2 季節模型

4.3 節假日模型

4.4 Prophet模型

5. 源碼

5.1 模型初始化

5.2 模型訓練

5.3 生成測試集

5.4 預測

6. 總結

7. 參考

1. 簡介

Prophet[1]是 Facebook 在2017年開源的時間序列預測算法，提供了 R 和 Python 語言的支持[2]。

Prophet容易上手，短短幾行代碼就能建立時序預測模型，算法的基本思想類似於時間序列分解，將時間序列分成為趨勢（Trend），季節性（seasonality）和節日（holiday）。

Prophet能自動識別一些簡單周期（年/周/日），並內置了節假日模塊，甚至支持了中國節日（包括農曆節日，如端午節、中秋節）。此外，它還能自動剔除異常值和處理缺失值。

Prophet提供了大量可配置的參數，使用者可根據具體需求調整模型，比如調整季節性的擬合度，添加自定義節假日，添加自定義變量等。不同的參數配置可能得到完全不一樣的模型，所以使用Prophet時需要對自定義參數有一些了解。

2. 安裝

通過Anaconda安裝Prophet非常簡單（用pip安裝可能會踩坑）：

# 1. 創建虛擬環境
conda create -n prophet python=3.7
# 2. 激活虛擬環境
source activate prophet
# 3. 通過conda-forge安裝prophet
conda install -c conda-forge fbprophet
3. 基本用法
Prophet的基本用法分以下幾步：
# 導入包
import pandas as pd
from fbprophet import Prophet

# 讀數據
# df(dataframe)需要包含兩列，一列 date，一列 y
df = pd.DataFrame({
    'ds': pd.date_range('20071210',periods=3700),
    'y': [np.random.rand() + np.cos(0.05 * x) for x in range(3700)]
})
# 初始化 prophet 模型
m = Prophet()
# 訓練 prophet 模型
m.fit(df)
# 生成預測數據（默認天粒度）
future = m.make_future_dataframe(periods=365)
# 使用訓練好的模型進行預測
forecast = m.predict(future)
# 預測繪圖
fig1 = m.plot(forecast)
# 分量繪圖
fig2 = m.plot_components(forecast)
數據示例
輸入和輸出數據格式如下：
輸入數據 df: pd.DataFrame
-dsy02007-12-109.59076112007-12-118.51959022007-12-128.18367732007-12-138.07246742007-12-147.893572
預測結果  forecast: pd.DataFrame

dsyhatyhat_loweryhat_upper32652017-01-158.2129427.4635608.93721532662017-01-168.5379937.7902599.26749232672017-01-178.3254287.5256759.05939132682017-01-188.1580597.4336348.88362732692017-01-198.1700467.4318018.840703
更多的參數配置請參考官方示例[3]。
4. 算法原理
Prophet基於時間序列的分解，把時間序列分解成四個部分，分別是趨勢項，季節項，節日項和剩餘項，將這幾項之和作為預測值：
其中
Prophet 算法就是通過擬合這幾項，然後最後把它們累加起來就得到了時間序列的預測值。
4.1 趨勢模型
在 Prophet 中，趨勢項可以有兩種選擇，一個是飽和增長模型（saturating model），也稱邏輯回歸模型（logistic function）；另一個是分段線性模型（piecewise linear model）。下面分別介紹這兩個模型：
飽和增長模型
在自然界中，很多趨勢是非線性增長的，比如人口，經濟增速等。它們通常都有一定的承載上界，比如facebook的用戶數量，這種增長常用 Logistic 增長模型進行建模。
Logistic 函數的基本形式是：
函數圖像如下圖所示：
在實際情況中，擬合趨勢
其中
當
漸進值
增長率變點（changepoint）的概念，變點前後時序的趨勢發生了變化。
Prophet默認會設置25個變點（n_changepoints=25），設置範圍是前80%的時序數據（changepoint_range=0.8），也就是在時間序列的前 80% 的區間內會設置變點，而且是均勻分布，如下圖所示。
假設在時間序列中設置
其中，
增長率
如圖所示，縱軸表示時序趨勢值，橫軸表示時間紫線所示，增長率紅線所示。但是在變點處，也就是
如下圖所示，當偏移量
在變點
最後logistic model表示為：
注意：在選擇邏輯回歸模型作為趨勢項時，需要人為指定
分段線性模型
對於那些不顯式為飽和增長的預測問題，Prophet還提供了另外一種簡潔且有用的趨勢模型 — 分段線性模型。變點之間的趨勢是線性函數，如下圖所示。
線性函數指的是
和邏輯回歸模型相同，分段線性模型同樣加入了變點，分段線性模型可表示為：
其中，
4.2 季節模型
通常時間序列包含多周期季節性（multi-period seasonality），比如寒暑假的影響以年為周期，工作日的影響以周為周期等。
Prophet採用傅立葉級數來擬合周期效應。傅立葉級數能將任何周期函數或周期信號分解成一個（可能由無窮個元素組成的）簡單振蕩函數的集合，即正弦函數和餘弦函數。
舉一個例子，以
可視化效果如下：

假設一個以
Prophet內置了三種周期：年、周和日。當
Prophet 作者按照經驗，設置了這些周期的默認級數，對於以年為周期的序列，
使用傅立葉級數擬合周期性需要估計
對於年周期（
時間
4.3 節假日模型
節假日和事件往往對時間序列預測有較大的影響，且通常不遵循周期模式，所以無法用季節性模型去建模節日。另外，節日和時間也有可能是非穩定周期的，比如中國的農曆節日，NBA的決賽時間，手機的新品發布時間等等。
不同的節假日可以看成相互獨立的模型，且不同的節假日有不同的前後窗口值，表示該節假日會影響前後一段時間的時間序列。比如春節有15天的窗口，而清明節只有3天窗口期。
對於第
用參數
用指示函數
時間
其中
4.4 Prophet模型
綜上，Prophet的預測模型[4]可以表示為：
若選擇邏輯回歸模型作為趨勢模型，上式可演變為：
若趨勢項為分段線性模型：
訓練時，採用L-BFGS做最大似然估計（maximum likelihood estimate）訓練出模型中的參數。
5. 源碼
若想完全掌握Prophet模型，還是需要去擼一遍源碼，一方面是在使用時心裡有個底，另一方面是為了更好的調參。
通過 debug 以下代碼走一遍 Prophet 的流程。
# 初始化 prophet 模型
m = Prophet()
# 訓練 prophet 模型
m.fit(df)
# 生成預測數據（默認天粒度）
future = m.make_future_dataframe(periods=365)
# 使用訓練好的模型進行預測
forecast = m.predict(future)
5.1 模型初始化
首先，初始化prophet模型：
m = Prophet()
初始化的時候可以指定一些參數，若不指定，Prophet會指定參數的默認值。具體的參數默認初始化如下：
# func __init__()
# 初始化參數（可設置）
growth='linear',              # 趨勢選擇 'linear'或'logistic'（分段線性 或 邏輯回歸）
changepoints=None,            # 變點時間（Series）
n_changepoints=25,            # changepoint個數:默認25
changepoint_range=0.8,        # changepoint範圍:默認前80%
yearly_seasonality='auto',    # 年周期性，可設 True or False
weekly_seasonality='auto',    # 周周期性，可設 True or False
daily_seasonality='auto',     # 日周期性，可設 True or False
holidays=None,                # 節假日，默認為None
seasonality_mode='additive',  # 周期疊加方式，'additive'或'multiplicative'
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,     # 

# 其他參數
start = None                   # 訓練集開始時間
y_scale = None                 # 訓練集 y 最大值
logistic_floor = False
t_scale = None                 # 訓練集時間跨度
changepoints_t = None          # 變點在時間軸的比例
seasonalities = OrderedDict({})
extra_regressors = OrderedDict({})
country_holidays = None       # 指定國家節日，中國為'CN'
stan_fit = None               # pystan擬合時的參數
params = {}                   # 訓練得到的參數
history = None                # 輸入數據
history_dates = None          # 歷史數據中的時間
train_component_cols = None   # 每列特徵對應的類別
component_modes = None        # {'additive':[],'multiplicative':[]}
train_holiday_names = None    # 配置的 holiday 名稱
其他參數都比較好理解，重點說一下後綴帶有prior_scale的參數。它們是控制變量的擬合程度的參數，可以理解為正則項係數，值越高，越趨向過擬合，值越低越趨向平滑。
初始化時，Prophet 會驗證以下參數：
func validate_inputs() # 驗證輸入
1. growth not in ('linear', 'logistic') # 驗證'growth'
2. (changepoint_range < 0) or (changepoint_range > 1) # 驗證changepoint_range區間
3. holidays lower_window/upper_window # 驗證holiday名稱和窗口值
4. seasonality_mode not in ['additive', 'multiplicative'] # 驗證'seasonality_mode'
以上，初始化的部分就結束了。
5.2 模型訓練
接下來便可以餵入數據訓練模型了：
m.fit(train_data)
其中，train_data是 pd.DataFrame 類型，必須包含兩列：
如果growth設為'logistic'，必須還有一列cap：
這裡的cap其實就是趨勢模型中的
餵入數據後，開始訓練Prophet模型，fit() 函數中的主要步驟如下：
初始化stan模型（prophet採用 pystan[5] 作為線性回歸模型）
簡單來說，這些步驟可以合併成最基本的機器學習範式：1-5步生成特徵矩陣（如下圖所示），6-7步訓練線性回歸模型（採用pystan）。
1. 擴充dataframe
首先通過 setup_dataframe() 函數為 train_data 增加三列，分別對時間和訓練數據進行歸一化。
df['floor']=0 # 最小值（默認）
df['t'] = (df['ds'] - self.start) / self.t_scale # 時間歸一化
df['y_scaled'] = (df['y'] - df['floor']) / self.y_scale # y 值歸一化
其中
t_scale = df['ds'].max() - df['ds'].min() # 最大時間和最小時間的差值
y_scale = (df['y'] - floor).abs() # y最大值的絕對值
train_data：
2. 設置季節性參數
若初始化模型時沒有設置季節性參數，Prophet 會通過 set_auto_seasonalities() 自動檢測訓練數據並設置季節性參數，主要有以下幾個參數：
yearly_disable：若訓練數據時間小於兩年，禁用年周期weekly_disable：若訓練數據時間小於兩周或min_dt>=一周，禁用周周期daily_disable：若訓練數據時間小於兩天或min_dt>=一天，禁用天周期parse_seasonality_args(): 獲取內置季節性的傅立葉級數（見下圖）
# prophet 內置的 [年/月/日] 周期季節性參數：
# period：周期
# fourier_order：傅立葉級數
# prior_scale：先驗值，用於調整擬合度
# mode：加法或乘法模型
self.seasonalities['yearly'] = {
        'period': 365.25,
        'fourier_order': fourier_order, # 默認為 10
        'prior_scale': self.seasonality_prior_scale, # 10
        'mode': self.seasonality_mode, # additive
        'condition_name': None
    }
self.seasonalities['weekly'] = {
        'period': 7,
        'fourier_order': fourier_order, # 默認為 3
        'prior_scale': self.seasonality_prior_scale, # 10
        'mode': self.seasonality_mode, # additive
        'condition_name': None
    }
self.seasonalities['daily'] = {
        'period': 1,
        'fourier_order': fourier_order, # 默認為 4
        'prior_scale': self.seasonality_prior_scale, # 10
        'mode': self.seasonality_mode, # additive
        'condition_name': None
    }
3. 生成季節性特徵
設定季節性參數後，通過 make_all_seasonality_features() 函數生成包含季節性特徵的dataframe，包括季節特徵、節日特徵和自定義特徵。make_all_seasonality_features的代碼結構如下所示：
make_seasonality_features(): 生成季節性特徵fourier_series()：生成傅立葉級數特徵make_holiday_features()：生成節日特徵make_holidays_df(): 根據國家和年份構建節假日dataframeconstruct_holiday_dataframe(): 構建節假日dataframemake_holiday_features()：構建節假日特徵additional regressors()：生成自定義特徵regressor_column_matrix()返回：seasonal_features, prior_scales, component_cols, modes
首先，通過 make_seasonality_features 構建季節性特徵，這裡主要採用傅立葉級數構建特徵（原理見第4節）。代碼如下：
# fourier_series: 生成傅立葉級數特徵
np.column_stack([
    fun((2.0 * (i + 1) * np.pi * t / period)) # t 為 1970-01-01 到現在的天數 array
    for i in range(series_order)
    for fun in (np.sin, np.cos)
])
下面兩張圖展示了年和周兩種周期的傅立葉級數特徵：
year_seasonality: 
week_seasonality: 
其次， make_holiday_features 函數根據節日配置參數構建節假日dataframe：holiday_df，如下表所示：
並根據節假日表生成對應的節假日特徵。
holiday_feature: 維的特徵，類似於one-hot表示
另外，Prophet還支持加入自定義特徵，源碼中 additional_regressors 函數綁定自定義特徵。
最後，合併seasonality、holiday、addtional特徵，生成總的季節性特徵seasonal_feature：
regressor_column_matrix 函數對 seasonal_feature 進行加工，生成每一維特徵對應的 one-hot 向量，方便之後的模型解釋和可視化：
modes：函數make_all_seasonality_features的返回結果，它是一個記錄加法變量和乘法變量的字典 (dict) 。
<class 'dict'>: 
{'additive': [
  'yearly',
  'weekly',
  'spring_festival',
  'double_11',
  'double_12',
  'yx_411',
  'yx_618',
  "New Year's Day",
  'Chinese New Year',
  'Tomb-Sweeping Day',
  'Labor Day',
  'Dragon Boat Festival',
  'Mid-Autumn Festival',
  'National Day',
  'additive_terms',
  'extra_regressors_additive',
  'holidays'],
'multiplicative': [
 'multiplicative_terms',
 'extra_regressors_multiplicative']
}
4. 設置變點
set_changepoint: 變點的設置有以下幾種情況。
如果自動設置變點，Prophet會默認在前80%的數據中選擇25個變點。
changepoints 示例如下：
5. 初始化變點參數
linear_growth_init：分段線性增長模型，初始化斜率
下圖為實際計算過程中
6. 初始化stan模型
Prophet預置了一份模型，初始化時會加載進來：
model = prophet_stan_model
以下是 stan 模型所用到的參數：
dat = {
    'T': history.shape[0],
    'K': seasonal_features.shape[1],
    'S': len(self.changepoints_t),
    'y': history['y_scaled'],
    't': history['t'],
    't_change': self.changepoints_t,
    'X': seasonal_features,
    'sigmas': prior_scales,
    'tau': self.changepoint_prior_scale,
    'trend_indicator': int(self.growth == 'logistic'),
    's_a': component_cols['additive_terms'],
    's_m': component_cols['multiplicative_terms'],
}
7. 訓練stan模型
這裡簡單介紹一下 Stan。Stan其實是一種用於統計推斷的概率程式語言。以線性回歸為例，下面是最簡單的線性回歸模型，具有單個預測變量，斜率和截距係數以及正態分布的噪聲。該模型可以編寫為：
上式等效於對於殘差的採樣：
並進一步簡化為：
Prophet的 stan 模型可直接寫為：
# stan分段線性模型
# y ~ 趨勢項 * （1 + 乘法項）+ 加法項（季節、節假日） + 殘差
y ~ normal(
  linear_trend(k, m, delta, t, A, t_change)
  .* (1 + X * (beta .* s_m))
  + X * (beta .* s_a),
  sigma_obs
)

vector linear_trend(
  real k, # 斜率
  real m, # 截距
  vector delta, # 增長率變化率向量
  vector t, # 時間
  matrix A, # 標誌位
  vector t_change 
) {
  return (k + A * delta) .* t + (m + A * (-t_change .* delta));
}

5.3 生成測試集
通過 make_future_dataframe 函數在訓練集基礎上往後延 periods 的時間作為測試集。
future = m.make_future_dataframe(periods=90)
5.4 預測
forecast = m.predict(future)
第一步，通過setup_dataframe()函數初始化測試集dataframe。
第二步，對測試集進行趨勢預測，這裡展示分段線法（piecewise_linear）的預測過程
再回顧一下分段線性模型的公式：
源碼中的趨勢計算如下：
# 梯度和截距累加
# Intercept changes
gammas = -changepoint_ts * deltas
# Get cumulative slope and intercept at each t
k_t = k * np.ones_like(t) # k 為初始斜率
m_t = m * np.ones_like(t) # m 為初始截距
for s, t_s in enumerate(changepoint_ts):
    indx = t >= t_s # 大於t_s的部分進行累加
    k_t[indx] += deltas[s] # 累加斜率
    m_t[indx] += gammas[s] # 累加截距
trend = k_t * t + m_t
trend = trend * self.y_scale + df['floor'] # scale還原
通過累加計算每一個變點之間的
第三步，周期性預測：通過 predict_seasonal_components 函數預測 seasonality, holidays 和 added regressors。
與訓練過程類似，通過 make_all_seasonality_features 函數構建預測日期的季節性特徵，生成類似於訓練集的特徵矩陣。
計算lower_p, upper_p，也就是預測結果的上下界。
for component in component_cols.columns:
    beta_c = self.params['beta'] * component_cols[component].values

    comp = np.matmul(X, beta_c.transpose())
    if component in self.component_modes['additive']:
        comp *= self.y_scale
    data[component] = np.nanmean(comp, axis=1)
    data[component + '_lower'] = np.nanpercentile(
        comp, lower_p, axis=1,
    )
    data[component + '_upper'] = np.nanpercentile(
        comp, upper_p, axis=1,
    )
計算每個周期性特徵分量 seasonal_components
第四步，計算趨勢和yhat的上下界（uncertainty）。
sample_posterior_predictive：
第五步，計算y_hat。
df2 = pd.concat((df[cols], intervals, seasonal_components), axis=1)
df2['yhat'] = (
        df2['trend'] * (1 + df2['multiplicative_terms'])
        + df2['additive_terms']
)
下圖為預測結果的可視化示例，分別展示了trend、holidays、weekly、yearly分量，提供了一定的可解釋性。
6. 總結
Prophet 針對一維時間序列進行建模，將時間序列分為趨勢、季節、節日三部分特徵，並輸出預測結果的上下界，同時保證了預測結果的可解釋性。
從原理上來說，Prophet是一個自動特徵提取的線性回歸模型，適用於少特徵、規律性高、單維度的時間序列場景。若業務場景有多維度特徵，多時間序列，複雜度高，則不適用Prophet。
7. 參考
[1]: Prophet官網：https://facebook.github.io/prophet/docs/quick_start.html#python-api
[2]: Prophet代碼地址：https://github.com/facebook/prophet
[3]: Prophet示例：https://github.com/facebook/prophet/tree/master/notebooks
[4]: Prophet預測模型論文：https://peerj.com/preprints/3190/
[5]: stan模型：https://pystan.readthedocs.io/en/latest/
[6]: Facebook 時間序列預測算法 Prophet 的研究. https://zhuanlan.zhihu.com/p/52330017
作者介紹
朱旭律，碩士，畢業於西安交通大學。2018年加入網易嚴選算法團隊負責銷量預測、廣告投放等項目的技術研發。