隨著得物業(yè)務(wù)的快速發(fā)展，積累了大量的時序數(shù)據(jù)，這些數(shù)據(jù)對精細化運營，提升效率、降低成本有著重要作用。在得物的時序數(shù)據(jù)挖掘場景中，時序預(yù)測Prophet模型使用頻繁，本文對Prophet的原理和源碼進行深入分析，歡迎閱讀和交流。

一、引入

時間序列是指按照時間先后順序收集或觀測的一系列數(shù)據(jù)點，這類數(shù)據(jù)通常都具有一定時間相關(guān)性，基于這種順序性，我們可以對時間序列進行多種數(shù)據(jù)挖掘任務(wù)，包括分類、聚類、異常檢測和預(yù)測等。

時序任務(wù)在許多領(lǐng)域都有廣泛的應(yīng)用，包括金融、商業(yè)、醫(yī)療、氣象、工業(yè)生產(chǎn)等。在得物的業(yè)務(wù)場景中，應(yīng)用最為廣泛的是時序預(yù)測問題，本文介紹的內(nèi)容主要和時序預(yù)測相關(guān)。

時序預(yù)測利用歷史時間序列數(shù)據(jù)構(gòu)建數(shù)學(xué)統(tǒng)計、機器學(xué)習(xí)或深度學(xué)習(xí)等模型，來預(yù)測未來的觀測值，其目的是為了對時間序列的趨勢、周期性、季節(jié)性、特殊事件等規(guī)律進行捕捉并預(yù)測，從而指導(dǎo)業(yè)務(wù)人員做出商業(yè)、運營決策。以下是時間序列預(yù)測在得物各業(yè)務(wù)域的需求情況，總之哪里有指標(biāo)，哪里就有預(yù)測。

解決時序預(yù)測問題的算法研究歷史非常悠久，是備受關(guān)注和持續(xù)發(fā)展的領(lǐng)域。根據(jù)算法原理和方法進行分類，時序預(yù)測模型可以分為以Holt-winters，ARIMA為代表的經(jīng)典統(tǒng)計模型，用單一時序變量進行參數(shù)擬合；以線性回歸、樹回歸為代表的傳統(tǒng)機器學(xué)習(xí)算法，在有監(jiān)督學(xué)習(xí)的框架下，構(gòu)建特征來預(yù)測目標(biāo)值；以及運用CNN、RNN、Transform等特征提取器，在Encoder-Decoder框架下進行預(yù)測的深度學(xué)習(xí)方法。

本文將介紹的Prophet模型和上述三種分類算法有所差異，但方法和原理上也結(jié)合了參數(shù)擬合、機器學(xué)習(xí)的思想，它將時序分解成趨勢變化、季節(jié)性、節(jié)假日、外部回歸因子之和，是結(jié)合時序分解的加法模型預(yù)測算法。

Prophet于2017年由Facebook’s Core Data Science team開源發(fā)布，盡管從時間上來看不是很新的模型，但是在得物實際的時序預(yù)測場景中取得了不俗的效果。目前網(wǎng)上的博客主要介紹了模型的基本原理、使用方式，在使用過程中筆者仍有一些疑問，例如：

• Prophet模型是如何進行訓(xùn)練和預(yù)測？
• 模型如何進行概率預(yù)測，得到預(yù)測的上界和下界？

筆者帶著這些問題閱讀了Prophet的源碼，發(fā)現(xiàn)了一些論文中沒有提到的細節(jié)和流程，對以上問題有了答案，這里分享給大家。代碼參考Prophet python語言實現(xiàn)版本v1.1.5。

二、準(zhǔn)備知識

如果讀者有機器學(xué)習(xí)基礎(chǔ)，熟悉最優(yōu)化計算方法、貝葉斯估計、MCMC采樣可以跳過這一小節(jié)，不熟悉的讀者可以按照這一小節(jié)提到的概念搜索相關(guān)資料。

參數(shù)估計

