Facebook 時間序列預測演算法 Prophet 的研究

Prophet 簡介

Facebook 去年開源了一個時間序列預測的演算法，叫做 fbprophet，它的官方網址與基本介紹來自於以下幾個網站：

  1. Github：github.com/facebook/pro
  2. 官方網址：facebook.github.io/prop
  3. 論文名字與網址：Forecasting at scale，peerj.com/preprints/319

從官網的介紹來看，Facebook 所提供的 prophet 演算法不僅可以處理時間序列存在一些異常值的情況，也可以處理部分缺失值的情形，還能夠幾乎全自動地預測時間序列未來的走勢。從論文上的描述來看，這個 prophet 演算法是基於時間序列分解和機器學習的擬合來做的，其中在擬合模型的時候使用了 pyStan 這個開源工具，因此能夠在較快的時間內得到需要預測的結果。除此之外，為了方便統計學家，機器學習從業者等人羣的使用，prophet 同時提供了 R 語言和 Python 語言的介面。從整體的介紹來看，如果是一般的商業分析或者數據分析的需求，都可以嘗試使用這個開源演算法來預測未來時間序列的走勢。

Prophet 的演算法原理

Prophet 數據的輸入和輸出

首先讓我們來看一個常見的時間序列場景，黑色表示原始的時間序列離散點，深藍色的線表示使用時間序列來擬合所得到的取值，而淺藍色的線表示時間序列的一個置信區間，也就是所謂的合理的上界和下界。prophet 所做的事情就是：

  1. 輸入已知的時間序列的時間戳和相應的值；
  2. 輸入需要預測的時間序列的長度；
  3. 輸出未來的時間序列走勢。
  4. 輸出結果可以提供必要的統計指標，包括擬合曲線，上界和下界等。

就一般情況而言，時間序列的離線存儲格式為時間戳和值這種格式，更多的話可以提供時間序列的 ID，標籤等內容。因此，離線存儲的時間序列通常都是以下的形式。其中 date 指的是具體的時間戳，category 指的是某條特定的時間序列 id，value 指的是在 date 下這個 category 時間序列的取值，label 指的是人工標記的標籤（0 表示異常，1『 表示正常，unknown 表示沒有標記或者人工判斷不清）。

而 fbprophet 所需要的時間序列也是這種格式的，根據官網的描述，只要用 csv 文件存儲兩列即可，第一列的名字是 ds, 第二列的名稱是 y。第一列表示時間序列的時間戳，第二列表示時間序列的取值。通過 prophet 的計算，可以計算出 yhat，yhat_lower，yhat_upper，分別表示時間序列的預測值，預測值的下界，預測值的上界。兩份表格如下面的兩幅圖表示。

Prophet 的演算法實現

在時間序列分析領域，有一種常見的分析方法叫做時間序列的分解（Decomposition of Time Series），它把時間序列 y_{t} 分成幾個部分，分別是季節項 S_{t} ，趨勢項 T_{t} ，剩餘項 R_{t} 。也就是說對所有的 tgeq 0 ，都有

y_{t} = S_{t} + T_{t} + R_{t}.

除了加法的形式，還有乘法的形式，也就是：

y_{t} = S_{t} imes T_{t} imes R_{t}.

以上式子等價於 ln y_{t} = ln S_{t} + ln T_{t} + ln R_{t} 。所以，有的時候在預測模型的時候，會先取對數，然後再進行時間序列的分解，就能得到乘法的形式。在 fbprophet 演算法中，作者們基於這種方法進行了必要的改進和優化。

一般來說，在實際生活和生產環節中，除了季節項，趨勢項，剩餘項之外，通常還有節假日的效應。所以，在 prophet 演算法裡面，作者同時考慮了以上四項，也就是：

y(t) = g(t) + s(t) + h(t) + epsilon_{t}.

