本章節分兩部分，一是ARMA模型預測，另一個是LSTM模型預測。過程中會結合數學原理探求本質，以便理解兩個模型擬合效果的比較。

一. ARMA模型預測

我們對模型arimatest=ARIMA(2,0,4)×(1,1,1,7)進行擬合后，做向前7步預測：

> # 預測 > fc<-forecast(arimatest,h=7,level=c(99.5)) > fcPoint Forecast Lo 99.5 Hi 99.5 162.1429 36285.83 -37121.33 109693.0 162.2857 49331.44 -38293.45 136956.3 162.4286 38016.79 -52834.00 128867.6 162.5714 37847.55 -54204.97 129900.1 162.7143 113477.48 20428.50 206526.5 162.8571 110839.43 17751.50 203927.4 163.0000 38841.40 -54396.23 132079.0

上述輸出結果左邊第一列Point值是將原序列按周排序，第一天是2015-9-30，位于2015年第39周的周三，樣本內最后一天是2018-2-5，位于第162周的周一。因此樣本外的第一天位于第162周的周二，即2018-2-6，所以根據frequency=7，系統自動將每一周切分成7份，看到的就是小數了。

注意，它和predict方法效果一樣：

> # 預測 > pred<-predict(arimatest,7) > pred $`pred` Time Series: Start = c(162, 2) End = c(163, 1) Frequency = 7 [1] 36285.83 49331.44 38016.79 37847.55 113477.48 110839.43 38841.40$se Time Series: Start = c(162, 2) End = c(163, 1) Frequency = 7 [1] 26151.15 31216.19 32365.40 32793.52 33148.51 33162.38 33215.72> accuracy(pred$pred,test_income[1:7])ME RMSE MAE MPE MAPE Test set -20273.42 32945.89 22699.64 -46.34546 51.47128

我們看看樣本外真實數據的實際對比：

模型擬合的值偏高，這是因為這周正好趕在春節前一周，可能對收入產生一定影響（負向擾動），均方根誤差為32945元，說明誤差很大，模型擬合的并不夠好，很多波動因素沒有考慮進去。我們再用圖形直觀地看下：

# 比較 ts.plot(test_income[1:7], pred$pred, lty = c(1,3)) abline(reg=lm(test_income[1:7]~pred$pred))

從上圖看出，向前4步的預測效果都比較好，第5/6步預測（周末2天）預測偏差較大，實際值明顯偏低，第7步預測也比較好。

我們再看下向前30步預測：

很明顯，超過7天的預測效果較差，這是由兩方面原因來解釋：其一是因為我們的ARIMA模型是一個季節模型，預測也會呈現明顯的7天周期性，同時根據ARMA預測模型的理論，向前預測的估計值就是數學期望（條件無偏最小方差估計），且預測的方差（置信區間）隨預測步長增加也會越來越大，所以只適合短期預測，不過可以通過動態（修正）預測來進行改善；其二是因為模型只提取了線性部分，對擾動沒有考慮，之前文章已驗證它并不是正態隨機波動，而是有明顯的ARCH效應，所以這樣的擾動不能被忽略，比如上圖中第二個7步預測正好是春節期間，所以收入有一個明顯增長（正向擾動）。

為了得到更強有力的預測結論，利用交叉驗證的方法，將原序列860天數據拆分成不同的訓練集和驗證集，其中訓練集包含90天數據，驗證集是緊隨其后的7天數據，然后不斷往前平移30天，得到25組數據集。R代碼如下：

# 樣本重新取樣，交叉驗證 library(rsample) periods_train <- 90 periods_test <- 7 skip_span <- 30ro_resamples <- rolling_origin(data,initial = periods_train,assess = periods_test,cumulative = FALSE,skip = skip_span)win.graph(width = 10,height = 10,pointsize = 8) par(mfcol=c(5,5)) for(i in 1:dim(ro_resamples)[1]){ro1_train<-analysis(ro_resamples$splits[[i]])startW <- as.numeric(strftime(head(ro1_train$賬期, 1), format = "%W"))startD <- as.numeric(strftime(head(ro1_train$賬期, 1), format =" %w"))ro1_train_data = ts(ro1_train$本日收入, frequency = 7,start = c(startW, startD))ro1_test<-assessment(ro_resamples$splits[[i]])startW_test<-(end(ro1_train_data)[1]*7+end(ro1_train_data)[2]+1)%/%7startD_test<-(end(ro1_train_data)[1]*7+end(ro1_train_data)[2]+1)%%7ro1_test_data = ts(ro1_test$本日收入, frequency = 7,start = c(startW_test, startD_test))# autoarima = auto.arima(ro1_train,trace=F)# print(autoarima$arma)# print(autoarima$aic)arimafit<-arima(ro1_train_data,order=c(1,0,2),seasonal=list(order=c(1, 1, 1),period=7))pred<-predict(arimafit,7)par(new=F)ts.plot(ro1_test_data, pred$pred, lty = c(1,3), col=1:8) }

