機器學習開放課程（終）：基於Facebook Prophet預測未來

資料探勘 · 發表 2018-11-19 11:40:12

摘要：作者：Egor Polusmak 編譯：weakish 時序預測在資料分析中有廣泛應用，例如：線上服務明年需要的伺服器數量指定日期超市的日用品需求量可交易的金融資產明天的收盤價有相當多預測未來趨勢的不同方法，例如，ARIMA、ARC...

作者：Egor Polusmak

編譯：weakish

時序預測在資料分析中有廣泛應用，例如：

線上服務明年需要的伺服器數量
指定日期超市的日用品需求量
可交易的金融資產明天的收盤價

有相當多預測未來趨勢的不同方法，例如，ARIMA、ARCH 、迴歸模型、神經網路。

本文將介紹Prophet，Facebook在2017年2月23日開源的時序預測庫。我們也將用它預測Medium每天發表的文章數。

概覽

導言

Prophet預測模型

Prophet實踐

安裝
資料集
探索性可視分析
進行預測
預測質量評估
視覺化

Box-Cox變換

總結

相關資源

導言

據Facebook研究院的文章所言，Prophet原本是為建立高質量的商業預測而研發的。Prophet嘗試處理許多商業時序資料中常見的困難：

人類行為導致的季節性效應：周、月、年迴圈，公共假期的峰谷；
新產品和市場事件導致的趨勢變動；
離群值。

據稱在許多情形下，預設配置的Prophet就可以生成媲美經驗豐富的分析師的預測。

此外，Prophet有一些直觀而易於解釋的定製功能，以供逐漸改善預測模型。特別重要的是，即使不是時序分析的專家，也可以理解這些引數。

據文章所言，Prophet的適用物件和場景很廣泛：

面向廣泛的分析師受眾，這些受眾可能在時序領域沒有很多經驗；
面向廣泛的預測問題；
自動估計大量預測的表現，包括標出可能的問題，以供分析師進一步調查。

Prophet預測模型

下面讓我們仔細看看Prophet是如何工作的。本質上，這個庫使用的是加性迴歸模型：y(t) = g(t) + s(t) + h(t) + ϵt

其中：

趨勢g(t)建模非週期性變動；
季節性s(t)建模週期性變動；
假日成分h(t)提供關於假日和事件的資訊。

下面我們將討論這些模型成分的一些重要性質。

趨勢

Prophet庫實現兩種趨勢模型。

第一種是 非線性飽和增長 。它可以用邏輯增長模型表示：

其中：

C是承載量（曲線的最大值）
k是增長率（曲線的“陡峭程度”）
m是偏置引數

這一邏輯迴歸等式可供建模非線性飽和增長，即數值的增長率隨增長而下降。一個典型的例子是應用或網站的使用者增長。

C和k實際上不一定是常量，可能隨著時間而發生變動。Prophet支援自動和手動調整這兩個引數，既可以通過擬合提供的歷史資料自行選擇最佳的趨勢改變點，也可以讓分析師手動指定增長率和承載量變動的時間點。

第二種趨勢模型是增長率恆定的簡單 分段線性模型 ，最適合不存在飽和的線性增長。

季節性

周季節性通過虛擬變數建模： monday 、 tuesday 、 wednesday 、 thursday 、 friday 、 saturday （週一到週六）。這些變數的值為0還是為1取決於是星期幾。 sunday （週日）變數沒有加入，因為上述六個變數都取0即可表示週日。相反，如果再引入週日變數，那麼每個變數都可以通過另外六個變數的線性組合表示，這種變數之間的相關性會對模型產生不利影響。

年季節性通過傅立葉級數建模。

ofollow,noindex">0.2版加入了新的 日季節性 特性，可以使用日以下尺度的時序資料並做出日以下尺度的預測。

假日和事件

h(t)表示每年的可預測反常日期，例如黑色星期五。

分析師需要提供定製的事件列表以利用這一特性。

誤差

最後的誤差項ϵt表示模型未反映的資訊，通常建模為高斯噪聲。

