編者按：Mail.Ru Search數(shù)據(jù)分析負責人Egor Polusmak介紹了Facebook出品的時序分析庫Prophet.

時序預(yù)測在數(shù)據(jù)分析中有廣泛應(yīng)用，例如：

在線服務(wù)明年需要的服務(wù)器數(shù)量

指定日期超市的日用品需求量

可交易的金融資產(chǎn)明天的收盤價

有相當多預(yù)測未來趨勢的不同方法，例如，ARIMA、ARCH、回歸模型、神經(jīng)網(wǎng)絡(luò)。

本文將介紹Prophet，F(xiàn)acebook在2017年2月23日開源的時序預(yù)測庫。我們也將用它預(yù)測Medium每天發(fā)表的文章數(shù)。

概覽

導(dǎo)言

Prophet預(yù)測模型

Prophet實踐

安裝

數(shù)據(jù)集

探索性可視分析

進行預(yù)測

預(yù)測質(zhì)量評估

可視化

Box-Cox變換

總結(jié)

相關(guān)資源

導(dǎo)言

據(jù)Facebook研究院的文章所言，Prophet原本是為創(chuàng)建高質(zhì)量的商業(yè)預(yù)測而研發(fā)的。Prophet嘗試處理許多商業(yè)時序數(shù)據(jù)中常見的困難：

人類行為導(dǎo)致的季節(jié)性效應(yīng)：周、月、年循環(huán)，公共假期的峰谷；

新產(chǎn)品和市場事件導(dǎo)致的趨勢變動；

離群值。

據(jù)稱在許多情形下，默認配置的Prophet就可以生成媲美經(jīng)驗豐富的分析師的預(yù)測。

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

據(jù)文章所言，Prophet的適用對象和場景很廣泛：

面向廣泛的分析師受眾，這些受眾可能在時序領(lǐng)域沒有很多經(jīng)驗；

面向廣泛的預(yù)測問題；

自動估計大量預(yù)測的表現(xiàn)，包括標出可能的問題，以供分析師進一步調(diào)查。

Prophet預(yù)測模型

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

其中：

趨勢g(t)建模非周期性變動；

季節(jié)性s(t)建模周期性變動；

假日成分h(t)提供關(guān)于假日和事件的信息。

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

趨勢

Prophet庫實現(xiàn)兩種趨勢模型。

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

其中：

C是承載量（曲線的最大值）

k是增長率（曲線的“陡峭程度”）

m是偏置參數(shù)

這一邏輯回歸等式可供建模非線性飽和增長，即數(shù)值的增長率隨增長而下降。一個典型的例子是應(yīng)用或網(wǎng)站的用戶增長。

C和k實際上不一定是常量，可能隨著時間而發(fā)生變動。Prophet支持自動和手動調(diào)整這兩個參數(shù)，既可以通過擬合提供的歷史數(shù)據(jù)自行選擇最佳的趨勢改變點，也可以讓分析師手動指定增長率和承載量變動的時間點。

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

季節(jié)性

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

年季節(jié)性通過傅里葉級數(shù)建模。

0.2版加入了新的日季節(jié)性特性，可以使用日以下尺度的時序數(shù)據(jù)并做出日以下尺度的預(yù)測。

假日和事件

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

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

誤差

最后的誤差項?t表示模型未反映的信息，通常建模為高斯噪聲。

Prophet評測

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

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

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

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

MAPE廣泛用于測量預(yù)測精確性，因為它將誤差表示為百分比，因此可以用于不同數(shù)據(jù)集上的模型評估。

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

稍微講下和Prophet作對比的算法。大多數(shù)都相當簡單，經(jīng)常用作基線：

naive是一個過度簡化的預(yù)測方法，僅僅根據(jù)上一時間點的觀測預(yù)測所有未來值。

snavie，類似naive，不過考慮了季節(jié)性因素。例如，在周季節(jié)性數(shù)據(jù)的情形下，用上周一的數(shù)據(jù)預(yù)測未來每周一的數(shù)據(jù)，用上周二的數(shù)據(jù)預(yù)測未來每周二的數(shù)據(jù)，以此類推。

