深入剖析時序Prophet模型：工作原理與源碼解析
隨著得物業務的快速發展，積累了大量的時序數據，這些數據對精細化運營，提升效率、降低成本有著重要作用。在得物的時序數據挖掘場景中，時序預測Prophet模型使用頻繁，本文對Prophet的原理和源碼進行深入分析，歡迎閱讀和交流。
一、引入
時間序列是指按照時間先后順序收集或觀測的一系列數據點，這類數據通常都具有一定時間相關性，基于這種順序性，我們可以對時間序列進行多種數據挖掘任務，包括分類、聚類、異常檢測和預測等。
時序任務在許多領域都有廣泛的應用，包括金融、商業、醫療、氣象、工業生產等。在得物的業務場景中，應用最為廣泛的是時序預測問題，本文介紹的內容主要和時序預測相關。
時序預測利用歷史時間序列數據構建數學統計、機器學習或深度學習等模型，來預測未來的觀測值，其目的是為了對時間序列的趨勢、周期性、季節性、特殊事件等規律進行捕捉并預測，從而指導業務人員做出商業、運營決策。以下是時間序列預測在得物各業務域的需求情況，總之哪里有指標，哪里就有預測。
解決時序預測問題的算法研究歷史非常悠久，是備受關注和持續發展的領域。根據算法原理和方法進行分類，時序預測模型可以分為以Holt-winters，ARIMA為代表的經典統計模型，用單一時序變量進行參數擬合；以線性回歸、樹回歸為代表的傳統機器學習算法，在有監督學習的框架下，構建特征來預測目標值；以及運用CNN、RNN、Transform等特征提取器，在Encoder-Decoder框架下進行預測的深度學習方法。
本文將介紹的Prophet模型和上述三種分類算法有所差異，但方法和原理上也結合了參數擬合、機器學習的思想，它將時序分解成趨勢變化、季節性、節假日、外部回歸因子之和，是結合時序分解的加法模型預測算法。
Prophet于2017年由Facebook’s Core Data Science team開源發布，盡管從時間上來看不是很新的模型，但是在得物實際的時序預測場景中取得了不俗的效果。目前網上的博客主要介紹了模型的基本原理、使用方式，在使用過程中筆者仍有一些疑問，例如：
- Prophet模型是如何進行訓練和預測？
- 模型如何進行概率預測，得到預測的上界和下界？
筆者帶著這些問題閱讀了Prophet的源碼，發現了一些論文中沒有提到的細節和流程，對以上問題有了答案，這里分享給大家。代碼參考Prophet python語言實現版本v1.1.5。
二、準備知識
如果讀者有機器學習基礎，熟悉最優化計算方法、貝葉斯估計、MCMC采樣可以跳過這一小節，不熟悉的讀者可以按照這一小節提到的概念搜索相關資料。
參數估計
通俗的說預測就是利用已知數據來推測產生該數據的模型和參數，然后用推測的模型和參數產生下一個結果。對于模型的參數估計方法，有頻率學派和貝葉斯學派之分。頻率學派認為模型參數是個固定的值。而貝葉斯學派是認為模型參數不是一個固定的值，而是源自某種潛在分布，希望從數據中推知該分布。
用大家最熟悉的線性回歸為例，假設噪聲服從正態分布：
在概率視角下，線性回歸可以表示為：
其中θ=(w0,w,σ^2)是模型所有的參數。
在頻率學派視角下θ是一個常量，可以用極大似然估計來估計參數值。其思想是對于N個觀測樣本來說使得其發生概率最大的參數就是最好的參數。整個觀測集發生的概率為：
由于連乘不容易參數求解，可以采用最大對數似然方法，線性回歸有解析解，我們可以直接求導計算得到最優的θ_hat：
在貝葉斯學派視角下θ是一個隨機變量而不是一個固定的值，由一個先驗分布進行約束，通過利用已知樣本進行統計推斷。貝葉斯學派認為一個好的θ不僅要考慮似然函數P(X|θ)還要考慮θ的先驗分布P(θ)，即最大化P(X|θ)P(θ)。由于X的先驗分布是固定的常數，最大化函數可以寫成P(X|θ)P(θ)/P(X)，根據貝葉斯公式可知P(X|θ)P(θ)/P(X)=P(θ|X)，即找到一個θ_hat能夠最大化后驗概率P(θ|X)。
求解參數θ可以選擇使用最大后驗估計（Maximum A Posterior, MAP）或者貝葉斯估計兩種方法，最大后驗估計直接求解讓后驗概率最大的θ。這樣問題就轉化還成一個約束問題，實際中可以使用梯度下降、牛頓法或者L-BFGS擬牛頓法等數值優化方法進行求解。
從最大似然估計和最大后驗估計來看求解的參數θ是一個確定值，但貝葉斯估計不是直接估計θ，而是估計θ的分布。在最大后驗估計中由于求θ極值過程中與P(X)無關，分母可以被忽略。但是在貝葉斯估計中是求整個后驗概率的分布，分母不能忽略。對于連續型隨機變量有：
則貝葉斯公式變為：
分母是積分形式一般沒有解析解，直接計算是非常困難的。于是引入馬爾可夫鏈蒙特卡洛算法（Markov Chain Monte Carlo, MCMC）進行采樣，從而得到各個參數的一系列采樣值。
Stan開源框架
https://mc-stan.org/ : Stan is a state-of-the-art platform for statistical modeling and high-performance statistical computation.
Stan是一個用于貝葉斯統計建模和推斷的開源框架，包含了定義概率模型的編程語言，以及高效執行貝葉斯推斷的算法庫，被廣泛用于統計建模和機器學習任務，它集成了上文談到的MAP、MCMC等求解算法，在Python中可以通過PyStan包來使用 Stan。
Prophet基本原理
Prophet基于時間序列分解和參數擬合來實現高效的訓練和預測。首先Prophet定義了一個關于時間步t的加法函數，函數值由趨勢項、季節項、節假日項、外部因子項以及觀測噪聲組成，核心公式如下：
其中g(t)代表趨勢項，s(t)代表季節項，h(t)代表節假日項（泛指外部回歸變量），ε_t代表誤差項。具體各項公式可以參考附錄中的文章以及官方論文。其中季節項、節假日項、外部因子項可以統一視為回歸因子，除了構造特征的方法不同以外，在模型訓練和預測階段都是一樣的處理方法。
Prophet其實是一種廣義線性模型，其數學形式是：
其中y是歷史label數據，x是Prophet中的加法回歸因子變量，α是趨勢項*（1+乘法回歸因子變量），β是回歸因子的權重，σ是噪聲服從高斯分布。模型訓練的就是公式中未知參數，Prophet的python代碼負責數據輸入、預處理、流程控制、可視化等部分功能，核心算法求解模塊調用Stan進行求解。Stan有自己定義的一套語言體系，定義了變量、輸入數據和模型，Prophet數學模型轉化成Stan的建模語言如下表示：
- 變量和數據
data {
int T; // Number of time periods
int<lower=1> K; // Number of regressors
vector[T] t; // Time
vector[T] cap; // Capacities for logistic trend
vector[T] y; // Time series
int S; // Number of changepoints
vector[S] t_change; // Times of trend changepoints
matrix[T,K] X; // Regressors
vector[K] sigmas; // Scale on seasonality prior
real<lower=0> tau; // Scale on changepoints prior
int trend_indicator; // 0 for linear, 1 for logistic, 2 for flat
vector[K] s_a; // Indicator of additive features
vector[K] s_m; // Indicator of multiplicative features
}- 模型
model {
//priors
k ~ normal(0, 5);
m ~ normal(0, 5);
delta ~ double_exponential(0, tau);
sigma_obs ~ normal(0, 0.5);
beta ~ normal(0, sigmas);
// Likelihood
y ~ normal_id_glm(
X_sa,
trend .* (1 + X_sm * beta),
beta,
sigma_obs
);
}論文中各項成分的參數預先設定了先驗分布，根據：
能夠得到后驗概率，根據第二部分提到的參數估計方法，運用最大后驗估計或MCMC采樣即可得到參數估計值。
三、Prophet源碼剖析
代碼結構
Prophet代碼目錄如下所示，核心邏輯主要是Prophet.stan、models.py、forecaster.py三個腳本，本文的分享主要圍繞這三個腳本。
|-- stan
|-- Prophet.stan # stan語言實現的Prophet模型
|-- setup.py
|-- Prophet
|-- serialize.py # Prophet訓練完成的模型的序列化、模型的存儲和加載，以便后續重用。
|-- plot.py # 可視化模塊
|-- models.py # stan模型定義
|-- diagnostics.py # 用于診斷 Prophet 模型性能的函數，比如計算誤差指標、繪制診斷圖等。
|-- tests # 封裝的測試腳本
|-- __init__.py
|-- __version__.py
|-- forecaster.py # Prophet 類的定義，模型的核心邏輯，數據的處理、模型的擬合和預測等功能
|-- utilities.py # 包含了一些輔助函數，用于處理時間序列數據、特征工程等。
|-- make_holidays.py # 節假日具體日期的構造- Prophet.stan：用Stan編程語言實現的Prophet模型腳本。
- models.py：python和Stan語言交互的模塊，定義使用Stan腳本語言時的輸入數據，擬合參數、輸出數據，控制Stan進行參數擬合、采樣。
- forecaster.py：Prophet模型核心代碼，定義了Prophet類，涉及整個模型的數據處理、模型訓練和模型推理等功能。
模型框架
Prophet模型訓練-預測流程如下所示：
初始參數
- Prophet類初始化
初始化模塊主要是讀取預先設定和Prophet有關的模型參數。分別有以下幾類：
趨勢項相關
growth: String 'linear', 'logistic' or 'flat' to specify a linear, logistic or
flat trend.
changepoints: List of dates at which to include potential changepoints. If
not specified, potential changepoints are selected automatically.
n_changepoints: Number of potential changepoints to include. Not used
if input `changepoints` is supplied. If `changepoints` is not supplied,
then n_changepoints potential changepoints are selected uniformly from
the first `changepoint_range` proportion of the history.
changepoint_range: Proportion of history in which trend changepoints will
be estimated. Defaults to 0.8 for the first 80%. Not used if
`changepoints` is specified.
changepoint_prior_scale: Parameter modulating the flexibility of the
automatic changepoint selection. Large values will allow many
changepoints, small values will allow few changepoints.- 季節性相關
yearly_seasonality: Fit yearly seasonality.
Can be 'auto', True, False, or a number of Fourier terms to generate.
weekly_seasonality: Fit weekly seasonality.
Can be 'auto', True, False, or a number of Fourier terms to generate.
daily_seasonality: Fit daily seasonality.
Can be 'auto', True, False, or a number of Fourier terms to generate.
holidays: pd.DataFrame with columns holiday (string) and ds (date type)
and optionally columns lower_window and upper_window which specify a
range of days around the date to be included as holidays.
lower_window=-2 will include 2 days prior to the date as holidays. Also
optionally can have a column prior_scale specifying the prior scale for
that holiday.
seasonality_mode: 'additive' (default) or 'multiplicative'.
seasonality_prior_scale: Parameter modulating the strength of the
seasonality model. Larger values allow the model to fit larger seasonal
fluctuations, smaller values dampen the seasonality. Can be specified
for individual seasonalities using add_seasonality.
holidays_prior_scale: Parameter modulating the strength of the holiday
components model, unless overridden in the holidays input.
holidays_mode: 'additive' or 'multiplicative'. Defaults to seasonality_mode.- 參數擬合相關
growth: String 'linear', 'logistic' or 'flat' to specify a linear, logistic or
flat trend.
changepoints: List of dates at which to include potential changepoints. If
not specified, potential changepoints are selected automatically.
n_changepoints: Number of potential changepoints to include. Not used
if input `changepoints` is supplied. If `changepoints` is not supplied,
then n_changepoints potential changepoints are selected uniformly from
the first `changepoint_range` proportion of the history.
changepoint_range: Proportion of history in which trend changepoints will
be estimated. Defaults to 0.8 for the first 80%. Not used if
`changepoints` is specified.- 不確定性相關
mcmc_samples: Integer, if greater than 0, will do full Bayesian inference
with the specified number of MCMC samples. If 0, will do MAP
estimation.
interval_width: Float, width of the uncertainty intervals provided
for the forecast. If mcmc_samples=0, this will be only the uncertainty
in the trend using the MAP estimate of the extrapolated generative
model. If mcmc.samples>0, this will be integrated over all model
parameters, which will include uncertainty in seasonality.
uncertainty_samples: Number of simulated draws used to estimate
uncertainty intervals. Settings this value to 0 or False will disable
uncertainty estimation and speed up the calculation.
stan_backend: str as defined in StanBackendEnum default: None - will try to訓練過程
數據預處理
- 數據檢查
python.Prophet.forecaster.Prophet.setup_dataframe
Prophet會檢查模型輸入參數、數據的合法性，幫助使用者快速定位問題。比如：
檢查輸入數據中是否有y列，ds列是否符合時間輸入規范，是否有缺失值。
檢查添加的額外回歸項，是否有缺失值，是否輸入數據中有添加的回歸項數據。
季節性項的condition列值是否是bool值等，用于控制模型選擇性的學習condition為True的樣本。
- 數據歸一化
python.Prophet.forecaster.Prophet.initialize_scales
對于y，有Absmax和Minmax兩種方式：現實場景中特征值或y，往往有很大的量綱差異，會降低模型的性能、預測的準確性。Prophet內置了對y和對外部回歸因子add regressors的歸一化。
AbsMax歸一化：
含義：AbsMax歸一化是將原始數據縮放到[-1, 1]的范圍內，使數據的絕對值最大值為1。
適用場景：AbsMax歸一化適用于數據中存在明顯的異常值或極端值的情況，可以保留數據的分布形狀并減少異常值對模型的影響。
MinMax歸一化：
含義：MinMax歸一化是將原始數據縮放到[0, 1]的范圍內，使數據的最小值對應0，最大值對應1。
適用場景：MinMax歸一化適用于需要將數據映射到一個固定范圍內的情況，保留了原始數據的相對關系，在大部分機器學習算法中經常用到。
if self.scaling == "absmax":
self.y_min = 0.
self.y_scale = float((df['y']).abs().max())
elif self.scaling == "minmax":
self.y_min = df['y'].min()
self.y_scale = float(df['y'].max() - self.y_min)
...
...對于add regressor特征則進行標準化，將特征數據縮放到均值為0，標準差為1的標準正態分布。
python/Prophet/forecaster.py:389
計算得到的y_scale和回歸項的均值和標準差，會更新到Prophet類中，供預測階段使用。
自動設置周期性
python.Prophet.forecaster.Prophet.set_auto_seasonalities
如果在初始化Prophet類時，沒有指定季節性相關的參數，則會根據數據長度和間隔自動增加季節性項，三種不同的周期性默認的傅立葉階數為10。
當歷史數據大于等于兩年，則增加yearly seasonality，默認周期為365.25；
當歷史數據大于等于兩周，且相鄰兩數據之間的時間之差小于1周，則打開weekly seasonality，默認周期為7；
當歷史數據大于等于2天，且相鄰兩數據之間的時間小于1天，則打開daily seasonality，默認周期為1；生成特征寬表
Prophet模型的成分中劃分了seasonality features、holiday features和add regressors，但其實是一樣的處理方法。根據Prophet類初始化時的參數，生成特征寬表，行為依次遞增的時間序列，列為每一個feature對應的值。
對于某一個seasonality特征，根據傳入的周期性和傅立葉階數，生成不同列，列數等于傅立葉階數，列值等于某一階的周期性函數值。然后結合該特征對應的condition_name的值，得到最終的特征值。condition_name是True、False的np.array類型，用來指導模型學習condition_name列中為True的時間段數據，忽略為False的時間段的數據。
對于holiday features和add regressors，則對于每一個特征，生成一列，列為每一個feature對應的值。特征寬表的特征總數為seasonality特征的傅立葉階數+holiday features數+add regressors數。
在生成特征寬表的同時，Prophet會定義component_cols變量，來維護哪些列同屬于同一個成分。比如把同一個季節性項的多個周期函數合并成一個成分，把同一個節假日不同時間的因子合并成一個成分，把加性或者乘性因子項合并成一個成分，這樣可以方便進行模型的成分分析和結果可視化。比如給定一個Prophet模型：
m = Prophet(
holidays=pd.DataFrame({
'holiday': '618',
'ds': pd.to_datetime(['2021-06-18', '2022-06-18', '2023-06-18', '2024-06-18']),
'lower_window': -1,
'upper_window': 0,
'mode': 'additive'
}))
model.add_seasonality({'name':'seasonal_1','period':7, 'fourier_order':1,'mode': 'additive'})
model.add_seasonality({'name':'seasonal_2','period':365, 'fourier_order':2,'mode': 'multiplicative'})假設yearly傅立葉階數為1的additive，weekly階數為2的additive函數，那么component_cols為：
設置間斷點
python.Prophet.forecaster.Prophet.set_changepoints
首先當有人工設置間斷點時，以設置的間斷點為主，間斷點必須在訓練數據之內，否則會報錯。
當沒有設置間斷點時，Prophet會根據初始化的參數n_changepoints間斷點數和changepoint_range間斷點篩選范圍，進行自動采樣。例如時間序列有100個點，n_changepoints=10，changepoint_range=0.8，那么Prophet會在前80個點中，等間隔的采樣10個點，作為候選間斷點。設置的間斷點數量不能超過changepoint_range范圍內的訓練集數量。
初始化stan參數
python.Prophet.forecaster.Prophet.calculate_initial_params
使用貝葉斯估計參數，需要給定參數的初始值，然后迭代至算法收斂，Prophet用全部訓練數據計算初始參數。每一項成分都需要進行參數初始化。
以線性趨勢為例，用標準化的y計算線性函數的斜率和偏置。其他回歸項因子β，突變點增長系數δ都設置為0。
參數擬合
在初始化階段得到的參數和數據，轉化成字典形式傳入Stan進行參數擬合，字典格式如下所示：
字典K-V值和Prophet/stan目錄下的Prophet.stan腳本相對應，然后通過Pystan調用Stan引擎求解參數，擬合后的參數存儲在Propeht.params中，求解參數有兩種方式：
- 牛頓法/LBFGS法進行最大后驗估計
python.Prophet.models.CmdStanPyBackend.fit
- 馬爾可夫鏈蒙特卡羅法算法進行參數估計
python.Prophet.models.CmdStanPyBackend.sampling
預測過程
構建預測dataframe
python.Prophet.forecaster.Prophet.make_future_dataframe
我們可以通過Prophet類中的make_future_dataframe函數構建未來預測的DataFrame，預測長度、預測頻率由初始化參數periods、freq設置。
構建預測的DataFrame之后，和在訓練過程一致，需經過數據檢查、歸一化，歸一化使用訓練過程中計算得到的y_scale和回歸項的均值和標準差。
點預測
點預測也就是在時間點上輸出單個預測值，預測值由趨勢預測和回歸項預測組成。
- 趨勢預測
python.Prophet.forecaster.Prophet.predict_trend
在貝葉斯回歸中，未知參數服從一個指定的先驗分布，Prophet使用Stan引擎計算得到的返回參數的期望作為趨勢項公式的帶入值。
def predict_trend(self, df):
k = np.nanmean(self.params['k'])
m = np.nanmean(self.params['m'])
deltas = np.nanmean(self.params['delta'], axis=0)
t = np.array(df['t'])
if self.growth == 'linear':
trend = self.piecewise_linear(t, deltas, k, m, self.changepoints_t)
elif self.growth == 'logistic':
cap = df['cap_scaled']
trend = self.piecewise_logistic(
t, cap, deltas, k, m, self.changepoints_t)
elif self.growth == 'flat':
# constant trend
trend = self.flat_trend(t, m)
return trend * self.y_scale + df['floor']對于flat、linear、logistic不同的趨勢分量，結合changepoints_t進行趨勢的預測，以線性趨勢為例。
python.Prophet.forecaster.Prophet.piecewise_linear
即論文中描述的公式：
最后根據訓練過程計算的y_scale縮放變量還原量級，作為最終的趨勢預測。
- 回歸項預測
python.Prophet.forecaster.Prophet.predict_seasonal_components
上文提到的Prophet的seasonality features、holiday features和add regressors項，在參數求解過程中是一樣處理方法。
β即為擬合出的權重系數，β乘以特征寬表對應的值，得到每個特征的效應值。根據component_cols維護的不同成分的特征歸屬。就可以得到seasonality、holiday features、add regressors、multiplicative_terms、additive_terms各個部分的效應值。這些在模型解釋性、結果可視化、偏差分析中都可以應用。
- 預測匯總
python/Prophet/forecaster.py:1287
預測成分匯總步驟把趨勢、季節性、節假日等預測效應值進行拼接，最后的預測結果公式如下：
預測結果值=趨勢*(1+乘法因素)+加法因素
預測的不確定性
Prophet通過觀測噪聲、參數不確定性以及未來趨勢不確定性三個部分輸出預測值的不確定性，最終通過yhat_lower、yhat_upper體現。
- 觀測噪聲不確定性
python/Prophet/forecaster.py:1556
一般時間序列分解模型，最后的一個分解項就是噪聲，Prophet也是這樣的，假設了輸入數據點含有服從正態分布的觀測噪聲，Stan模型腳本中的sigma_obs就是擬合了這部分噪聲。觀測噪聲不確定性為：
sigma = self.params['sigma_obs'][iteration]
noise_terms = np.random.normal(0, sigma, trends.shape) * self.y_scale變量sigma是sigma_obs的估計值。在每個時間步驟上，從正態分布中采樣n_samples次，得到噪聲不確定項目noise_terms。
- 趨勢不確定性
python.Prophet.forecaster.Prophet._make_trend_shift_matrix
對于趨勢的不確定性，官方文檔中解釋：預測中最大的不確定性來源在于未來趨勢變化的可能性。因此我們盡最大可能采取了一個合理的假設，假設未來的趨勢變化會與歷史上觀察到的相似。特別是，我們假設未來的趨勢變化的平均頻率和幅度將與歷史上觀察到的相同。我們將這些趨勢變化進行投影，并通過計算它們的分布來獲得不確定性區間。
具體而言，趨勢不確性由突變點出現的位置和突變的比例確定。首先計算歷史上突變點出現的間隔的均值*歷史上突變點的個數，得到每個時間點上產生突變點概率likelihood。然后計算歷史突變點的變化率的絕對值的均值mean_delta。最后進行采樣來模擬不確定性。均勻分布采樣輸出矩陣（n_samples, future_length），對每一個元素判斷是否小于似然值likelihood，來判斷是否發生了突變，如果小于則發生了突變，然后輸出均值為mean_delta的拉普拉斯分布，得到突變點變化率的采樣值。
- 參數的不確定性
python.Prophet.forecaster.Prophet.predict_uncertainty
當初始化參數MCMC=0，即調用Stan用MAP算法，此時算法得出的是收斂時的最終估計值，不確定性由觀測不確定性和趨勢不確定性組成，參數沒有不確定性的估計。
MCMC>0時候，除了觀測噪聲、趨勢不確定性外，由于Stan采用MCMC采樣方法，得到的返回值是待估計參數值的若干個采樣點，根據每個參數的采樣點得到不同的預測值。
- 不確定性結果輸出
進一步，由于每一個時間步有以上不確定性成分的n_samples個采樣數據。根據預設的置信度寬度，得到上、下分位數值：
lower_p = 100 * (1.0 - self.interval_width) / 2
upper_p = 100 * (1.0 + self.interval_width) / 2對于每個時間步驟，求lower_p分位數、upper_p分位數，即得到yhat_lower、yhat_upper。
四、總結
Prophet框架在有明顯規律的單變量時序預測場景中有著非常不錯的表現。在得物的多個業務場景中已經得到了驗證。閱讀源碼能夠幫助我們更好的了解模型細節，在模型優化時幫助我們有的放矢，在模型選型時幫助我們掌握模型的適應性和局限性。
參考資料
??https://peerj.com/preprints/3190/??
??https://github.com/facebook/prophet??
??https://zhuanlan.zhihu.com/p/37543542??
本文轉載自 ??得物技術??，作者： 探流