其中 g(t) 表示趨勢項，它表示時間序列在非週期上面的變化趨勢； s(t) 表示週期項，或者稱為季節項，一般來說是以周或者年為單位； h(t) 表示節假日項，表示在當天是否存在節假日；epsilon_{t} 表示誤差項或者稱為剩餘項。Prophet 演算法就是通過擬合這幾項，然後最後把它們累加起來就得到了時間序列的預測值。

趨勢項模型 g(t)

在 Prophet 演算法裡面，趨勢項有兩個重要的函數，一個是基於邏輯回歸函數（logistic function）的，另一個是基於分段線性函數（piecewise linear function）的。

首先，我們來介紹一下基於邏輯回歸的趨勢項是怎麼做的。

如果回顧邏輯回歸函數的話，一般都會想起這樣的形式： sigma(x) = 1/(1+e^{-x}), 它的導數是 sigma(x) = sigma(x) cdot(1-sigma(x)), 並且 lim_{x ightarrow +infty} sigma(x) = 1, lim_{x ightarrow -infty} sigma(x) = 0. 如果增加一些參數的話，那麼邏輯回歸就可以改寫成：

f(x) = C / (1 + e^{-k(x-m)}),

這裡的 C 稱為曲線的最大漸近值， k 表示曲線的增長率， m 表示曲線的中點。當 C=1, k = 1, m =0 時，恰好就是大家常見的 sigmoid 函數的形式。從 sigmoid 的函數表達式來看，它滿足以下的微分方程： y=y(1-y) 。

那麼，如果使用分離變數法來求解微分方程 y=y(1-y) 就可以得到：

frac{y}{y} + frac{y}{1-y} = 1 Rightarrow lnfrac{y}{1-y} = 1 Rightarrow y = 1/(1+K e^{-x}).

但是在現實環境中，函數 f(x) = C / (1+e^{-k(x-m)}) 的三個參數 C, k, m 不可能都是常數，而很有可能是隨著時間的遷移而變化的，因此，在 Prophet 裡面，作者考慮把這三個參數全部換成了隨著時間而變化的函數，也就是 C = C(t), k = k(t), m = m(t).

除此之外，在現實的時間序列中，曲線的走勢肯定不會一直保持不變，在某些特定的時候或者有著某種潛在的週期曲線會發生變化，這種時候，就有學者會去研究變點檢測，也就是所謂 change point detection。例如下面的這幅圖的 t_{1}^{*}, t_{2}^{*} 就是時間序列的兩個變點。

在 Prophet 裡面，是需要設置變點的位置的，而每一段的趨勢和走勢也是會根據變點的情況而改變的。在程序裡面有兩種方法，一種是通過人工指定的方式指定變點的位置；另外一種是通過演算法來自動選擇。在默認的函數裡面，Prophet 會選擇 n_changepoints = 25 個變點，然後設置變點的範圍是前 80%（changepoint_range），也就是在時間序列的前 80% 的區間內會設置變點。通過 forecaster.py 裡面的 set_changepoints 函數可以知道，首先要看一些邊界條件是否合理，例如時間序列的點數是否少於 n_changepoints 等內容；其次如果邊界條件符合，那變點的位置就是均勻分佈的，這一點可以通過 np.linspace 這個函數看出來。

下面假設已經放置了 S 個變點了，並且變點的位置是在時間戳 s_{j}, 1leq jleq S 上，那麼在這些時間戳上，我們就需要給出增長率的變化，也就是在時間戳 s_{j} 上發生的 change in rate。可以假設有這樣一個向量： oldsymbol{delta}inmathbb{R}^{S}, 其中 delta_{j} 表示在時間戳 s_{j} 上的增長率的變化量。如果一開始的增長率我們使用 k 來代替的話，那麼在時間戳 t 上的增長率就是 k + sum_{j:t>s_{j}} delta_{j} ，通過一個指示函數 mathbf{a}(t)in {0,1}^{S} 就是

a_{j}(t) = egin{cases} 1, ext{ if } tgeq s_{j},\ 0, ext{ otherwise.} end{cases}