25組數據用同一ARIMA季節模型的7天預測效果如下：

設置平移10天，得到70組，7天預測效果如下：

設置平移45天，得到17組，7天預測效果如下：

紅色曲線表示向前7步預測，黑色曲線為驗證集7天數據，可看出不論設置的平移天數是大是小，大部分情形下模型都能捕捉預測到其趨勢走向，有的7天預測效果非常好，有的則稍差一些。總的來說，通過這些分析，我們可得出以下結論：

門店日收入是一個以 7 天為周期的季節性生意，是非平穩的，經季節差分調整的序列平穩，它包含一個線性平穩模型，以及一個不穩定的波動（條件異方差）模型，且波動不是一個線性模式，這樣的波動或擾動從業務上解釋就是通常發生在節假日，當然也包括其它情形，比如門店人為因素、天氣因素、經濟環境因素等等，而且它還很可能是非對稱或時間不可逆的，即增長和下降幅度并不相同。

二. 非線性模型

以上，讓我們開始重新思考，考察原序列的分布是否滿足獨立同分布（IID，iid或i.i.d.，Independently and Identically Distributed），IID是機器學習領域的重要假設，它要求訓練數據和測試數據是滿足相同分布的，這是使得訓練數據獲得的模型能夠在測試集獲得好的效果的一個基本保障。

主要通過散點圖和BDS檢驗兩種方法。散點圖是看序列和過去延遲階序列的二維聯合分布是否為正態（橢圓形云狀），BDS檢驗是將序列和過去延遲階序列作為一個點投影到向量空間，延遲一階就是二維，延遲兩階就是三維，如果獨立同分布，則滿足高維空間的相關性求和應該是低維空間相關性求和的冪乘關系，這里定義的相關性求和實質是兩個點的距離小于給定參數??的點對數目的期望極限。

我們繪制日收入序列的自相關12階延遲的散點圖，選12階是因為運行ar函數系統根據AIC給出的自回歸階數：

# ar函數查看滯后階 ar(income_data,method = 'mle')$order# IID獨立同分布檢驗 win.graph(width = 10,height = 10,pointsize = 8) lag.plot(income_data,lags = 12)

如果是獨立同分布，尤其針對二元正態分布，其散點圖形狀應該呈現一個橢圓，此處很明顯數據是非正態的，則相關過程也必然是非線性的。

再看BDS檢驗，直接調用R中的bds.test函數：

> iid_result<-bds.test(income_data) > iid_resultBDS Test data: income_data Embedding dimension = 2 3 Epsilon for close points = 24472.58 48945.16 73417.74 97890.32 Standard Normal = [ 24472.5803 ] [ 48945.1605 ] [ 73417.7408 ] [ 97890.3211 ] [ 2 ] 29.9547 24.5187 20.8156 19.5068 [ 3 ] 32.6452 23.3057 18.5383 17.3194p-value = [ 24472.5803 ] [ 48945.1605 ] [ 73417.7408 ] [ 97890.3211 ] [ 2 ] 0 0 0 0 [ 3 ] 0 0 0 0> iid_result$p.value24472.5802675822 48945.1605351644 73417.7408027466 97890.3210703288 2 3.824768e-197 9.343533e-133 3.128022e-96 9.606095e-85 3 9.378948e-234 3.883166e-120 1.014683e-76 3.357185e-67

檢驗結果在2歷史、3歷史下p值幾乎均為0，拒絕原假設，故不是獨立同分布，說明很難用線性模型來擬合，這也與之前做條件異方差分析的結果相合。

為了更清楚的說明非線性模型，這里有必要介紹下面幾種：

（一）自激發??模型

?體制表示有??個分段函數，每個分段函數都是??階自回歸??模式的線性表達式，它最大特點是必須要有一個??階延遲的時刻變量??作為門限控制，?是門限參數。其缺點是條件均值方程不連續，畢竟它是分段函數，很容易產生間斷點。

（二）平滑轉移??模型

?實質是對??標準化，然后以它作為自變量，施加一個??連續函數，這個連續函數特點是值域為??的區間，所以適用的形式可以是Sigmoid函數或分布函數，所以相當于分配了一個0到1之間的權重系數，這就是光滑的來歷。其優點是可以直接運用單位根是否在圓內明確時序??的平穩性，缺點是參數??和??很難估計。

（三）神經網絡模型（多層感知機）