通俗的說預(yù)測就是利用已知數(shù)據(jù)來推測產(chǎn)生該數(shù)據(jù)的模型和參數(shù)，然后用推測的模型和參數(shù)產(chǎn)生下一個結(jié)果。對于模型的參數(shù)估計方法，有頻率學(xué)派和貝葉斯學(xué)派之分。頻率學(xué)派認(rèn)為模型參數(shù)是個固定的值。而貝葉斯學(xué)派是認(rèn)為模型參數(shù)不是一個固定的值，而是源自某種潛在分布，希望從數(shù)據(jù)中推知該分布。

用大家最熟悉的線性回歸為例，假設(shè)噪聲服從正態(tài)分布：

在概率視角下，線性回歸可以表示為：

其中θ=(w0,w,σ^2)是模型所有的參數(shù)。

在頻率學(xué)派視角下θ是一個常量，可以用極大似然估計來估計參數(shù)值。其思想是對于N個觀測樣本來說使得其發(fā)生概率最大的參數(shù)就是最好的參數(shù)。整個觀測集發(fā)生的概率為：

由于連乘不容易參數(shù)求解，可以采用最大對數(shù)似然方法，線性回歸有解析解，我們可以直接求導(dǎo)計算得到最優(yōu)的θ_hat：

在貝葉斯學(xué)派視角下θ是一個隨機變量而不是一個固定的值，由一個先驗分布進行約束，通過利用已知樣本進行統(tǒng)計推斷。貝葉斯學(xué)派認(rèn)為一個好的θ不僅要考慮似然函數(shù)P(X|θ)還要考慮θ的先驗分布P(θ)，即最大化P(X|θ)P(θ)。由于X的先驗分布是固定的常數(shù)，最大化函數(shù)可以寫成P(X|θ)P(θ)/P(X)，根據(jù)貝葉斯公式可知P(X|θ)P(θ)/P(X)=P(θ|X)，即找到一個θ_hat能夠最大化后驗概率P(θ|X)。

求解參數(shù)θ可以選擇使用最大后驗估計（Maximum A Posterior, MAP）或者貝葉斯估計兩種方法，最大后驗估計直接求解讓后驗概率最大的θ。這樣問題就轉(zhuǎn)化還成一個約束問題，實際中可以使用梯度下降、牛頓法或者L-BFGS擬牛頓法等數(shù)值優(yōu)化方法進行求解。

從最大似然估計和最大后驗估計來看求解的參數(shù)θ是一個確定值，但貝葉斯估計不是直接估計θ，而是估計θ的分布。在最大后驗估計中由于求θ極值過程中與P(X)無關(guān)，分母可以被忽略。但是在貝葉斯估計中是求整個后驗概率的分布，分母不能忽略。對于連續(xù)型隨機變量有：

則貝葉斯公式變?yōu)椋?/p>

分母是積分形式一般沒有解析解，直接計算是非常困難的。于是引入馬爾可夫鏈蒙特卡洛算法（Markov Chain Monte Carlo, MCMC）進行采樣，從而得到各個參數(shù)的一系列采樣值。

Stan開源框架

https://mc-stan.org/ : Stan is a state-of-the-art platform for statistical modeling and high-performance statistical computation.

Stan是一個用于貝葉斯統(tǒng)計建模和推斷的開源框架，包含了定義概率模型的編程語言，以及高效執(zhí)行貝葉斯推斷的算法庫，被廣泛用于統(tǒng)計建模和機器學(xué)習(xí)任務(wù)，它集成了上文談到的MAP、MCMC等求解算法，在Python中可以通過PyStan包來使用 Stan。

Prophet基本原理

Prophet基于時間序列分解和參數(shù)擬合來實現(xiàn)高效的訓(xùn)練和預(yù)測。首先Prophet定義了一個關(guān)于時間步t的加法函數(shù)，函數(shù)值由趨勢項、季節(jié)項、節(jié)假日項、外部因子項以及觀測噪聲組成，核心公式如下：

其中g(shù)(t)代表趨勢項，s(t)代表季節(jié)項，h(t)代表節(jié)假日項（泛指外部回歸變量），ε_t代表誤差項。具體各項公式可以參考附錄中的文章以及官方論文。其中季節(jié)項、節(jié)假日項、外部因子項可以統(tǒng)一視為回歸因子，除了構(gòu)造特征的方法不同以外，在模型訓(xùn)練和預(yù)測階段都是一樣的處理方法。

