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

一、引入

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

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

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

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

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

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

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

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

二、準(zhǔn)備知識(shí)

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

參數(shù)估計(jì)

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

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

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

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

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

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

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

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

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

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

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

Stan開源框架

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

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

Prophet基本原理

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

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

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

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

• 變量和數(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
  );
}

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

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

三、Prophet源碼剖析

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

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

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

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

模型框架

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

初始參數(shù)

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

趨勢(shì)項(xiàng)相關(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)練過(guò)程

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

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

python.Prophet.forecaster.Prophet.setup_dataframe

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

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

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

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

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

python.Prophet.forecaster.Prophet.initialize_scales

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

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

MinMax歸一化：
含義：MinMax歸一化是將原始數(shù)據(jù)縮放到[0, 1]的范圍內(nèi)，使數(shù)據(jù)的最小值對(duì)應(yīng)0，最大值對(duì)應(yīng)1。
適用場(chǎng)景：MinMax歸一化適用于需要將數(shù)據(jù)映射到一個(gè)固定范圍內(nèi)的情況，保留了原始數(shù)據(jù)的相對(duì)關(guān)系，在大部分機(jī)器學(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)
    ...
    ...

對(duì)于add regressor特征則進(jìn)行標(biāo)準(zhǔn)化，將特征數(shù)據(jù)縮放到均值為0，標(biāo)準(zhǔn)差為1的標(biāo)準(zhǔn)正態(tài)分布。

python/Prophet/forecaster.py:389

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

自動(dòng)設(shè)置周期性

python.Prophet.forecaster.Prophet.set_auto_seasonalities

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

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

生成特征寬表

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

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

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

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

python.Prophet.forecaster.Prophet.set_changepoints

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

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

初始化stan參數(shù)

python.Prophet.forecaster.Prophet.calculate_initial_params

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

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

參數(shù)擬合

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

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

• 牛頓法/LBFGS法進(jìn)行最大后驗(yàn)估計(jì)

python.Prophet.models.CmdStanPyBackend.fit

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

python.Prophet.models.CmdStanPyBackend.sampling

預(yù)測(cè)過(guò)程

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

python.Prophet.forecaster.Prophet.make_future_dataframe

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

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

點(diǎn)預(yù)測(cè)

點(diǎn)預(yù)測(cè)也就是在時(shí)間點(diǎn)上輸出單個(gè)預(yù)測(cè)值，預(yù)測(cè)值由趨勢(shì)預(yù)測(cè)和回歸項(xiàng)預(yù)測(cè)組成。

• 趨勢(shì)預(yù)測(cè)

python.Prophet.forecaster.Prophet.predict_trend

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

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ù)制

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

python.Prophet.forecaster.Prophet.piecewise_linear

即論文中描述的公式：

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

• 回歸項(xiàng)預(yù)測(cè)

python.Prophet.forecaster.Prophet.predict_seasonal_components

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

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

• 預(yù)測(cè)匯總

python/Prophet/forecaster.py:1287

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

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

預(yù)測(cè)的不確定性

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

• 觀測(cè)噪聲不確定性

python/Prophet/forecaster.py:1556

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

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

復(fù)制

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

• 趨勢(shì)不確定性

python.Prophet.forecaster.Prophet._make_trend_shift_matrix

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

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

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

python.Prophet.forecaster.Prophet.predict_uncertainty

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

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

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

進(jìn)一步，由于每一個(gè)時(shí)間步有以上不確定性成分的n_samples個(gè)采樣數(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ù)制

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

四、總結(jié)

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

參考資料

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

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

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

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