每個隱藏層（hidden layer）的一個節點其實都是一個自回歸線性表達式的函數，一般稱激活函數，比如最常見的Sigmoid函數。在機器學習領域，習慣稱??為偏移bias，稱??為權重，它們都是向量或矩陣。

神經網絡模型是半參數模型，激活函數和自回歸表達式已知，只需要算出各層的權重和偏移這些參數即可。計算機方式與數學方式最大區別在于：計算機喜歡用迭代方式解決復雜運算問題，數學喜歡用推理方式解決本質問題。所以最好的解決問題方法就是先理解數學原理，然后用計算機解決。比如在本例中，我們發現日收入是一個非線性問題，用傳統的線性方式預測的精度較低，于是就可以用計算機擬合非線性模型去解決，是一個較好的選擇。

三. LSTM模型深入理解

理解LSTM需要先清晰的理解RNN原理，先看下面這個極簡的例子，是keras發明者在自己的書中通過一個例子闡述RNN原理，它實現了只有一個隱藏層、一個隱層節點的神經網絡：

import numpy as np timesteps = 100 input_features = 32 output_features = 64 inputs = np.random.random((timesteps, input_features)) state_t = np.zeros((output_features,)) W = np.random.random((output_features, input_features)) U = np.random.random((output_features, output_features)) b = np.random.random((output_features,)) successive_outputs = [] for input_t in inputs:output_t = np.tanh(np.dot(W, input_t) + np.dot(U, state_t) + b)successive_outputs.append(output_t)state_t = output_t final_output_sequence = np.stack(successive_outputs, axis=0)

上述程序只是示例單個序列，不難看到下述張量及它們的shape：

inputs：[100, 32]，輸入張量的樣本時序集合，shape[0]=100是時序步長，shape[1]=32是輸入特征維度

state_t：[64, ]，狀態張量，是一個長度為64的array

W：[64, 32]，輸入張量權重，本例中只有正向傳播，沒有更新

U： [64, 64]，狀態張量權重，本例中只有正向傳播，沒有更新

b：[64, ]，偏移bias，也是一個長度為64的array

input_t：[32, ]，輸入張量，是inputs在每一個時刻的輸入

output_t：[64, ]，輸出張量，shape為[64, ]

對輸入張量和狀態張量作運算??，相當于??，運算的張量結果為??，再進行??激活函數，輸出范圍為??的值域區間，作為output_t，將它作為下一時刻的state_t，然后重復迭代下一個時刻的運算，如此循環，直到input_t時序步長結束。此例中，相當于迭代了100次。

針對這個例子我自己畫了一個對應的RNN原理圖，對比可看出LSTM原理圖的雛形。

一個簡單的RNN原理圖

LSTM原理圖

上述LSTM原理圖是一種邏輯關系圖，其實并不有助于理解，因此我按照神經網絡模型的輸入層、隱藏層、輸出層的形式，重新梳理畫了如下圖，會更容易理解來龍去脈及性質：

我們作如下深入解釋：

（1）三個門控??、??、??分別代表遺忘、輸入、輸出，它們都是由當前時刻的輸入數據??和上一時刻的輸出也即hidden state隱藏狀態??，分別作權重??和偏移??的線性組合，然后通過??激活函數形成，它們的值域均為?。門控可以理解成開關系數，接近??表示信息不通過，接近??表示信息通過。

（2）一條信息傳送帶??代表cell state單元狀態，也有翻譯成細胞狀態、神經元狀態，這個cell state實際就是整個 LSTM cell 的狀態變量。其中與之相關的??表示當前信息，它也是由??和??分別作權重??和偏移??的線性組合，不過它是通過??函數形成，值域為?。注意??和??它們都不是門控，它們是??和??兩個門控所作用的載體，?控制上一個信息??，?控制本次輸入的信息??，綜合作用疊加形成新的載體??，所以??。上圖中只有這條信息傳送帶我畫了一條實線，貫穿整個過程。

（3）最終輸出??代表hidden state隱狀態。某種程度上雖然也可以將??作為過程的一個輸出，但更為重要的、真正意義上的輸出只有??，它由輸出門控??所控制。因為??是由兩部分信息迭代累加，很可能超出一定范圍，所以需要用??函數約束一下值域范圍，因此最終的輸出也就是新的隱狀態??。

（4）上述??表示的是向量內積，而不要簡單看成數量乘關系。實際情況下，輸入、狀態、輸出都是向量形式，權重是矩陣形式，可以結合下面（8）參數一起看以便于理解。

（5）??和??性質完全不同，?是激活函數，值域為?，所以完全充當門控的作用，而??函數雖然也常作為激活函數，值域為?，但在此處LSTM中它的作用是控制信息的范圍不會溢出。