Prophet其實是一種廣義線性模型，其數(shù)學(xué)形式是：

其中y是歷史label數(shù)據(jù)，x是Prophet中的加法回歸因子變量，α是趨勢項*（1+乘法回歸因子變量），β是回歸因子的權(quán)重，σ是噪聲服從高斯分布。模型訓(xùn)練的就是公式中未知參數(shù)，Prophet的python代碼負(fù)責(zé)數(shù)據(jù)輸入、預(yù)處理、流程控制、可視化等部分功能，核心算法求解模塊調(diào)用Stan進行求解。Stan有自己定義的一套語言體系，定義了變量、輸入數(shù)據(jù)和模型，Prophet數(shù)學(xué)模型轉(zhuǎn)化成Stan的建模語言如下表示：

• 變量和數(shù)據(jù)

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
  );
}

論文中各項成分的參數(shù)預(yù)先設(shè)定了先驗分布，根據(jù)：

能夠得到后驗概率，根據(jù)第二部分提到的參數(shù)估計方法，運用最大后驗估計或MCMC采樣即可得到參數(shù)估計值。

三、Prophet源碼剖析

代碼結(jié)構(gòu)

Prophet代碼目錄如下所示，核心邏輯主要是Prophet.stan、models.py、forecaster.py三個腳本，本文的分享主要圍繞這三個腳本。

|-- stan
    |-- Prophet.stan # stan語言實現(xiàn)的Prophet模型
|-- setup.py
|-- Prophet
    |-- serialize.py # Prophet訓(xùn)練完成的模型的序列化、模型的存儲和加載，以便后續(xù)重用。
    |-- plot.py # 可視化模塊
    |-- models.py # stan模型定義
    |-- diagnostics.py # 用于診斷 Prophet 模型性能的函數(shù)，比如計算誤差指標(biāo)、繪制診斷圖等。
    |-- tests # 封裝的測試腳本
    |-- __init__.py
    |-- __version__.py
    |-- forecaster.py #  Prophet 類的定義，模型的核心邏輯，數(shù)據(jù)的處理、模型的擬合和預(yù)測等功能
    |-- utilities.py # 包含了一些輔助函數(shù)，用于處理時間序列數(shù)據(jù)、特征工程等。
    |-- make_holidays.py # 節(jié)假日具體日期的構(gòu)造

• Prophet.stan：用Stan編程語言實現(xiàn)的Prophet模型腳本。
• models.py：python和Stan語言交互的模塊，定義使用Stan腳本語言時的輸入數(shù)據(jù)，擬合參數(shù)、輸出數(shù)據(jù)，控制Stan進行參數(shù)擬合、采樣。
• forecaster.py：Prophet模型核心代碼，定義了Prophet類，涉及整個模型的數(shù)據(jù)處理、模型訓(xùn)練和模型推理等功能。

模型框架

Prophet模型訓(xùn)練-預(yù)測流程如下所示：

初始參數(shù)

• Prophet類初始化
  初始化模塊主要是讀取預(yù)先設(shè)定和Prophet有關(guān)的模型參數(shù)。分別有以下幾類：

趨勢項相關(guān)

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.

• 季節(jié)性相關(guān)

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.

• 參數(shù)擬合相關(guān)

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.

• 不確定性相關(guān)

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

訓(xùn)練過程

數(shù)據(jù)預(yù)處理

• 數(shù)據(jù)檢查

python.Prophet.forecaster.Prophet.setup_dataframe

Prophet會檢查模型輸入?yún)?shù)、數(shù)據(jù)的合法性，幫助使用者快速定位問題。比如：

檢查輸入數(shù)據(jù)中是否有y列，ds列是否符合時間輸入規(guī)范，是否有缺失值。

檢查添加的額外回歸項，是否有缺失值，是否輸入數(shù)據(jù)中有添加的回歸項數(shù)據(jù)。

季節(jié)性項的condition列值是否是bool值等，用于控制模型選擇性的學(xué)習(xí)condition為True的樣本。

• 數(shù)據(jù)歸一化

python.Prophet.forecaster.Prophet.initialize_scales

