時(shí)序預測的三種方式：統計學(xué)模型、機器學(xué)習、循環(huán)神經(jīng)網(wǎng)絡(luò )

作者 | luanhz來(lái)源 | 小數志

導讀

時(shí)序預測是一類(lèi)經(jīng)典的問(wèn)題，在學(xué)術(shù)界和工業(yè)界都有著(zhù)廣泛的研究和應用。甚至說(shuō)，世間萬(wàn)物加上時(shí)間維度后都可抽象為時(shí)間序列問(wèn)題，例如股****價(jià)格、天氣變化等等。關(guān)于時(shí)序預測問(wèn)題的相關(guān)理論也極為廣泛，除了經(jīng)典的各種統計學(xué)模型外，當下火熱的機器學(xué)習以及深度學(xué)習中的循環(huán)神經(jīng)網(wǎng)絡(luò )也都可以用于時(shí)序預測問(wèn)題的建模。今天，本文就來(lái)介紹三種方式的簡(jiǎn)單應用，并在一個(gè)真實(shí)的時(shí)序數據集上加以驗證。

時(shí)間序列預測，其主要任務(wù)是基于某一指標的歷史數據來(lái)預測其在未來(lái)的取值，例如上圖中的曲線(xiàn)記錄了1949年至1960年共12年144個(gè)月份的每月航班乘客數（具體單位未經(jīng)考證），那么時(shí)序預測要解決的問(wèn)題就是：給定前9年的歷史數據，例如1949-1957，那么能否預測出1958-1960兩年間的乘客數量的問(wèn)題。

為了解決這一問(wèn)題，大概當前主流的解決方式有4種：

• 統計學(xué)模型，較為經(jīng)典的AR系列，包括AR、MA、ARMA以及ARIMA等，另外Facebook（準確的講，現在應該叫Meta了）推出的Prophet模型，其實(shí)本質(zhì)上也是一種統計學(xué)模型，只不過(guò)是傳統的趨勢、周期性成分的基礎上，進(jìn)一步細化考慮了節假日、時(shí)序拐點(diǎn)等因素的影響，以期帶來(lái)更為精準的時(shí)序規律刻畫(huà)；

• 機器學(xué)習模型，在有監督機器學(xué)習中，回歸問(wèn)題主要解決的是基于一系列Feature來(lái)預測某一Label的可能取值的問(wèn)題，那么當以歷史數據作為Feature時(shí)其實(shí)自然也就可以將時(shí)序預測問(wèn)題抽象為回歸問(wèn)題，從這一角度講，所有回歸模型都可用于解決時(shí)序預測。關(guān)于用機器學(xué)習抽象時(shí)序預測，推薦查看這篇論文《Machine Learning Strategies for Time Series Forecasting》；

• 深度學(xué)習模型，深度學(xué)習主流的應用場(chǎng)景當屬CV和NLP兩大領(lǐng)域，其中后者就是專(zhuān)門(mén)用于解決序列問(wèn)題建模的問(wèn)題，而時(shí)間序列當然屬于序列數據的一種特殊形式，所以自然可以運用循環(huán)神經(jīng)網(wǎng)絡(luò )來(lái)建模時(shí)序預測；

• 隱馬爾科夫模型，馬爾科夫模型是用于刻畫(huà)相鄰狀態(tài)轉換間的經(jīng)典抽象，而隱馬爾科夫模型則在其基礎上進(jìn)一步增加了隱藏狀態(tài)，來(lái)以此豐富模型的表達能力。但其一大假設條件是未來(lái)狀態(tài)僅與當前狀態(tài)有關(guān)，而不利于利用多個(gè)歷史狀態(tài)來(lái)共同參與預測，較為常用的可能就是天氣預報的例子了。

本文主要考慮前三種時(shí)序預測建模方法，并分別選?。?）Prophet模型，2）RandomForest回歸模型，3）LSTM三種方案加以測試。

首先在這個(gè)航班乘客真實(shí)數據集上進(jìn)行測試，依次對比三個(gè)所選模型的預測精度。該數據集共有12年間每個(gè)月的乘客數量，以1958年1月作為切分界面劃分訓練集和測試集，即前9年的數據作為訓練集，后3年的數據作為測試集驗證模型效果。數據集切分后的示意圖如下：

df = pd.read_csv("AirPassengers.csv", parse_dates=["date"]).rename(columns={"date":"ds", "value":"y"})X_train = df[df.ds<"19580101"]X_test = df[df.ds>="19580101"]
plt.plot(X_train['ds'], X_train['y'])plt.plot(X_test['ds'], X_test['y'])

1.Prophet模型預測。Prophet是一個(gè)高度封裝好的時(shí)序預測模型，接受一個(gè)DataFrame作為訓練集（要求有ds和y兩個(gè)字段列），在預測時(shí)也接受一個(gè)DataFrame，但此時(shí)只需有ds列即可，關(guān)于模型的詳細介紹可參考其官方文檔：https://facebook.github.io/prophet/。模型訓練及預測部分核心代碼如下：

from prophet import Prophetpro = Prophet()pro.fit(X_train)pred = pro.predict(X_test)
pro.plot(pred)

訓練后的結果示意圖如下：