（6）信息傳送帶??和輸出??的關系就好比運載火箭和飛船（著陸艙）的關系，運載火箭在整個過程中提供動力，飛船和著陸艙才是我們需要的最終結果，飛船執行太空實驗任務，著陸艙將宇航員和實驗結果帶回地球。

（7）紅色虛線框代表一個LSTM單元，即LSTM Cell，它接受上一時刻的??、??和當前時刻的??作為輸入，經LSTM單元處理后，得到當前時刻的??、??作為輸出，再和下一時刻的??一起再次作為輸入，重復迭代這樣的循環過程，所以說LSTM就是將同一LSTM Cell沿時間展開的RNN，循環迭代多少次呢，就是下面要提到的時間步timesteps，在一個timesteps循環迭代過程中，記憶傳送穿過該條樣本，并得到記錄。

（8）幾個重要參數含義：

（8.1）在Tensorflow里，有一個BasicLSTMCell函數，其調用形式為cell = tf.contrib.rnn.BasicLSTMCell(num_units, state_is_tuple)，返回的是LSTM Cell所對應的兩個輸出??和??：

參數num_units字面解釋為數量單元，有的地方喜歡用size或hidden_size，字面解釋為隱層大小，其實都是一個意思，表示一個LSTM單元（LSTM Cell）的輸出維度，結合上面RNN程序示例看出這個維度和權重??和??、偏移??、狀態??、輸出??的維度都是相一致的，即num_units=64。這里順便提一下共享權重，針對三個門控和一個??，先將?concat連接為一個新的向量，其維度等于輸入維度加上輸出維度，結合上述RNN案例即為32+64=96，然后分別作用于這 4 類共享權重??，對應矩陣維度均為??。

參數state_is_tuple表示是否將兩個內部狀態??和??作為二元組，設置為True表示接受和返回的狀態是c_state和m_state的2-tuples，這是官方建議的。若設置為False表示這兩個內部狀態將按列連接起來返回，盡管目前False是默認值，但它很快要被deprecated。

（8.2）keras進一步封裝成LSTM函數，其調用形式為lstm=keras.layers.recurrent.LSTM(units, input_dim, return_sequences,?stateful)，返回視參數return_sequences而定：

參數units和上面num_units一樣都是表示輸出維度；

參數input_dim表示輸入維度，當使用該層為模型首層時，應指定該值或等價的指定input_shape，例如shape為(sample, timesteps, input_dim) 或 (batch_size, timesteps, input_dim)的3D張量 ；

參數return_sequences表示是否返回每個時間步連續輸出的完整序列，設置為True返回shape為 (batch_size, timesteps, output_dim) 的三維張量，默認設置為False僅返回最后一步的輸出序列，即shape為 (batch_size, output_dim) 的二維張量。該參數用于多個LSTM層之間的傳遞；

參數stateful?表示索引 i 處每個樣本的最后狀態將被用作下一次批處理中索引 i 處樣本的初始狀態，而這對于存在自相關性的時間序列必須設置為True，否則默認值False會在正常（無狀態）模式下，對樣本重新洗牌，時間序列與其滯后項之間的依賴關系將會丟失。結合上面的LSTM模型其實就是??是否在下一個batch重新初始化，對于一個完整的時間序列，我們當然不希望也不能重新初始化，否則這條信息傳送帶將失去應有的作用。

（8.3）如果需要運用多個LSTM Cell，Tensorflow可以用MultiRNNCell(cells)函數，cells是通過BasicLSTMCell實例生成的列表；keras可以用model.add(LSTM(units))來添加多個LSTM層，方法更為簡便靈活。

四. LSTM模型預測（一）：框架說明

我們先梳理一下框架：

1. 數據集：我們在原序列860天數據基礎上，又補充了135天的時序數據，合并作為原序列。

2. 測試集：上述995天數據拆分為兩部分，取20%作為測試集，僅用于預測模型的最終評估。

3. 訓練集和驗證集：995天數據取80%為訓練集和驗證集，這其中再取80%用于訓練，20%給驗證集用于超參數（hyperparameter）調試。

4. 參數stateful?設置：因為我們在ARMA分析時得知，收入時序的ACF呈現7天自相關性，需要長期記憶，故我們需要設置stateful=True。

5. 輸入：輸入的3D張量(batch_size, timesteps, input_dim)設置為(256,3,1)，著重說一下：

input_dim=1：因為只有收入數據，故輸入feature維度為1；

timesteps=3：時間步長的設置應該和序列AR階保持一致，本例中經7步差分后序列的ARMA模型仍存在7天周期的季節AR(1)，以及一個普通的AR(2)或AR(1)，它們之間的作用是交錯的，因此時間步長無法簡單設置，我們經實際調優后設置為3；