Prophet評測

Facebook的論文（見相關資源）對比了Prophet和其他幾種時序預測方法。根據他們的研究，相比其他模型，Prophet顯著降低了預測誤差（使用平均絕對百分誤差測量預測精確度）。

為了便於理解上面的評測結果，下面簡單溫習下平均絕對百分誤差（MAPE）的概念。

將 實際（歷史）值 記為yi，模型給出的 預測值 記為ŷi。那麼 預測誤差 ei = yi - ŷi， 相對預測誤差 pi = ei/yi.

由此，我們定義MAPE = mean(|pi|)

MAPE廣泛用於測量預測精確性，因為它將誤差表示為百分比，因此可以用於不同資料集上的模型評估。

此外，評估預測演算法時，為了瞭解誤差的絕對值，可以計算平均絕對誤差（MAE）：MAE = mean(|ei|)

稍微講下和Prophet作對比的演算法。大多數都相當簡單，經常用作基線：

naive 是一個過度簡化的預測方法，僅僅根據上一時間點的觀測預測所有未來值。
snavie ，類似 naive ，不過考慮了季節性因素。例如，在周季節性資料的情形下，用上週一的資料預測未來每週一的資料，用上週二的資料預測未來每週二的資料，以此類推。
mean ，使用資料的平均值作為預測。
arima 是 自迴歸整合移動平均 的簡稱，參見第9課瞭解這一演算法的細節。
ets 是 指數平滑 的簡稱，參見第9課瞭解詳情。

Prophet實踐

安裝

首先安裝Prophet庫。Prophet支援Python和R語言。選擇哪種語言取決於個人偏好和專案需求。我們這裡將使用Python。

Python下可以通過pip安裝：

pip install fbprophet

R也有對應的CRAN包。

引入所需模組並初始化環境：

import warnings
warnings.filterwarnings('ignore')

import numpy as np
import pandas as pd
from scipy import stats
import statsmodels.api as sm
import matplotlib.pyplot as plt

%matplotlib inline

資料集

我們將預測Medium上每天釋出的文章數。

首先載入資料集（資料集下載地址見文末）：

df = pd.read_csv('../../data/medium_posts.csv', sep='\t')

接著，我們丟棄除了 published （釋出日期）和 url （可以作為文章的唯一標識）之外的所有特徵，順便移除可能存在的重複值和缺失值。

df = df[['published', 'url']].dropna().drop_duplicates()

接下來我們需要轉換 published 為時間日期格式，因為pandas預設視為字串。

df['published'] = pd.to_datetime(df['published'])

根據時間排序dataframe，然後檢視前三條記錄。

df.sort_values(by=['published']).head(n=3)

Medium是從2012年8月15日起開放的。然而，如你所見，至少有一些行的釋出日期更早，這些極可能是不合法的。我們將清洗時序資料，限制一下時間範圍：

df = df[(df['published'] > '2012-08-15') & 
(df['published'] < '2017-06-26')].
sort_values(by=['published'])
df.head(n=3)

最後3條：

df.tail(n=3)

由於需要預測釋出文章數，我們將聚合、計數給定時間點的不同文章數。我們將相應的新增列命名為 posts （帖子）：

aggr_df = df.groupby('published')[['url']].count()
aggr_df.columns = ['posts']

實踐中我們感興趣的是一天釋出的文章數，但是當前資料屬於 日尺度以下的時序 ：

aggr_df.head(n=3)

為了修正這一點，我們需要根據時間“箱”聚合文章數。在時序分析中，這一過程稱為 重取樣 。並且，如果我們降低了資料的取樣率，那麼這稱為 降取樣 。

幸運的是，pandas內建這一功能：

daily_df = aggr_df.resample('D').apply(sum)
daily_df.head(n=3)

探索性可視分析

一般來說，檢視資料的圖形表示可能提供幫助和指引。我們將在整個時間區間上建立時序圖，這可能提供季節性和明顯的異常偏離的線索。

首先我們引入並初始化Plotly庫，以建立美觀的互動圖：