當然，這是通過(guò)Prophet內置的可視化函數給出的結果，也可通過(guò)手動(dòng)繪制測試集真實(shí)標簽與預測結果間的對比：

易見(jiàn)，雖然序列的整體****上具有良好的擬合結果，但在具體取值上其實(shí)差距還是比較大的。

2.機器學(xué)習模型，這里選用常常用作各種baseline的RandomForest模型。在使用機器學(xué)習實(shí)現時(shí)序預測時(shí)，通常需要通過(guò)滑動(dòng)窗口的方式來(lái)提取特征和標簽，而后在實(shí)現預測時(shí)實(shí)際上也需滑動(dòng)的截取測試集特征實(shí)現單步預測，參考論文《Machine Learning Strategies for Time Series Forecasting》中的做法，該問(wèn)題可大致描述如下：

據此，設置特征提取窗口長(cháng)度為12，構建訓練集和測試集的方式如下：

data = df.copy()n = 12for i in range(1, n+1):    data['ypre_'+str(i)] = data['y'].shift(i)data = data[['ds']+['ypre_'+str(i) for i in range(n, 0, -1)]+['y']]
# 提取訓練集和測試集X_train = data[data['ds']<"19580101"].dropna()[['ypre_'+str(i) for i in range(n, 0, -1)]]y_train = data[data['ds']<"19580101"].dropna()[['y']]X_test = data[data['ds']>="19580101"].dropna()[['ypre_'+str(i) for i in range(n, 0, -1)]]y_test = data[data['ds']>="19580101"].dropna()[['y']]
# 模型訓練和預測rf = RandomForestRegressor(n_estimators=10, max_depth=5)rf.fit(X_train, y_train)y_pred = rf.predict(X_test)
# 結果對比繪圖y_test.assign(yhat=y_pred).plot()

可見(jiàn)，預測效果也較為一般，尤其是對于最后兩年的預測結果，與真實(shí)值差距還是比較大的。用機器學(xué)習模型的思維很容易解釋這一現象：隨機森林模型實(shí)際上是在根據訓練數據集來(lái)學(xué)習曲線(xiàn)之間的規律，由于該時(shí)序整體呈現隨時(shí)間增長(cháng)的趨勢，所以歷史數據中的最高點(diǎn)也不足以cover住未來(lái)的較大值，因而在測試集中超過(guò)歷史數據的所有標簽其實(shí)都是無(wú)法擬合的。

3.深度學(xué)習中的循環(huán)神經(jīng)網(wǎng)絡(luò )，其實(shí)深度學(xué)習一般要求數據集較大時(shí)才能發(fā)揮其優(yōu)勢，而這里的數據集顯然是非常小的，所以?xún)H設計一個(gè)最為簡(jiǎn)單的模型：1層LSTM+1層Linear。模型搭建如下：

class Model(nn.Module):    def __init__(self):        super().__init__()        self.rnn = nn.LSTM(input_size=1, hidden_size=10, batch_first=True)        self.linear = nn.Linear(10, 1)
    def forward(self, x):        x, _ = self.rnn(x)        x = x[:, -1, :]        x = self.linear(x)        return x

數據集構建思路整體同前述的機器學(xué)習部分，而后，按照進(jìn)行模型訓練煉丹，部分結果如下：

# 數據集轉化為3DX_train_3d = torch.Tensor(X_train.values).reshape(*X_train.shape, 1)y_train_2d = torch.Tensor(y_train.values).reshape(*y_train.shape, 1)X_test_3d = torch.Tensor(X_test.values).reshape(*X_test.shape, 1)y_test_2d = torch.Tensor(y_test.values).reshape(*y_test.shape, 1)
# 模型、優(yōu)化器、評估準則model = Model()creterion = nn.MSELoss()optimizer = optim.Adam(model.parameters())
# 訓練過(guò)程for i in range(1000):    out = model(X_train_3d)    loss = creterion(out, y_train_2d)    optimizer.zero_grad()    loss.backward()    optimizer.step()
    if (i+1)%100 == 0:        y_pred = model(X_test_3d)        loss_test = creterion(y_pred, y_test_2d)        print(i, loss.item(), loss_test.item())
# 訓練結果99 65492.08984375 188633.796875199 64814.4375 187436.4375299 64462.09765625 186815.5399 64142.70703125 186251.125499 63835.5 185707.46875599 63535.15234375 185175.1875699 63239.39453125 184650.46875799 62947.08203125 184131.21875899 62657.484375 183616.203125999 62370.171875 183104.671875

通過(guò)上述1000個(gè)epoch，大體可以推斷該模型不會(huì )很好的擬合了，所以果斷放棄吧！

當然必須指出的是，上述測試效果只能說(shuō)明3種方案在該數據集上的表現，而不能代表這一類(lèi)模型在用于時(shí)序預測問(wèn)題時(shí)的性能。實(shí)際上，時(shí)序預測問(wèn)題本身就是一個(gè)需要具體問(wèn)題具體分析的場(chǎng)景，沒(méi)有放之四海而皆準的好模型，就像“No Free Lunch”一樣！

本文僅是作為時(shí)序預測系列推文的一個(gè)牛刀小試，后續將不定期更新其他相關(guān)心得和總結。