batch_size=256：對于存在條件異方差的序列，batch_size如果選的特別小比如1，則模型更新權重參數的頻率很高，振蕩嚴重，效果會很差。經調試后設置為256效果較好，每批次相當于輸入8-9個月的樣本時序數據，結合前面stateful?的設置，模型會在每個batch的對應位置樣本之間傳遞??狀態信息。

6. LSTM層數：鑒于數據量不算大，一層足夠，經調優確實一層效果最好。

7. LSTM狀態單元數量：即輸出維度，決定模型復雜度，經調優設置為?最佳。

8. 訓練迭代次數：調用EarlyStopping在200輪以內穩定，故設置為??。

9. 訓練的參數數量：我們只添加一個LSTM層和一個全鏈接dense（密集）層。

LSTM層：對于每個門控??、??、??及傳送帶狀態??，都有??，??輸入特征維度是??，??輸出維度是??，偏移??維度同輸出一致為??，故所需要訓練的參數共有?，這里要注意的是參數數量和?、??均無關系，但每一次??會重新計算本輪訓練的全部參數，每一批次??內會根據輸入的樣本不斷更新這些參數。

dense層：作為最終輸出，即預測次日的收入，設置的??，它也是經過權值矩陣及偏移的線性計算，再通過激活函數輸出，即：?，?就是上一層LSTM的輸出，其維度為??，??是權重，本層最終輸出為?，因此dense層訓練的參數數量為??。

四. LSTM模型預測（二）：原序列

以上是超參數調優結果，似乎我們可以開始編寫代碼了，但從我之前的文章已獲知收入序列是非平穩序列，一個航空模型的特例——沒有趨勢但有季節特征，如果直接拿這樣的時間序列交給模型作預測，無疑增加它的學習負擔。實踐結果也是如此，經調參后驗證集loss值最好成績為0.05946，相應的測試集RMSE為45154.894，而且在同樣的參數下模型表現不穩定，通常情形RMSE值還要更高，下面這個實驗結果還是我較早之前做的：

RMSE值高達64268.95，loss值一直趨于平緩不再下降，甚至后一段還上升，產生過擬合：

測試集擬合的預測值與真實值效果比對圖：

從圖看出，擬合效果欠佳。所以我們需要采取其它手段進一步優化，比如：進行7步季節差分，同時為避免預測值出現負值，對原序列取對數，這樣做還有一個好處是能減弱序列的異方差性質，增加正態性，在程序結束部分當然還要對序列進行還原操作。

四. LSTM模型預測（三）：原序列取對數、季節差分

下面開始正式構建代碼，我將使用keras框架高效搭建，用Fran?ois Chollet自己的話說，他寫keras的初衷就像樂高玩具一樣，提供一個模型架構，組裝一下就可以了。

讀入數據，先取對數，然后季節差分，最后標準化，這里為了和LSTM中的狀態變量值域保持一致，選擇的是MinMaxScaler函數將值域控制在[-1,1]之間：

train_data = read_file('A店-每日收入指標跟蹤表-訓練.csv') test_data = read_file('A店-每日收入指標跟蹤表-測試.csv')train_data = train_data.iloc[:,1:2] train_data.columns = ['本日收入']test_data = test_data.iloc[:,1:2] test_data.columns = ['本日收入']# 訓練和測試兩個數據集合并 data = concat((train_data,test_data), axis=0)# 取對數 data = np.log(data)# 7天季節差分 diff_data = data.diff(7) diff_data.dropna(inplace=True)# 標準化 scaler = MinMaxScaler(feature_range=(-1, 1)) scaled = scaler.fit_transform(diff_data)

定義超參數，以下是不斷調優后的結果：

# 定義超參數 batch_size = 256 # 批數據大小 epochs = 150 # 實驗訓練次數 rnn_unit = 50 # 狀態的輸出維度 data_dim = 1 # 最終輸出維度，即預測天數 timesteps = 3 # 時間步長

根據timesteps準備訓練數據格式，series_to_supervised()函數按照timesteps提供輸入的序列，然后劃分能被batch整除的數據段，最后按照LSTM輸入格式要求reshape張量。這里需要了解張量變換的重要性，相當于為了模型所需要的維度，去做相應的變換：

# 按batch拆分數據 reframed = series_to_supervised(scaled, timesteps, 1) batch_len = len(reframed) // batch_size batch_end = batch_len*8 // 10 train_scaled, test_scaled = reframed.values[0:batch_size*batch_end, :], reframed.values[batch_size*batch_end:batch_size*batch_len, :] test_X, test_y = split_xy(test_scaled, timesteps, 1) # 訓練集中再分20%驗證集 train_batch_len = len(train_scaled) // batch_size train_batch_end = train_batch_len*8 // 10 train_scaled, val_scaled = train_scaled[0:batch_size*train_batch_end, :], train_scaled[batch_size*train_batch_end:batch_size*train_batch_len:, :] train_X, train_y = split_xy(train_scaled, timesteps, 1) val_X, val_y = split_xy(val_scaled, timesteps, 1) # reshape train_X_3D = train_X.reshape(train_X.shape[0], timesteps, data_dim) val_X_3D = val_X.reshape(val_X.shape[0], timesteps, data_dim) test_X_3D = test_X.reshape(test_X.shape[0], timesteps, data_dim)