mean，使用數(shù)據(jù)的平均值作為預(yù)測。

arima是自回歸集成移動平均的簡稱，參見第9課了解這一算法的細節(jié)。

ets是指數(shù)平滑的簡稱，參見第9課了解詳情。

Prophet實踐

安裝

首先安裝Prophet庫。Prophet支持Python和R語言。選擇哪種語言取決于個人偏好和項目需求。我們這里將使用Python。

Python下可以通過pip安裝：

pip install fbprophet

R也有對應(yīng)的CRAN包。

引入所需模塊并初始化環(huán)境：

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

數(shù)據(jù)集

我們將預(yù)測Medium上每天發(fā)布的文章數(shù)。

首先載入數(shù)據(jù)集（數(shù)據(jù)集下載地址見文末）：

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

接著，我們丟棄除了published（發(fā)布日期）和url（可以作為文章的唯一標識）之外的所有特征，順便移除可能存在的重復(fù)值和缺失值。

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

接下來我們需要轉(zhuǎn)換published為時間日期格式，因為pandas默認視為字符串。

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

根據(jù)時間排序dataframe，然后查看前三條記錄。

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

Medium是從2012年8月15日起開放的。然而，如你所見，至少有一些行的發(fā)布日期更早，這些極可能是不合法的。我們將清洗時序數(shù)據(jù)，限制一下時間范圍：

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)

由于需要預(yù)測發(fā)布文章數(shù)，我們將聚合、計數(shù)給定時間點的不同文章數(shù)。我們將相應(yīng)的新增列命名為posts（帖子）：

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

aggr_df.columns = ['posts']

實踐中我們感興趣的是一天發(fā)布的文章數(shù)，但是當前數(shù)據(jù)屬于日尺度以下的時序：

aggr_df.head(n=3)

為了修正這一點，我們需要根據(jù)時間“箱”聚合文章數(shù)。在時序分析中，這一過程稱為重采樣。并且，如果我們降低了數(shù)據(jù)的采樣率，那么這稱為降采樣。

幸運的是，pandas內(nèi)置這一功能：

daily_df = aggr_df.resample('D').apply(sum)

daily_df.head(n=3)

探索性可視分析

一般來說，查看數(shù)據(jù)的圖形表示可能提供幫助和指引。我們將在整個時間區(qū)間上創(chuàng)建時序圖，這可能提供季節(jié)性和明顯的異常偏離的線索。

首先我們引入并初始化Plotly庫，以創(chuàng)建美觀的交互圖：

from plotly.offline import init_notebook_mode, iplot

from plotly import graph_objs as go

# 初始化plotly

init_notebook_mode(connected=True)

我們還將定義一個用于繪圖的幫助函數(shù)：

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)

讓我們試著直接繪制數(shù)據(jù)集：

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

即使Plotly提供了縮放功能，這樣的高頻數(shù)據(jù)仍舊很難分析。除了明顯的增長和加速趨勢外，很難從這樣的圖形中推斷出什么有意義的結(jié)論。

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

我們將降采樣的dataframe保存到另一個變量中，因為之后我們將按日處理數(shù)據(jù)：

weekly_df = daily_df.resample('W').apply(sum)

plotly_df(weekly_df, title='Posts on Medium (weekly)')

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

Plotly提供的最有用的功能之一是快速深入不同時間段，這有助于更好地理解數(shù)據(jù)以及找到關(guān)于可能的趨勢、周期性、反常效應(yīng)的線索。

例如，放大連續(xù)幾年會顯示對應(yīng)基督教節(jié)日的時間點，這些對人類行為有重大影響。

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

daily_df = daily_df.loc[daily_df.index >= '2015-01-01']

daily_df.head(n=3)

基于可視化分析，我們可以看到我們的數(shù)據(jù)集呈現(xiàn)出一個顯著的不斷增長的趨勢。同時也展示了周積極性和年季節(jié)性，還有每年中的一些異常日期。

進行預(yù)測

Prophet的API和sklearn很相似。首先創(chuàng)建一個模型，接著調(diào)用fit方法，最后做出預(yù)測。fit方法的輸入是一個包含兩列的DataFrame：