對于y，有Absmax和Minmax兩種方式：現(xiàn)實場景中特征值或y，往往有很大的量綱差異，會降低模型的性能、預(yù)測的準(zhǔn)確性。Prophet內(nèi)置了對y和對外部回歸因子add regressors的歸一化。

AbsMax歸一化：
含義：AbsMax歸一化是將原始數(shù)據(jù)縮放到[-1, 1]的范圍內(nèi)，使數(shù)據(jù)的絕對值最大值為1。
適用場景：AbsMax歸一化適用于數(shù)據(jù)中存在明顯的異常值或極端值的情況，可以保留數(shù)據(jù)的分布形狀并減少異常值對模型的影響。

MinMax歸一化：
含義：MinMax歸一化是將原始數(shù)據(jù)縮放到[0, 1]的范圍內(nèi)，使數(shù)據(jù)的最小值對應(yīng)0，最大值對應(yīng)1。
適用場景：MinMax歸一化適用于需要將數(shù)據(jù)映射到一個固定范圍內(nèi)的情況，保留了原始數(shù)據(jù)的相對關(guān)系，在大部分機器學(xué)習(xí)算法中經(jīng)常用到。

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特征則進行標(biāo)準(zhǔn)化，將特征數(shù)據(jù)縮放到均值為0，標(biāo)準(zhǔn)差為1的標(biāo)準(zhǔn)正態(tài)分布。

python/Prophet/forecaster.py:389

計算得到的y_scale和回歸項的均值和標(biāo)準(zhǔn)差，會更新到Prophet類中，供預(yù)測階段使用。

自動設(shè)置周期性

python.Prophet.forecaster.Prophet.set_auto_seasonalities

如果在初始化Prophet類時，沒有指定季節(jié)性相關(guān)的參數(shù)，則會根據(jù)數(shù)據(jù)長度和間隔自動增加季節(jié)性項，三種不同的周期性默認(rèn)的傅立葉階數(shù)為10。

當(dāng)歷史數(shù)據(jù)大于等于兩年，則增加yearly seasonality，默認(rèn)周期為365.25；
當(dāng)歷史數(shù)據(jù)大于等于兩周，且相鄰兩數(shù)據(jù)之間的時間之差小于1周，則打開weekly seasonality，默認(rèn)周期為7；
當(dāng)歷史數(shù)據(jù)大于等于2天，且相鄰兩數(shù)據(jù)之間的時間小于1天，則打開daily seasonality，默認(rèn)周期為1；

生成特征寬表

Prophet模型的成分中劃分了seasonality features、holiday features和add regressors，但其實是一樣的處理方法。根據(jù)Prophet類初始化時的參數(shù)，生成特征寬表，行為依次遞增的時間序列，列為每一個feature對應(yīng)的值。

對于某一個seasonality特征，根據(jù)傳入的周期性和傅立葉階數(shù)，生成不同列，列數(shù)等于傅立葉階數(shù)，列值等于某一階的周期性函數(shù)值。然后結(jié)合該特征對應(yīng)的condition_name的值，得到最終的特征值。condition_name是True、False的np.array類型，用來指導(dǎo)模型學(xué)習(xí)condition_name列中為True的時間段數(shù)據(jù)，忽略為False的時間段的數(shù)據(jù)。

對于holiday features和add regressors，則對于每一個特征，生成一列，列為每一個feature對應(yīng)的值。特征寬表的特征總數(shù)為seasonality特征的傅立葉階數(shù)+holiday features數(shù)+add regressors數(shù)。

在生成特征寬表的同時，Prophet會定義component_cols變量，來維護哪些列同屬于同一個成分。比如把同一個季節(jié)性項的多個周期函數(shù)合并成一個成分，把同一個節(jié)假日不同時間的因子合并成一個成分，把加性或者乘性因子項合并成一個成分，這樣可以方便進行模型的成分分析和結(jié)果可視化。比如給定一個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'})

假設(shè)yearly傅立葉階數(shù)為1的additive，weekly階數(shù)為2的additive函數(shù)，那么component_cols為：

設(shè)置間斷點

python.Prophet.forecaster.Prophet.set_changepoints

首先當(dāng)有人工設(shè)置間斷點時，以設(shè)置的間斷點為主，間斷點必須在訓(xùn)練數(shù)據(jù)之內(nèi)，否則會報錯。