開始準備模型，rnn_unit=50，stateful=True，一個LSTM層，一個dense層，這里試過加一個正則化dropout層，但效果沒有提升，所以注釋掉了，回調callbacks_list用于在驗證集上調試，定義loss='mse'和優化器optimizer='rmsprop'：

model = Sequential() model.add(LSTM(rnn_unit, stateful=True, batch_input_shape=(batch_size, timesteps, data_dim))) # model.add(Dropout(0.2)) model.add(Dense(data_dim)) model.summary() callbacks_list = [keras.callbacks.EarlyStopping(monitor='val_loss',patience=2,),keras.callbacks.ModelCheckpoint(filepath='my_model.h5', monitor='val_loss',save_best_only=True,),keras.callbacks.ReduceLROnPlateau(monitor='val_loss', factor = 0.01, patience = 2,)] model.compile(loss='mse', #maeoptimizer='rmsprop') #rmsprop，adam

keras層結構如下，模型不能設置過于復雜，總共需要訓練50*(50+1+1)*4+(50+1)=10451個參數：

_________________________________________________________________ Layer (type) Output Shape Param # ================================================================= lstm_1 (LSTM) (256, 50) 10400 _________________________________________________________________ dense_1 (Dense) (256, 1) 51 ================================================================= Total params: 10,451 Trainable params: 10,451 Non-trainable params: 0 _________________________________________________________________

調優完超參后，在迭代epochs=150次過程中，在第126次后中止：

history = model.fit(train_X_3D, train_y, batch_size=batch_size, epochs=epochs, verbose=1,validation_data=(val_X_3D, val_y), shuffle=False, callbacks=callbacks_list) # 繪畫訓練和測試數據集的損失曲線 plt.plot(history.history['loss'], label='train') plt.plot(history.history['val_loss'], label='validate') plt.legend() plt.show()256/256 [==============================] - 0s 70us/step - loss: 0.0312 - val_loss: 0.0388 Epoch 125/150256/256 [==============================] - 0s 143us/step - loss: 0.0312 - val_loss: 0.0388 Epoch 126/150256/256 [==============================] - 0s 37us/step - loss: 0.0312 - val_loss: 0.0389

訓練集和驗證集損失圖：

開始訓練，喂入數據，循環epochs次，因為之前stateful=True，所以每次新的迭代需要用reset_states()重置LSTM狀態變量：

for i in range(epochs):history = model.fit(train_X_3D, train_y, batch_size=batch_size, epochs=1, verbose=1,validation_data=(val_X_3D, val_y), shuffle=False, callbacks=callbacks_list)model.reset_states()

讓我們看看訓練的效果，150次迭代后，驗證集上的loss值穩定在0.0379左右：

cost = model.evaluate(val_X_3D, val_y, batch_size=batch_size) print('validate cost:', cost)256/256 [==============================] - 0s 33us/step - loss: 0.0309 - val_loss: 0.0378 Train on 256 samples, validate on 256 samples Epoch 1/1256/256 [==============================] - 0s 31us/step - loss: 0.0309 - val_loss: 0.0376 Train on 256 samples, validate on 256 samples Epoch 1/1256/256 [==============================] - 0s 27us/step - loss: 0.0309 - val_loss: 0.0378256/256 [==============================] - 0s 8us/step validate cost: 0.03793208673596382

前面訓練、驗證、確定超參后，一旦運行測試集，就不能再根據運行結果作修改以防信息泄露。下面開始用這個model直接進行預測，按batch送入數據，按照之前數據處理順序的逆過程，先還原標準化，再還原差分，最后還原對數：