ds（datestamp），類型必須是date或datetime。

y是想要預(yù)測的數(shù)值。

我們先引入庫并關(guān)閉不重要的診斷信息：

from fbprophet importProphet

import logging

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

轉(zhuǎn)換dataframe至Prophet要求的格式：

df = daily_df.reset_index()

df.columns = ['ds', 'y']

df.tail(n=3)

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

為了測量預(yù)測的質(zhì)量，我們需要將數(shù)據(jù)集分為歷史部分（數(shù)據(jù)的前部，也是最大部分）和預(yù)測部分（時間線的最后部分）。我們從數(shù)據(jù)集中移除最后一個月的數(shù)據(jù)，作為測試目標：

prediction_size = 30

train_df = df[:-prediction_size]

train_df.tail(n=3)

現(xiàn)在我們需要創(chuàng)建一個新的Prophet對象（模型），我們可以傳入模型的參數(shù)，但是這里我們將使用默認值。接著在訓(xùn)練數(shù)據(jù)集上調(diào)用fit方法訓(xùn)練模型。

m = Prophet()

m.fit(train_df);

然后我們使用輔助函數(shù)Prophet.make_future_dataframe創(chuàng)建一個dataframe，其中包括所有歷史日期，以及之前留置的未來30天。

future = m.make_future_dataframe(periods=prediction_size)

future.tail(n=3)

調(diào)用predict方法，傳入我們想要創(chuàng)建預(yù)測的日期后就可以預(yù)測新值了。如果像這里一樣同時提供歷史數(shù)據(jù)，那除了預(yù)測之外我們還能得到歷史的內(nèi)擬合。

forecast = m.predict(future)

forecast.tail(n=3)

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

Prophet庫內(nèi)置可視化工具，方便迅速評估結(jié)果。

首先，Prophet.plot方法可以繪制所有預(yù)測的數(shù)據(jù)點：

m.plot(forecast);

這張圖沒有提供多少信息。我們唯一能得出的結(jié)論是模型將許多數(shù)據(jù)點視為離群值。

在我們的情形中，Prophet.plot_components函數(shù)也許要有用得多。它讓我們可以分別觀察模型的不同成分：趨勢、年季節(jié)性、周季節(jié)性。如果給模型提供了關(guān)于假日和事件的信息，圖形也會顯示相關(guān)信息。

m.plot_components(forecast);

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

預(yù)測質(zhì)量評估

讓我們計算預(yù)測的最后30天的誤差測度，以便評估算法的質(zhì)量。為此，我們需要觀測yi和相應(yīng)的預(yù)測值?i。

看下Prophet為我們創(chuàng)建的forecast對象：

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

結(jié)果：

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

我們看到，其中沒有歷史值。我們需要從原數(shù)據(jù)集df中獲取實際值y，然后和forecast對象中的預(yù)測值比較。為此我們定義一個輔助函數(shù)：

def make_comparison_dataframe(historical, forecast):

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

應(yīng)用這一函數(shù)：

cmp_df = make_comparison_dataframe(df, forecast)

cmp_df.tail(n=3)

我們再定義一個計算MAPE和MAE的輔助函數(shù)：

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)

結(jié)果：

MAPE 22.7243579814

MAE 70.452857085

可視化

現(xiàn)在讓我們自行可視化Prophet創(chuàng)建的模型。它將包括實際值、預(yù)測值和置信區(qū)間。

首先，我們繪制較短時期內(nèi)的圖形，這樣數(shù)據(jù)點更容易分辨。接著，我們只顯示模型在預(yù)測期間內(nèi)的表現(xiàn)（最后30天）。這兩個測度看起來能帶來更清晰的圖形。

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

我們將定義一個輔助函數(shù)show_forecast并調(diào)用它：

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

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

我們還可以從圖中得出的結(jié)論是，許多實際值位于置信區(qū)間之外。Prophet也許不太適合方差不穩(wěn)定的時序，至少在默認設(shè)置下是如此。我們將通過轉(zhuǎn)換數(shù)據(jù)來嘗試修正這一點。