當(dāng)沒有設(shè)置間斷點時，Prophet會根據(jù)初始化的參數(shù)n_changepoints間斷點數(shù)和changepoint_range間斷點篩選范圍，進行自動采樣。例如時間序列有100個點，n_changepoints=10，changepoint_range=0.8，那么Prophet會在前80個點中，等間隔的采樣10個點，作為候選間斷點。設(shè)置的間斷點數(shù)量不能超過changepoint_range范圍內(nèi)的訓(xùn)練集數(shù)量。

初始化stan參數(shù)

python.Prophet.forecaster.Prophet.calculate_initial_params

使用貝葉斯估計參數(shù)，需要給定參數(shù)的初始值，然后迭代至算法收斂，Prophet用全部訓(xùn)練數(shù)據(jù)計算初始參數(shù)。每一項成分都需要進行參數(shù)初始化。

以線性趨勢為例，用標(biāo)準(zhǔn)化的y計算線性函數(shù)的斜率和偏置。其他回歸項因子β，突變點增長系數(shù)δ都設(shè)置為0。

參數(shù)擬合

在初始化階段得到的參數(shù)和數(shù)據(jù)，轉(zhuǎn)化成字典形式傳入Stan進行參數(shù)擬合，字典格式如下所示：

字典K-V值和Prophet/stan目錄下的Prophet.stan腳本相對應(yīng)，然后通過Pystan調(diào)用Stan引擎求解參數(shù)，擬合后的參數(shù)存儲在Propeht.params中，求解參數(shù)有兩種方式：

• 牛頓法/LBFGS法進行最大后驗估計

python.Prophet.models.CmdStanPyBackend.fit

• 馬爾可夫鏈蒙特卡羅法算法進行參數(shù)估計

python.Prophet.models.CmdStanPyBackend.sampling

預(yù)測過程

構(gòu)建預(yù)測dataframe

python.Prophet.forecaster.Prophet.make_future_dataframe

我們可以通過Prophet類中的make_future_dataframe函數(shù)構(gòu)建未來預(yù)測的DataFrame，預(yù)測長度、預(yù)測頻率由初始化參數(shù)periods、freq設(shè)置。

構(gòu)建預(yù)測的DataFrame之后，和在訓(xùn)練過程一致，需經(jīng)過數(shù)據(jù)檢查、歸一化，歸一化使用訓(xùn)練過程中計算得到的y_scale和回歸項的均值和標(biāo)準(zhǔn)差。

點預(yù)測

點預(yù)測也就是在時間點上輸出單個預(yù)測值，預(yù)測值由趨勢預(yù)測和回歸項預(yù)測組成。

• 趨勢預(yù)測

python.Prophet.forecaster.Prophet.predict_trend

在貝葉斯回歸中，未知參數(shù)服從一個指定的先驗分布，Prophet使用Stan引擎計算得到的返回參數(shù)的期望作為趨勢項公式的帶入值。

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']

復(fù)制

對于flat、linear、logistic不同的趨勢分量，結(jié)合changepoints_t進行趨勢的預(yù)測，以線性趨勢為例。

python.Prophet.forecaster.Prophet.piecewise_linear

即論文中描述的公式：

最后根據(jù)訓(xùn)練過程計算的y_scale縮放變量還原量級，作為最終的趨勢預(yù)測。

• 回歸項預(yù)測

python.Prophet.forecaster.Prophet.predict_seasonal_components

上文提到的Prophet的seasonality features、holiday features和add regressors項，在參數(shù)求解過程中是一樣處理方法。

β即為擬合出的權(quán)重系數(shù)，β乘以特征寬表對應(yīng)的值，得到每個特征的效應(yīng)值。根據(jù)component_cols維護的不同成分的特征歸屬。就可以得到seasonality、holiday features、add regressors、multiplicative_terms、additive_terms各個部分的效應(yīng)值。這些在模型解釋性、結(jié)果可視化、偏差分析中都可以應(yīng)用。

• 預(yù)測匯總

python/Prophet/forecaster.py:1287

預(yù)測成分匯總步驟把趨勢、季節(jié)性、節(jié)假日等預(yù)測效應(yīng)值進行拼接，最后的預(yù)測結(jié)果公式如下：