from plotly.offline import init_notebook_mode, iplot
from plotly import graph_objs as go

# 初始化plotly
init_notebook_mode(connected=True)

我們還將定義一個用於繪圖的幫助函式：

def plotly_df(df, title=''):
"""視覺化dataframe所有列為線圖。"""
common_kw = dict(x=df.index, mode='lines')
data = [go.Scatter(y=df[c], name=c, **common_kw) for c in df.columns]
layout = dict(title=title)
fig = dict(data=data, layout=layout)
iplot(fig, show_link=False)

讓我們試著直接繪製資料集：

plotly_df(daily_df, title='Posts on Medium (daily)')

即使Plotly提供了縮放功能，這樣的高頻資料仍舊很難分析。除了明顯的增長和加速趨勢外，很難從這樣的圖形中推斷出什麼有意義的結論。

為了減少噪聲，我們將按周重取樣文章數。順便提下，分箱之外，其他降噪的技術還包括移動平均平滑和指數平滑等。

我們將降取樣的dataframe儲存到另一個變數中，因為之後我們將按日處理資料：

weekly_df = daily_df.resample('W').apply(sum)
plotly_df(weekly_df, title='Posts on Medium (weekly)')

從幫助分析師預測的角度來說，這張圖的效果更好。

Plotly提供的最有用的功能之一是快速深入不同時間段，這有助於更好地理解資料以及找到關於可能的趨勢、週期性、反常效應的線索。

例如，放大連續幾年會顯示對應基督教節日的時間點，這些對人類行為有重大影響。

根據上圖，我們將省略2015年之前的觀測。頭幾年每天釋出的文章數很低，可能給預測增加噪聲，因為模型可能不得不擬合這些異常的歷史資料而不能更好地利用最近幾年更相關、更具指示性的資料。

daily_df = daily_df.loc[daily_df.index >= '2015-01-01']
daily_df.head(n=3)

基於視覺化分析，我們可以看到我們的資料集呈現出一個顯著的不斷增長的趨勢。同時也展示了周積極性和年季節性，還有每年中的一些異常日期。

進行預測

Prophet的API和 sklearn 很相似。首先建立一個模型，接著呼叫 fit 方法，最後做出預測。 fit 方法的輸入是一個包含兩列的 DataFrame ：

ds （datestamp），型別必須是 date 或 datetime 。
y 是想要預測的數值。

我們先引入庫並關閉不重要的診斷資訊：

from fbprophet import Prophet
import logging

logging.getLogger().setLevel(logging.ERROR)

轉換dataframe至Prophet要求的格式：

df = daily_df.reset_index()
df.columns = ['ds', 'y']
df.tail(n=3)

Prophet的作者建議至少根據幾個月的資料進行預測，如果有超過一年的歷史資料就最理想了。很幸運，我們這裡有好幾年的資料可供模型擬合。

為了測量預測的質量，我們需要將資料集分為 歷史部分 （資料的前部，也是最大部分）和 預測部分 （時間線的最後部分）。我們從資料集中移除最後一個月的資料，作為測試目標：

prediction_size = 30
train_df = df[:-prediction_size]
train_df.tail(n=3)

現在我們需要建立一個新的Prophet物件（模型），我們可以傳入模型的引數，但是這裡我們將使用預設值。接著在訓練資料集上呼叫 fit 方法訓練模型。

m = Prophet()
m.fit(train_df);

然後我們使用輔助函式 Prophet.make_future_dataframe 建立一個dataframe，其中包括所有歷史日期，以及之前留置的未來30天。

future = m.make_future_dataframe(periods=prediction_size)
future.tail(n=3)

呼叫 predict 方法，傳入我們想要建立預測的日期後就可以預測新值了。如果像這裡一樣同時提供歷史資料，那除了預測之外我們還能得到歷史的內擬合。

forecast = m.predict(future)
forecast.tail(n=3)

在結果dataframe中，我們可以看到很多描繪預測的列，包括趨勢成分、季節性成分，還有它們的置信區間。預測本身儲存於 yhat 列。