Box-Cox變換

目前為止，我們使用的都是Prophet的默認設(shè)置，數(shù)據(jù)也基本上是原始數(shù)據(jù)。我們這里不打算討論如何調(diào)整模型的參數(shù)，不過，即使在不動默認參數(shù)的情況下，還是有提升的空間。在這一節(jié)，我們將在原始時序上應(yīng)用Box-Cox變換，看看會有什么效果。

簡單介紹下這一變換。這一單調(diào)數(shù)據(jù)變換可以穩(wěn)定方差。我們將使用單參數(shù)Box-Cox變換：

使用上式的逆函數(shù)可以還原至原數(shù)據(jù)的尺度：

相應(yīng)的Python函數(shù)：

def inverse_boxcox(y, lambda_):

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

首先，我們準備數(shù)據(jù)，設(shè)置索引：

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

接著，我們應(yīng)用Scipy的stats.boxcox函數(shù)（Box-Cox變換）。這里它將返回兩個值，第一個值是轉(zhuǎn)換后的序列，第二個值是找到的最優(yōu)λ值（最大似然）：

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

train_df2.reset_index(inplace=True)

創(chuàng)建一個新Prophet模型，并重復(fù)之前的擬合-預(yù)測流程：

m2 = Prophet()

m2.fit(train_df2)

future2 = m2.make_future_dataframe(periods=prediction_size)

forecast2 = m2.predict(future2)

然后通過逆函數(shù)和已知的λ值反轉(zhuǎn)Box-Cox變換：

for column in ['yhat', 'yhat_lower', 'yhat_upper']:

forecast2[column] = inverse_boxcox(forecast2[column],

lambda_prophet)

復(fù)用之前創(chuàng)建對比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)

結(jié)果：

MAPE 11.5921879552

MAE 39.072031256

毫無疑問，模型的質(zhì)量提升了。

最后，讓我們繪制最新結(jié)果的可視化圖形，和之前的放一起對比下。

show_forecast(cmp_df, prediction_size, 100, 'No transformations')

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

轉(zhuǎn)換前

轉(zhuǎn)換后

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

總結(jié)

我們介紹了Prophet這一開源的時序預(yù)測庫，并進行了一些時序預(yù)測的實踐。

如我們所見，Prophet并不神奇，開箱即用的預(yù)測并不理想。它仍然需要數(shù)據(jù)科學(xué)家在必要的時候探索預(yù)測結(jié)果，調(diào)整模型參數(shù)，轉(zhuǎn)換數(shù)據(jù)。

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

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

相關(guān)資源

GitHub上的Prophet官方倉庫：https://github.com/facebook/prophet

Prophet官方文檔：https://facebookincubator.github.io/prophet/docs/quick_start.html

Sean J. Taylor和Benjamin Letham的論文Forecasting at scale解釋了作為Prophet基礎(chǔ)的算法。

Chris Moffitt寫的Prophet概覽（以預(yù)測網(wǎng)站流量為例）：http://pbpython.com/prophet-overview.html

Rob J. Hyndman和George Athanasopoulos寫的Forecasting: principles and practice——很好的關(guān)于時序預(yù)測的一本書（有在線版）

本文配套的Jupyter Notebook: git.io/fpslo

Medium數(shù)據(jù)集：https://drive.google.com/file/d/1G3YjM6mR32iPnQ6O3f6rE9BVbhiTiLyU/view

課程回顧

機器學(xué)習(xí)開放課程系列至此告一段落，所以這里列下之前的十課：

Pandas探索性分析

可視化

分類、決策樹、K近鄰

線性分類、線性回歸

Bagging、隨機森林

特征工程

PCA、聚類

Vowpal Wabbit

時序數(shù)據(jù)

梯度提升

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

連喵星人也被課程吸引（來源：Yulia Kameneva）

配套Kaggle競賽：

基于網(wǎng)頁會話檢測惡意用戶：https://www.kaggle.com/c/catch-me-if-you-can-intruder-detection-through-webpage-session-tracking2

預(yù)測Medium文章推薦數(shù)：https://www.kaggle.com/c/how-good-is-your-medium-article/