預(yù)測結(jié)果值=趨勢*(1+乘法因素)+加法因素

預(yù)測的不確定性

Prophet通過觀測噪聲、參數(shù)不確定性以及未來趨勢不確定性三個部分輸出預(yù)測值的不確定性，最終通過yhat_lower、yhat_upper體現(xiàn)。

• 觀測噪聲不確定性

python/Prophet/forecaster.py:1556

一般時間序列分解模型，最后的一個分解項就是噪聲，Prophet也是這樣的，假設(shè)了輸入數(shù)據(jù)點含有服從正態(tài)分布的觀測噪聲，Stan模型腳本中的sigma_obs就是擬合了這部分噪聲。觀測噪聲不確定性為：

sigma = self.params['sigma_obs'][iteration]
noise_terms = np.random.normal(0, sigma, trends.shape) * self.y_scale

復(fù)制

變量sigma是sigma_obs的估計值。在每個時間步驟上，從正態(tài)分布中采樣n_samples次，得到噪聲不確定項目noise_terms。

• 趨勢不確定性

python.Prophet.forecaster.Prophet._make_trend_shift_matrix

對于趨勢的不確定性，官方文檔中解釋：預(yù)測中最大的不確定性來源在于未來趨勢變化的可能性。因此我們盡最大可能采取了一個合理的假設(shè)，假設(shè)未來的趨勢變化會與歷史上觀察到的相似。特別是，我們假設(shè)未來的趨勢變化的平均頻率和幅度將與歷史上觀察到的相同。我們將這些趨勢變化進行投影，并通過計算它們的分布來獲得不確定性區(qū)間。

具體而言，趨勢不確性由突變點出現(xiàn)的位置和突變的比例確定。首先計算歷史上突變點出現(xiàn)的間隔的均值*歷史上突變點的個數(shù)，得到每個時間點上產(chǎn)生突變點概率likelihood。然后計算歷史突變點的變化率的絕對值的均值mean_delta。最后進行采樣來模擬不確定性。均勻分布采樣輸出矩陣（n_samples, future_length），對每一個元素判斷是否小于似然值likelihood，來判斷是否發(fā)生了突變，如果小于則發(fā)生了突變，然后輸出均值為mean_delta的拉普拉斯分布，得到突變點變化率的采樣值。

• 參數(shù)的不確定性

python.Prophet.forecaster.Prophet.predict_uncertainty

當(dāng)初始化參數(shù)MCMC=0，即調(diào)用Stan用MAP算法，此時算法得出的是收斂時的最終估計值，不確定性由觀測不確定性和趨勢不確定性組成，參數(shù)沒有不確定性的估計。

MCMC>0時候，除了觀測噪聲、趨勢不確定性外，由于Stan采用MCMC采樣方法，得到的返回值是待估計參數(shù)值的若干個采樣點，根據(jù)每個參數(shù)的采樣點得到不同的預(yù)測值。

• 不確定性結(jié)果輸出

進一步，由于每一個時間步有以上不確定性成分的n_samples個采樣數(shù)據(jù)。根據(jù)預(yù)設(shè)的置信度寬度，得到上、下分位數(shù)值：

lower_p = 100 * (1.0 - self.interval_width) / 2
upper_p = 100 * (1.0 + self.interval_width) / 2

復(fù)制

對于每個時間步驟，求lower_p分位數(shù)、upper_p分位數(shù)，即得到y(tǒng)hat_lower、yhat_upper。

四、總結(jié)

Prophet框架在有明顯規(guī)律的單變量時序預(yù)測場景中有著非常不錯的表現(xiàn)。在得物的多個業(yè)務(wù)場景中已經(jīng)得到了驗證。閱讀源碼能夠幫助我們更好的了解模型細節(jié)，在模型優(yōu)化時幫助我們有的放矢，在模型選型時幫助我們掌握模型的適應(yīng)性和局限性。

參考資料

??https://peerj.com/preprints/3190/??

??https://github.com/facebook/prophet??

??https://zhuanlan.zhihu.com/p/37543542??

本文轉(zhuǎn)載自 ??得物技術(shù)??，作者： 探流