Prophet庫內建視覺化工具，方便迅速評估結果。

首先， Prophet.plot 方法可以繪製所有預測的資料點：

m.plot(forecast);

這張圖沒有提供多少資訊。我們唯一能得出的結論是模型將許多資料點視為離群值。

在我們的情形中， Prophet.plot_components 函式也許要有用得多。它讓我們可以分別觀察模型的不同成分：趨勢、年季節性、周季節性。如果給模型提供了關於假日和事件的資訊，圖形也會顯示相關資訊。

m.plot_components(forecast);

從趨勢圖中，我們可以看到，Prophet很好地擬合了2016年末後的加速增長趨勢。從周季節性圖中，我們可以得出結論，和工作日相比，週六、週日的新文章較少。而年季節性圖清楚地顯示了聖誕節那一天的迅猛下跌。

預測質量評估

讓我們計算預測的最後30天的誤差測度，以便評估演算法的質量。為此，我們需要觀測yi和相應的預測值ŷi。

看下Prophet為我們建立的 forecast 物件：

print(', '.join(forecast.columns))

結果：

ds, trend, trend_lower, trend_upper, yhat_lower, yhat_upper, seasonal, seasonal_lower, seasonal_upper, seasonalities, seasonalities_lower, seasonalities_upper, weekly, weekly_lower, weekly_upper, yearly, yearly_lower, yearly_upper, yhat

我們看到，其中沒有歷史值。我們需要從原資料集 df 中獲取實際值 y ，然後和 forecast 物件中的預測值比較。為此我們定義一個輔助函式：

def make_comparison_dataframe(historical, forecast):
return forecast.set_index('ds')[['yhat', 'yhat_lower', 'yhat_upper']].join(historical.set_index('ds'))

應用這一函式：

cmp_df = make_comparison_dataframe(df, forecast)
cmp_df.tail(n=3)

我們再定義一個計算MAPE和MAE的輔助函式：

def calculate_forecast_errors(df, prediction_size):
df = df.copy()

df['e'] = df['y'] - df['yhat']
df['p'] = 100 * df['e'] / df['y']

predicted_part = df[-prediction_size:]

error_mean = lambda error_name: np.mean(np.abs(predicted_part[error_name]))

return {'MAPE': error_mean('p'), 'MAE': error_mean('e')}

然後用它計算MAPE和MAE：

for err_name, err_value in 
 calculate_forecast_errors(cmp_df, prediction_size).items():
print(err_name, err_value)

結果：

MAPE 22.7243579814
MAE 70.452857085

視覺化

現在讓我們自行視覺化Prophet建立的模型。它將包括實際值、預測值和置信區間。

首先，我們繪製較短時期內的圖形，這樣資料點更容易分辨。接著，我們只顯示模型在預測期間內的表現（最後30天）。這兩個測度看起來能帶來更清晰的圖形。

最後，我們將使用Plotly讓圖形可以互動，這對探索很重要。

我們將定義一個輔助函式 show_forecast 並呼叫它：

def show_forecast(cmp_df, num_predictions, num_values, title):
def create_go(name, column, num, **kwargs):
points = cmp_df.tail(num)
args = dict(name=name, x=points.index, y=points[column], mode='lines')
args.update(kwargs)
return go.Scatter(**args)

lower_bound = create_go('Lower Bound', 'yhat_lower', num_predictions,
line=dict(width=0),
marker=dict(color="444"))
upper_bound = create_go('Upper Bound', 'yhat_upper', num_predictions,
line=dict(width=0),
marker=dict(color="444"),
fillcolor='rgba(68, 68, 68, 0.3)', 
fill='tonexty')
forecast = create_go('Forecast', 'yhat', num_predictions,
line=dict(color='rgb(31, 119, 180)'))
actual = create_go('Actual', 'y', num_values,
marker=dict(color="red"))

data = [lower_bound, upper_bound, forecast, actual]

layout = go.Layout(yaxis=dict(title='Posts'), title=title, showlegend = False)
fig = go.Figure(data=data, layout=layout)
iplot(fig, show_link=False)