# 預測 forecasts = list() test_batch = batch_len - batch_end for i in range(test_batch):test_batch_X, test_batch_y = test_X[i*batch_size:(i+1)*batch_size, :], test_y[i*batch_size:(i+1)*batch_size, :]predict_input = test_batch_X.reshape(batch_size, timesteps, data_dim)forecast = model.predict(predict_input, batch_size=batch_size)for j in range(len(forecast)): forecasts.append(forecast[j,:]) # store the forecast forecasts = np.array(forecasts)# 定義還原為收入數值的函數 def diff_inverse(series, pred_data, scaler):inverted = list()for i in range(len(pred_data)):# create array from forecastforecast = array(pred_data[i])fc = forecast.reshape(len(forecast),1)# invert scalinginv_scale = scaler.inverse_transform(fc)# 當日預測差分數據 + 7天前的原始數據 = 當日預測收入if len(forecast) == 1:index = batch_end * batch_size + i + timestepslast_diff = series.values[index]inverted.append(inv_scale[0, :] + last_diff)else:for j in range(len(forecast)):index = batch_end*batch_size + i*batch_size + j + timestepslast_diff = series.values[index]inverted.append(inv_scale[j,:] + last_diff)return inverted# 還原真實 real = diff_inverse(data, test_y, scaler)# 還原預測 predicted = diff_inverse(data, forecasts, scaler) predicted = np.array(predicted)# 對數還原 real = np.exp(real) predicted = np.exp(predicted)

計算測試集上的RMSE均方根誤差，其值為36817.065：

# calculate RMSE rmse = sqrt(mean_squared_error(real, predicted)) print('Test RMSE: %.3f' % rmse)Test RMSE: 36817.065

保存文件，并繪制預測效果圖：

# 保存預測結果 result = concat((DataFrame(real), DataFrame(predicted)), axis=1) result = DataFrame(result) result.columns = ['測試集收入真實值','測試集收入預測值'] # plt.plot(result) # plt.show() result.to_csv('output\\LSTM_predict_result.csv', encoding='utf_8_sig')# 繪圖 plt.rcParams['font.sans-serif'] = ['SimHei'] # 用來正常顯示中文標簽 plt.rcParams['axes.unicode_minus'] = False # 用來正常顯示負號fig = plt.figure() fig.add_subplot() plt.plot(real, color='blue', label='真實值') plt.plot(predicted, color='red', label='預測值') plt.legend(loc='best') plt.title('測試集預測效果比對') plt.show()

圖像如下，藍色是真實值，紅色是預測值，大部分情況擬合的很不錯，部分擬合的不太好：

最后看數據，基本預測的都還不錯，就是節日預測的很不好。我將測試集當日收入的真實值、指標值和預測值做了比對，預測數據比原先手工分解指標的效果要好：

最后，讓我們看看模型的穩定性，連續運行幾次：

validate cost: 0.037923164665699005 Test RMSE: 35885.605

validate cost: 0.03791944310069084 Test RMSE: 36851.139

畢竟這家門店的日收入標準差為49799，如果日均3萬多元的誤差在接受范圍內，那么就可以將這個模型的預測值作為日收入指標，即門店的次日收入期望，當然這個模型本身我認為還有進一步優化的余地。

四. LSTM模型預測（四）：增加特征維度 - 日期屬性

既然模型擬合節日效果不好，我試圖增加日期屬性，程序沒有太大變動，只是利用OneHot將日期屬性進行編碼，并作為輸入特征：

# 特征“日期屬性”離散化 data = pd.get_dummies(data,columns=['日期屬性'])

OneHot編碼效果如下圖所示：

簡單調參后，設置batch_size = 128、timesteps = 7、input_dim = 4，模型在50輪左右就中止了：

設置epochs = 50，訓練速度很快，一兩秒后得出結果：

128/512 [======>.......................] - ETA: 0s - loss: 0.0145 512/512 [==============================] - 0s 83us/step - loss: 0.0295 - val_loss: 0.0322 Train on 512 samples, validate on 128 samples Epoch 1/1128/512 [======>.......................] - ETA: 0s - loss: 0.0145 512/512 [==============================] - 0s 44us/step - loss: 0.0295 - val_loss: 0.0324128/128 [==============================] - 0s 12us/step validate cost: 0.03258267790079117 Test RMSE: 35737.020

驗證集loss值為0.03258，測試集的RMSE為35737，沒有太多改觀。比對圖如下：

預測比對結果如下，節日預測并沒有很大改觀，說明增加日期屬性的特征并不能提升次日收入的預測能力：

為測試模型的泛化能力，調整測試集比例占比為10%后，驗證集loss值減少為0.02877，但測試集的RMSE增加為49980.582：

調整測試集比例占比為30%后，驗證集loss值為0.02935，測試集的RMSE值為37677.604：

四. LSTM模型預測（五）：增加特征維度 - 增長率

這是我之前做過的實驗，源于朋友Adam曾經給過我的建議，提示添加諸如增長率等新的特征會有助于模型預測，于是我在“本日收入”維度之外，增加了一個收入增長率的人造特征，表達預測趨勢。這個??時刻的增長率用??表示，實際上有?，數據如下：