那麼在時間戳 t 上面的增長率就是 k + mathbf{a}^{T}oldsymbol{delta}. 一旦變化量 k 確定了，另外一個參數 m 也要隨之確定。在這裡需要把線段的邊界處理好，因此通過數學計算可以得到：

gamma_{j} = igg(s_{j} - m - sum_{ell <j} gamma_{ell} igg) cdot igg( 1- frac{k + sum_{ell < j} delta_{ell}}{k + sum_{ellleq j}delta_{ell}} igg).

所以，分段的邏輯回歸增長模型就是：

g(t) = frac{C(t)}{1+exp(-(k+oldsymbol{a}(t)^{t}oldsymbol{delta}) cdot (t - (m+oldsymbol{a}(t)^{T}oldsymbol{gamma})},

其中，

oldsymbol{a}(t) = (a_{1}(t),cdots,a_{S}(t))^{T}, oldsymbol{delta} = (delta_{1},cdots,delta_{S})^{T}, oldsymbol{gamma} = (gamma_{1},cdots,gamma_{S})^{T}.

在邏輯回歸函數裡面，有一個參數是需要提前設置的，那就是 Capacity，也就是所謂的 C(t) ，在使用 Prophet 的 growth = logistic 的時候，需要提前設置好 C(t) 的取值纔行。

再次，我們來介紹一下基於分段線性函數的趨勢項是怎麼做的。眾所周知，線性函數指的是 y=kx+b, 而分段線性函數指的是在每一個子區間上，函數都是線性函數，但是在整段區間上，函數並不完全是線性的。正如下圖所示，分段線性函數就是一個折線的形狀。

因此，基於分段線性函數的模型形如：

g(t)=(k+oldsymbol{a}(t)oldsymbol{delta})cdot t+(m+oldsymbol{a}(t)^{T}oldsymbol{gamma}),

其中 k 表示增長率（growth rate）， oldsymbol{delta} 表示增長率的變化量， m 表示 offset parameter。而這兩種方法（分段線性函數與邏輯回歸函數）最大的區別就是 oldsymbol{gamma} 的設置不一樣，在分段線性函數中， oldsymbol{gamma}=(gamma_{1},cdots,gamma_{S})^{T}, gamma_{j}=-s_{j}delta_{j}. 注意：這與之前邏輯回歸函數中的設置是不一樣的。

在 Prophet 的源代碼中，forecast.py 這個函數裡麪包含了最關鍵的步驟，其中 piecewise_logistic 函數表示了前面所說的基於邏輯回歸的增長函數，它的輸入包含了 cap 這個指標，因此需要用戶事先指定 capacity。而在 piecewise_linear 這個函數中，是不需要 capacity 這個指標的，因此 m = Prophet() 這個函數默認的使用 growth = linear 這個增長函數，也可以寫作 m = Prophet(growth = linear)；如果想用 growth = logistic，就要這樣寫：

m = Prophet(growth=logistic)
df[cap] = 6
m.fit(df)
future = m.make_future_dataframe(periods=prediction_length, freq=min)
future[cap] = 6

變點的選擇（Changepoint Selection）

在介紹變點之前，先要介紹一下 Laplace 分佈，它的概率密度函數為：

f(x|mu, b) = expigg(-|x-mu|/bigg)/2b,

其中 mu 表示位置參數， b>0 表示尺度參數。Laplace 分佈與正態分佈有一定的差異。

在 Prophet 演算法中，是需要給出變點的位置，個數，以及增長的變化率的。因此，有三個比較重要的指標，那就是

  1. changepoint_range，
  2. n_changepoint，
  3. changepoint_prior_scale。

changepoint_range 指的是百分比，需要在前 changepoint_range 那麼長的時間序列中設置變點，在默認的函數中是 changepoint_range = 0.8。n_changepoint 表示變點的個數，在默認的函數中是 n_changepoint = 25。changepoint_prior_scale 表示變點增長率的分佈情況，在論文中， delta_{j} sim Laplace(0, au) ，這裡的 au 就是 change_point_scale。

在整個開源框架裡面，在默認的場景下，變點的選擇是基於時間序列的前 80% 的歷史數據，然後通過等分的方法找到 25 個變點（change points），而變點的增長率是滿足 Laplace 分佈 delta_{j} sim Laplace (0,0.05) 的。因此，當 au 趨近於零的時候， delta_{j} 也是趨向於零的，此時的增長函數將變成全段的邏輯回歸函數或者線性函數。這一點從 g(t) 的定義可以輕易地看出。

對未來的預估（Trend Forecast Uncertainty）

從歷史上長度為 T 的數據中，我們可以選擇出 S 個變點，它們所對應的增長率的變化量是 delta_{j} sim Laplace(0, au) 。此時我們需要預測未來，因此也需要設置相應的變點的位置，從代碼中看，在 forecaster.py 的 sample_predictive_trend 函數中，通過 Poisson 分佈等概率分佈方法找到新增的 changepoint_ts_new 的位置，然後與 changepoint_t 拼接在一起就得到了整段序列的 changepoint_ts。

changepoint_ts_new = 1 + np.random.rand(n_changes) * (T - 1)
changepoint_ts = np.concatenate((self.changepoints_t, changepoint_ts_new))

第一行代碼的 1 保證了 changepoint_ts_new 裡面的元素都大於 change_ts 裡面的元素。除了變點的位置之外，也需要考慮 delta 的情況。這裡令 lambda = sum_{j=1}^{S}|delta_{j}|/S ，於是新的增長率的變化量就是按照下面的規則來選擇的：當 j>T 時，

delta_{j}=egin{cases} 0 ext{, with probability } (T-S)/T \ sim Laplace(0,lambda) ext{, with probability } S/T end{cases}.

季節性趨勢

幾乎所有的時間序列預測模型都會考慮這個因素，因為時間序列通常會隨著天，周，月，年等季節性的變化而呈現季節性的變化，也稱為週期性的變化。對於週期函數而言，大家能夠馬上聯想到的就是正弦餘弦函數。而在數學分析中，區間內的週期性函數是可以通過正弦和餘弦的函數來表示的：假設 f(x) 是以 2pi 為週期的函數，那麼它的傅立葉級數就是 a_{0} + sum_{n=1}^{infty}(a_{n}cos(nx) + b_{n}sin(nx)) 。

在論文中，作者使用傅立葉級數來模擬時間序列的週期性。假設 P 表示時間序列的週期， P = 365.25 表示以年為週期， P = 7 表示以周為週期。它的傅立葉級數的形式都是：

s(t) = sum_{n=1}^{N}igg( a_{n}cosigg(frac{2pi n t}{P}igg) + b_{n}sinigg(frac{2pi n t}{P}igg)igg).

就作者的經驗而言，對於以年為週期的序列（ P = 365.25 ）而言， N = 10 ；對於以周為週期的序列（P = 7 ）而言， N = 3 。這裡的參數可以形成列向量：

oldsymbol{eta} = (a_{1},b_{1},cdots,a_{N},b_{N})^{T}.

當 N = 10 時，

X(t) = igg[cosigg(frac{2pi(1)t}{365.25}igg),cdots,sinigg(frac{2pi(10)t}{365.25}igg)igg]

當 N = 3 時，

X(t) = igg[cosigg(frac{2pi(1)t}{7}igg),cdots,sinigg(frac{2pi(3)t}{7}igg)igg].

因此，時間序列的季節項就是： s(t) = X(t) oldsymbol{eta}, 而 oldsymbol{eta} 的初始化是 oldsymbol{eta} sim Normal(0,sigma^{2})。這裡的 sigma 是通過 seasonality_prior_scale 來控制的，也就是說 sigma= seasonality_prior_scale。這個值越大，表示季節的效應越明顯；這個值越小，表示季節的效應越不明顯。同時，在代碼裡面，seasonality_mode 也對應著兩種模式，分別是加法和乘法，默認是加法的形式。在開源代碼中， X(t) 函數是通過 fourier_series 來構建的。

節假日效應（holidays and events）

在現實環境中，除了週末，同樣有很多節假日，而且不同的國家有著不同的假期。在 Prophet 裡面，通過維基百科裡面對各個國家的節假日的描述，hdays.py 收集了各個國家的特殊節假日。除了節假日之外，用戶還可以根據自身的情況來設置必要的假期，例如 The Super Bowl，雙十一等。

由於每個節假日對時間序列的影響程度不一樣，例如春節，國慶節則是七天的假期，對於勞動節等假期來說則假日較短。因此，不同的節假日可以看成相互獨立的模型，並且可以為不同的節假日設置不同的前後窗口值，表示該節假日會影響前後一段時間的時間序列。用數學語言來說，對與第 i 個節假日來說， D_{i} 表示該節假日的前後一段時間。為了表示節假日效應，我們需要一個相應的指示函數（indicator function），同時需要一個參數 kappa_{i} 來表示節假日的影響範圍。假設我們有 L 個節假日，那麼

h(t)=Z(t) oldsymbol{kappa}=sum_{i=1}^{L} kappa_{i}cdot 1_{{tin D_{i}}},

其中 Z(t)=(1_{{tin D_{1}}},cdots,1_{{tin D_{L}}}), ext{ } oldsymbol{kappa}=(kappa_{1},cdots,kappa_{L})^{T}.

其中 oldsymbol{kappa}sim Normal(0,v^{2}) 並且該正態分佈是受到 v = holidays_prior_scale 這個指標影響的。默認值是 10，當值越大時，表示節假日對模型的影響越大；當值越小時，表示節假日對模型的效果越小。用戶可以根據自己的情況自行調整。

模型擬合（Model Fitting）

按照以上的解釋，我們的時間序列已經可以通過增長項，季節項，節假日項來構建了，i.e.

y(t)=g(t)+s(t)+h(t)+epsilon

下一步我們只需要擬合函數就可以了，在 Prophet 裡面，作者使用了 pyStan 這個開源工具中的 L-BFGS 方法來進行函數的擬合。具體可以參考 forecast.py 裡面的 stan_init 函數。

Prophet 中可以設置的參數

在 Prophet 中，用戶一般可以設置以下四種參數：

  1. Capacity：在增量函數是邏輯回歸函數的時候，需要設置的容量值。
  2. Change Points：可以通過 n_changepoints 和 changepoint_range 來進行等距的變點設置，也可以通過人工設置的方式來指定時間序列的變點。
  3. 季節性和節假日：可以根據實際的業務需求來指定相應的節假日。
  4. 光滑參數： au= changepoint_prior_scale 可以用來控制趨勢的靈活度， sigma= seasonality_prior_scale 用來控制季節項的靈活度， v= holidays prior scale 用來控制節假日的靈活度。

如果不想設置的話，使用 Prophet 默認的參數即可。

Prophet 的實際使用

Prophet 的簡單使用

因為 Prophet 所需要的兩列名稱是 ds 和 y，其中，ds 表示時間戳，y 表示時間序列的值，因此通常來說都需要修改 pd.dataframe 的列名字。如果原來的兩列名字是 timestamp 和 value 的話，只需要這樣寫：

df = df.rename(columns={timestamp:ds, value:y})

如果 timestamp 是使用 unixtime 來記錄的，需要修改成 YYYY-MM-DD hh:mm:ss 的形式：

df[ds] = pd.to_datetime(df[ds],unit=s)

在一般情況下，時間序列需要進行歸一化的操作，而 pd.dataframe 的歸一化操作也十分簡單：

df[y] = (df[y] - df[y].mean()) / (df[y].std())

然後就可以初始化模型，然後擬合模型，並且進行時間序列的預測了。

初始化模型：m = Prophet()
擬合模型：m.fit(df)
計算預測值：periods 表示需要預測的點數，freq 表示時間序列的頻率。
future = m.make_future_dataframe(periods=30, freq=min)
future.tail()
forecast = m.predict(future)

而 freq 指的是 pd.dataframe 裡面的一個指標，min 表示按分鐘來收集的時間序列。具體參見文檔：pandas.pydata.org/panda

在進行了預測操作之後，通常都希望把時間序列的預測趨勢畫出來：

畫出預測圖：
m.plot(forecast)
畫出時間序列的分量：
m.plot_components(forecast)

如果要畫出更詳細的指標，例如中間線，上下界，那麼可以這樣寫：

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()

其實 Prophet 預測的結果都放在了變數 forecast 裡面，列印結果的話可以這樣寫：第一行是列印所有時間戳的預測結果，第二行是列印最後五個時間戳的預測結果。

print(forecast[[ds, yhat, yhat_lower, yhat_upper]])
print(forecast[[ds, yhat, yhat_lower, yhat_upper]].tail())

Prophet 的參數設置

Prophet 的默認參數可以在 forecaster.py 中看到：

def __init__(
self,
growth=linear,
changepoints=None,
n_changepoints=25,
changepoint_range=0.8,
yearly_seasonality=auto,
weekly_seasonality=auto,
daily_seasonality=auto,
holidays=None,
seasonality_mode=additive,
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,
):

增長函數的設置

在 Prophet 裡面，有兩個增長函數，分別是分段線性函數（linear）和邏輯回歸函數（logistic）。而 m = Prophet() 默認使用的是分段線性函數（linear），並且如果要是用邏輯回歸函數的時候，需要設置 capacity 的值，i.e. df[cap] = 100，否則會出錯。

m = Prophet()
m = Prophet(growth=linear)
m = Prophet(growth=logistic)

變點的設置

在 Prophet 裡面，變點默認的選擇方法是前 80% 的點中等距選擇 25 個點作為變點，也可以通過以下方法來自行設置變點，甚至可以人為設置某些點。

m = Prophet(n_changepoints=25)
m = Prophet(changepoint_range=0.8)
m = Prophet(changepoint_prior_scale=0.05)
m = Prophet(changepoints=[2014-01-01])

而變點的作圖可以使用：

from fbprophet.plot import add_changepoints_to_plot
fig = m.plot(forecast)
a = add_changepoints_to_plot(fig.gca(), m, forecast)

週期性的設置

通常來說，可以在 Prophet 裡面設置週期性，無論是按月還是周其實都是可以設置的，例如：

m = Prophet(weekly_seasonality=False)
m.add_seasonality(name=monthly, period=30.5, fourier_order=5)
m = Prophet(weekly_seasonality=True)
m.add_seasonality(name=weekly, period=7, fourier_order=3, prior_scale=0.1)

節假日的設置

有的時候，由於雙十一或者一些特殊節假日，我們可以設置某些天數是節假日，並且設置它的前後影響範圍，也就是 lower_window 和 upper_window。

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

m = Prophet(holidays=holidays, holidays_prior_scale=10.0)

結束語

對於商業分析等領域的時間序列，Prophet 可以進行很好的擬合和預測，但是對於一些週期性或者趨勢性不是很強的時間序列，用 Prophet 可能就不合適了。但是，Prophet 提供了一種時序預測的方法，在用戶不是很懂時間序列的前提下都可以使用這個工具得到一個能接受的結果。具體是否用 Prophet 則需要根據具體的時間序列來確定。

參考文獻：

  1. otexts.org/fpp2/compone
  2. en.wikipedia.org/wiki/D
  3. A review of change point detection methods, CTruong, L. Oudre, N.Vayatis
  4. github.com/facebook/pro
  5. facebook.github.io/prop

推薦閱讀：

相關文章