show_forecast(cmp_df, prediction_size, 100, 'New posts on Medium')

初看起來模型預測的均值還挺合理的。從圖上看，之所以之前算出來的MAPE值比較高，可能是因為模型沒能捕捉周季節性增長高峰的放大效應。

我們還可以從圖中得出的結論是，許多實際值位於置信區間之外。Prophet也許不太適合方差不穩定的時序，至少在預設設定下是如此。我們將通過轉換資料來嘗試修正這一點。

Box-Cox變換

目前為止，我們使用的都是Prophet的預設設定，資料也基本上是原始資料。我們這裡不打算討論如何調整模型的引數，不過，即使在不動預設引數的情況下，還是有提升的空間。在這一節，我們將在原始時序上應用Box-Cox變換，看看會有什麼效果。

簡單介紹下這一變換。這一單調資料變換可以穩定方差。我們將使用單引數Box-Cox變換：

使用上式的逆函式可以還原至原資料的尺度：

相應的Python函式：

def inverse_boxcox(y, lambda_):
return np.exp(y) if lambda_ == 0 else np.exp(np.log(lambda_ * y + 1) / lambda_)

首先，我們準備資料，設定索引：

train_df2 = train_df.copy().set_index('ds')

接著，我們應用 Scipy 的 stats.boxcox 函式（Box-Cox變換）。這裡它將返回兩個值，第一個值是轉換後的序列，第二個值是找到的最優λ值（最大似然）：

train_df2['y'], lambda_prophet = stats.boxcox(train_df2['y'])
train_df2.reset_index(inplace=True)

建立一個新Prophet模型，並重復之前的擬合-預測流程：

m2 = Prophet()
m2.fit(train_df2)
future2 = m2.make_future_dataframe(periods=prediction_size)
forecast2 = m2.predict(future2)

然後通過逆函式和已知的λ值反轉Box-Cox變換：

for column in ['yhat', 'yhat_lower', 'yhat_upper']:
forecast2[column] = inverse_boxcox(forecast2[column],
lambda_prophet)

複用之前建立對比dataframe和計算誤差的程式碼：

cmp_df2 = make_comparison_dataframe(df, forecast2)
for err_name, err_value in calculate_forecast_errors(cmp_df2, prediction_size).items():
print(err_name, err_value)

結果：

MAPE 11.5921879552
MAE 39.072031256

毫無疑問，模型的質量提升了。

最後，讓我們繪製最新結果的視覺化圖形，和之前的放一起對比下。

show_forecast(cmp_df, prediction_size, 100, 'No transformations')
show_forecast(cmp_df2, prediction_size, 100, 'Box–Cox transformation')

很明顯，第二張圖中的預測值更加接近真實值。

總結

我們介紹了Prophet這一開源的時序預測庫，並進行了一些時序預測的實踐。

如我們所見，Prophet並不神奇，開箱即用的預測並不理想。它仍然需要資料科學家在必要的時候探索預測結果，調整模型引數，轉換資料。

不過，這一個使用者友好、易於定製的庫。在某些情形下，單單將分析師事先知道的異常日期納入考慮這一功能就會帶來很大的不同。

總的來說，Prophet庫值得收入你的分析工具箱。

課程回顧

機器學習開放課程系列至此告一段落，所以這裡列下之前的十課：

配套視訊（目前更新到第6課）： https:// youtu.be/QKTuw4PNOsU

配套Kaggle競賽：

基於網頁會話檢測惡意使用者： https://www. kaggle.com/c/catch-me-i f-you-can-intruder-detection-through-webpage-session-tracking2
預測Medium文章推薦數： https://www. kaggle.com/c/how-good-i s-your-medium-article/

課程最新進展： https:// mlcourse.ai/news

機器學習開放課程（終）：基於Facebook Prophet預測未來

概覽

導言

Prophet預測模型

Prophet實踐

Box-Cox變換

總結

相關資源

課程回顧

您可能也會喜歡…