有人可能會問，按正常邏輯上面第一行增長率1.657應該出現在第二行上。但如果這樣的話，當前增長率數據表達的是對前一天的增長趨勢，而我們是要充分獲取輸入信息來預測后一天的收入，趨勢正好弄反了，對我們預測毫無益處。

先繪制收入和增長率合在一起的圖形：

下面是分開顯示的效果圖：

同樣的模型訓練數據，經過在測試集上的驗證，得到的效果非常好，RMSE值只有17474.04，loss值一直趨于下降，最終平緩至接近訓練集誤差水平。

測試集上預測值與真實值擬合的效果圖如下：

看的出來，擬合效果接近完美。為什么？請參見我后面的數學解釋。

四. LSTM模型預測（六）：增加高相關系數的特征

這里簡單提及一下還是我早些時候做過的工作：結合之前線性相關性分析結果，在保持本日收入和增長率兩個特征的基礎上，又增加了本日園內自營游樂收入、園內游樂合計、園內游樂自營、耗幣收入、耗幣數量這5個相關系數為0.9以上的特征，事實證明擬合效果較差，和僅保留本日收入一個維度的結果差不多，RMSE為55888.902。

這充分說明，維度特征的選擇對于擬合結果至關重要，即便在有增長率的情形下，由于增加了對時序預測無關的特征，只會增加數據冗余（噪聲）和復雜度，預測效果很可能不升反降。

五. 兩種模型的數學解釋

我們分析一下上述兩種模型（三）和（五），用??表示輸出（預測），用??、?... 表示不同維度的輸入（特征）：

第一種模型（三）：?

函數??顯然是未知的，對于一維輸入收入??，模型需要在一個線性關系上采用激活函數來近似擬合這個即便存在也極為復雜的函數關系。根據前章闡述的結論，比如我們可以用ARIMA??模型來擬合，那么這個函數關系??可以通過一個隱函數??間接表示為?，其中??，，，??還是一個存在條件異方差的殘差序列，因此這個模型挑戰難度非常之大。從訓練結果來看，在迭代超過100次后，梯度不再下降反而有所反彈，試過增加patience和dropout正則化均沒有改觀，說明很可能由于??的復雜性，僅在局部擬合的還不錯，就好像模型走進入了一個深淺毫無規律的海底世界，它只沿某一梯度下降到一個局部低洼處，但究竟哪里才是最深的位置？這種探索難度可想而知。

第二種模型（五）：?

函數??在這里是顯式明確的，因為對于二維輸入收入??和增長率??，我們很容易得出??，這個三維圖形長如下樣子：

為便于看清這個三維圖形，我調整角度旋轉了一下，這個曲面很像一個馬鞍形狀，取值范圍也很客觀地從樣本數據提取。有了這樣的顯式函數內在形式，模型會容易找到梯度方向，從任何一點都很容易走到馬鞍的中心，所以訓練起來很輕松。我們的二維數據實際可以看作是上述這個二元函數在一定范圍的樣本取樣，因此用LSTM模型擬合的效果奇好。

類似還有預測正弦函數，因為函數都是顯式的，通過在一個很小的區間上（timesteps時序）取樣，正弦光滑曲線某一段幾乎都可以用一個線性線段近似，所以擬合效果也非常好。綜上，對于存在簡單顯式關系的函數形式，或者說對于簡單光滑的曲線或曲面，梯度下降模型通常都會表現的很好。

但問題來了，第二種模型的預測由于最后一天缺少次日收入（也正是我們要預測的），因此增長率無法獲知，所以最后一天的次日收入仍然是無法預測的。如果要達到預測目的，就必須要先行預測出次日增長率，這當然也是一個思路，對于這家門店增長率也存在7天自相關性，因此增長率的預測同樣存在挑戰。

六. 總結

1. 從模型表現來看，ARMA和LSTM旗鼓相當，其實ARMA更適合探尋數據內在規律本質，LSTM更適合做靈活、廣泛的預測，包括線性的、非線性的，輸入特征也可以是數值的、非數值的。

2. ARMA和LSTM都只能做短期預測，它們內部都存在最基本的線性擬合原理，嚴格來說，它們更適合處理平穩時序，盡管對非平穩時序也有類似差分這樣相應的措施，但仍然無法完全解釋清楚這些不規則擾動的規律。

3. 對于存在自相關性的時序，使用stateful能優化效果，從驗證集loss值來看，能提升8%左右。

總的來說，在數據集規模較小的情況下（不到1000條），獲得以上預測效果已經算不錯，可以從以下三個方面進行優化：

1. 預測增長率，如果增長率預測效果好，可直接計算次日收入。

2. 增加對預測有直接幫助的特征，比如次日天氣預報、次日日期屬性。

3. 更為全面、深入地調優超參數。

總結

以上是生活随笔為你收集整理的收入时间序列——之预测总结篇的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。