近期整理了一下 Facebook 的 Prophet，個人感覺這是一個非常不錯的時間序列預測工具。
Prophet 簡介
Facebook 去年開源了一個時間序列預測的算法，叫做 fbprophet，它的官方網址與基本介紹來自于以下幾個網站：
從官網的介紹來看，Facebook 所提供的 prophet 算法不僅可以處理時間序列存在一些異常值的情況，也可以處理部分缺失值的情形，還能夠幾乎全自動地預測時間序列未來的走勢。從論文上的描述來看，這個 prophet 算法是基于時間序列分解和機器學習的拟合來做的，其中在拟合模型的時候使用了 pyStan 這個開源工具，是以能夠在較快的時間内得到需要預測的結果。除此之外，為了友善統計學家，機器學習從業者等人群的使用，prophet 同時提供了 R 語言和 Python 語言的接口。從整體的介紹來看，如果是一般的商業分析或者資料分析的需求，都可以嘗試使用這個開源算法來預測未來時間序列的走勢。
Prophet 的算法原理
Prophet 資料的輸入和輸出
首先讓我們來看一個常見的時間序列場景，黑色表示原始的時間序列離散點，深藍色的線表示使用時間序列來拟合所得到的取值，而淺藍色的線表示時間序列的一個置信區間，也就是所謂的合理的上界和下界。prophet 所做的事情就是：輸入已知的時間序列的時間戳和相應的值；
輸入需要預測的時間序列的長度；
輸出未來的時間序列走勢。
輸出結果可以提供必要的統計名額，包括拟合曲線，上界和下界等。
就一般情況而言，時間序列的離線存儲格式為時間戳和值這種格式，更多的話可以提供時間序列的 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)，它把時間序列
分成幾個部分，分别是季節項
，趨勢項
，剩餘項
。也就是說對所有的
，都有
除了加法的形式，還有乘法的形式，也就是：
以上式子等價于
。是以，有的時候在預測模型的時候，會先取對數，然後再進行時間序列的分解，就能得到乘法的形式。在 fbprophet 算法中，作者們基于這種方法進行了必要的改進和優化。
一般來說，在實際生活和生産環節中，除了季節項，趨勢項，剩餘項之外，通常還有節假日的效應。是以，在 prophet 算法裡面，作者同時考慮了以上四項，也就是：
其中
表示趨勢項，它表示時間序列在非周期上面的變化趨勢；
表示周期項，或者稱為季節項，一般來說是以周或者年為機關；
表示節假日項，表示在當天是否存在節假日；
表示誤差項或者稱為剩餘項。Prophet 算法就是通過拟合這幾項，然後最後把它們累加起來就得到了時間序列的預測值。
趨勢項模型
在 Prophet 算法裡面，趨勢項有兩個重要的函數，一個是基于邏輯回歸函數(logistic function)的，另一個是基于分段線性函數(piecewise linear function)的。
首先，我們來介紹一下基于邏輯回歸的趨勢項是怎麼做的。
如果回顧邏輯回歸函數的話，一般都會想起這樣的形式：
它的導數是
并且
如果增加一些參數的話，那麼邏輯回歸就可以改寫成：
這裡的
稱為曲線的最大漸近值，
表示曲線的增長率，
表示曲線的中點。當
時，恰好就是大家常見的 sigmoid 函數的形式。從 sigmoid 的函數表達式來看，它滿足以下的微分方程：
。
那麼，如果使用分離變量法來求解微分方程
就可以得到：
但是在現實環境中，函數
的三個參數
不可能都是常數，而很有可能是随着時間的遷移而變化的，是以，在 Prophet 裡面，作者考慮把這三個參數全部換成了随着時間而變化的函數，也就是
除此之外，在現實的時間序列中，曲線的走勢肯定不會一直保持不變，在某些特定的時候或者有着某種潛在的周期曲線會發生變化，這種時候，就有學者會去研究變點檢測，也就是所謂 change point detection。例如下面的這幅圖的
就是時間序列的兩個變點。
在 Prophet 裡面，是需要設定變點的位置的，而每一段的趨勢和走勢也是會根據變點的情況而改變的。在程式裡面有兩種方法，一種是通過人工指定的方式指定變點的位置；另外一種是通過算法來自動選擇。在預設的函數裡面，Prophet 會選擇 n_changepoints = 25 個變點，然後設定變點的範圍是前 80%(changepoint_range)，也就是在時間序列的前 80% 的區間内會設定變點。通過 forecaster.py 裡面的 set_changepoints 函數可以知道，首先要看一些邊界條件是否合理，例如時間序列的點數是否少于 n_changepoints 等内容；其次如果邊界條件符合，那變點的位置就是均勻分布的，這一點可以通過 np.linspace 這個函數看出來。
下面假設已經放置了
個變點了，并且變點的位置是在時間戳
上，那麼在這些時間戳上，我們就需要給出增長率的變化，也就是在時間戳
上發生的 change in rate。可以假設有這樣一個向量：
其中
表示在時間戳
上的增長率的變化量。如果一開始的增長率我們使用
來代替的話，那麼在時間戳
上的增長率就是
，通過一個訓示函數
就是
那麼在時間戳
上面的增長率就是
一旦變化量
确定了，另外一個參數
也要随之确定。在這裡需要把線段的邊界處理好，是以通過數學計算可以得到：
是以，分段的邏輯回歸增長模型就是：
其中，
在邏輯回歸函數裡面，有一個參數是需要提前設定的，那就是 Capacity，也就是所謂的
，在使用 Prophet 的 growth = 'logistic' 的時候，需要提前設定好
的取值才行。
再次，我們來介紹一下基于分段線性函數的趨勢項是怎麼做的。衆所周知，線性函數指的是
而分段線性函數指的是在每一個子區間上，函數都是線性函數，但是在整段區間上，函數并不完全是線性的。正如下圖所示，分段線性函數就是一個折線的形狀。
是以，基于分段線性函數的模型形如：
其中
表示增長率(growth rate)，
表示增長率的變化量，
表示 offset parameter。而這兩種方法(分段線性函數與邏輯回歸函數)最大的差別就是
的設定不一樣，在分段線性函數中，
注意：這與之前邏輯回歸函數中的設定是不一樣的。
在 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 分布，它的機率密度函數為：
其中
表示位置參數，
表示尺度參數。Laplace 分布與正态分布有一定的差異。
在 Prophet 算法中，是需要給出變點的位置，個數，以及增長的變化率的。是以，有三個比較重要的名額，那就是changepoint_range，
n_changepoint，
changepoint_prior_scale。
changepoint_range 指的是百分比，需要在前 changepoint_range 那麼長的時間序列中設定變點，在預設的函數中是 changepoint_range = 0.8。n_changepoint 表示變點的個數，在預設的函數中是 n_changepoint = 25。changepoint_prior_scale 表示變點增長率的分布情況，在論文中，
，這裡的
就是 change_point_scale。
在整個開源架構裡面，在預設的場景下，變點的選擇是基于時間序列的前 80% 的曆史資料，然後通過等分的方法找到 25 個變點(change points)，而變點的增長率是滿足 Laplace 分布
的。是以，當
趨近于零的時候，
也是趨向于零的，此時的增長函數将變成全段的邏輯回歸函數或者線性函數。這一點從
的定義可以輕易地看出。
對未來的預估(Trend Forecast Uncertainty)
從曆史上長度為
的資料中，我們可以選擇出
個變點，它們所對應的增長率的變化量是
。此時我們需要預測未來，是以也需要設定相應的變點的位置，從代碼中看，在 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 裡面的元素。除了變點的位置之外，也需要考慮
的情況。這裡令
，于是新的增長率的變化量就是按照下面的規則來選擇的：當
時，
季節性趨勢
幾乎所有的時間序列預測模型都會考慮這個因素，因為時間序列通常會随着天，周，月，年等季節性的變化而呈現季節性的變化，也稱為周期性的變化。對于周期函數而言，大家能夠馬上聯想到的就是正弦餘弦函數。而在數學分析中，區間内的周期性函數是可以通過正弦和餘弦的函數來表示的：假設
是以
為周期的函數，那麼它的傅立葉級數就是
。
在論文中，作者使用傅立葉級數來模拟時間序列的周期性。假設
表示時間序列的周期，
表示以年為周期，
表示以周為周期。它的傅立葉級數的形式都是：
就作者的經驗而言，對于以年為周期的序列(
)而言，
；對于以周為周期的序列(
)而言，
。這裡的參數可以形成列向量：
當
時，
當
時，
是以，時間序列的季節項就是：
而
的初始化是
。這裡的
是通過 seasonality_prior_scale 來控制的，也就是說
seasonality_prior_scale。這個值越大，表示季節的效應越明顯；這個值越小，表示季節的效應越不明顯。同時，在代碼裡面，seasonality_mode 也對應着兩種模式，分别是加法和乘法，預設是加法的形式。在開源代碼中，
函數是通過 fourier_series 來建構的。
節假日效應(holidays and events)
在現實環境中，除了周末，同樣有很多節假日，而且不同的國家有着不同的假期。在 Prophet 裡面，通過維基百科裡面對各個國家的節假日的描述，hdays.py 收集了各個國家的特殊節假日。除了節假日之外，使用者還可以根據自身的情況來設定必要的假期，例如 The Super Bowl，雙十一等。
由于每個節假日對時間序列的影響程度不一樣，例如春節，國慶節則是七天的假期，對于勞動節等假期來說則假日較短。是以，不同的節假日可以看成互相獨立的模型，并且可以為不同的節假日設定不同的前後視窗值，表示該節假日會影響前後一段時間的時間序列。用數學語言來說，對與第
個節假日來說，
表示該節假日的前後一段時間。為了表示節假日效應，我們需要一個相應的訓示函數(indicator function)，同時需要一個參數
來表示節假日的影響範圍。假設我們有
個節假日，那麼
其中
其中
并且該正态分布是受到
holidays_prior_scale 這個名額影響的。預設值是 10，當值越大時，表示節假日對模型的影響越大；當值越小時，表示節假日對模型的效果越小。使用者可以根據自己的情況自行調整。
模型拟合(Model Fitting)
按照以上的解釋，我們的時間序列已經可以通過增長項，季節項，節假日項來建構了，i.e.
下一步我們隻需要拟合函數就可以了，在 Prophet 裡面，作者使用了 pyStan 這個開源工具中的 L-BFGS 方法來進行函數的拟合。具體可以參考 forecast.py 裡面的 stan_init 函數。
Prophet 中可以設定的參數
在 Prophet 中，使用者一般可以設定以下四種參數：Capacity：在增量函數是邏輯回歸函數的時候，需要設定的容量值。
Change Points：可以通過 n_changepoints 和 changepoint_range 來進行等距的變點設定，也可以通過人工設定的方式來指定時間序列的變點。
季節性和節假日：可以根據實際的業務需求來指定相應的節假日。
光滑參數：
changepoint_prior_scale 可以用來控制趨勢的靈活度，
seasonality_prior_scale 用來控制季節項的靈活度，
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)
在進行了預測操作之後，通常都希望把時間序列的預測趨勢畫出來：
畫出預測圖：
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 則需要根據具體的時間序列來确定。
參考文獻：
![python 時間序列prophet 模型分析_時間序列預測模型Prophet[Facebook][圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)