理论加实践，终于把时间序列预测ARIMA模型讲明白了
上篇我們一起學習了一些關于時間序列預測的知識。而本文將通過一段時間內電力負荷波動的數據集來實戰演示完整的ARIMA模型的建模及參數選擇過程，其中包括數據準備、隨機性、穩定性檢驗。本文旨在實踐中學習，在實戰過程中穿插理論知識梳理和學習，相信大家一定有所收獲。
本文主要內容
時間序列建模基本步驟
獲取被觀測系統時間序列數據。
對數據繪圖，觀測是否為平穩時間序列；對于非平穩時間序列要先進行??階差分運算，化為平穩時間序列。
經過第二步處理，已經得到平穩時間序列。要對平穩時間序列分別求得其自相關系數ACF 和偏自相關系數PACF ，通過對自相關圖和偏自相關圖的分析，得到最佳的階層??和階數?。
由以上得到的?、、，得到ARIMA模型。然后開始對得到的模型進行模型檢驗。
ARIMA模型介紹
ARIMA 模型[1]是一種流行且廣泛使用的時間序列預測統計方法。
ARIMA 是代表autoRegressive I integrated Moving a average[2]自回歸綜合移動平均線的首字母縮寫詞，它是一類在時間序列數據中捕獲一組不同標準時間結構的模型。預測方程中平穩序列的滯后稱為“自回歸”項，預測誤差的滯后稱為“移動平均”項，需要差分才能使其平穩的時間序列被稱為平穩序列的“綜合”版本。隨機游走和隨機趨勢模型、自回歸模型和指數平滑模型都是 ARIMA 模型的特例。
ARIMA 模型可以被視為一個“過濾器”，它試圖將信號與噪聲分開，然后將信號外推到未來以獲得預測。ARIMA模型特別適合于擬合顯示非平穩性的數據。
一般概念
為了能夠使用ARIMA，你需要了解一些概念。
平穩性
從統計學的角度來看，平穩性是指數據的分布在時間上平移時不發生變化。因此，非平穩數據顯示了由于趨勢而產生的波動，必須對其進行轉換才能進行分析。例如，季節性會導致數據的波動，并可以通過“季節性差異”過程消除。
差分
從統計學的角度來看，數據差分是指將非平穩數據轉換為平穩的過程，去除其非恒定的趨勢。“差分消除了時間序列水平的變化，消除了趨勢和季節性，從而穩定了時間序列的平均值。”??季節性差分應用于季節性時間序列以去除季節性成分。
ARIMA模型拆解
剖析ARIMA的各個部分，以便更好地理解它如何幫助我們時間序列建模，并對其進行預測。
AR - 自回歸
自回歸模型，顧名思義，就是及時地“回顧”過去，分析數據中先前的值，并對它們做出假設。這些先前的值稱為“滯后”。一個例子是顯示每月鉛筆銷售的數據。每個月的銷售總額將被認為是數據集中的一個“進化變量”。這個模型是作為“利益的演化變量根據其自身的滯后值(即先驗值)進行回歸”而建立的。
I - 表示綜合
與類似的“ARMA”模型相反，ARIMA中的“I”指的是它的綜合方面。當應用差分步驟時，數據是“綜合”的，以消除非平穩性。表示原始觀測值的差異，以允許時間序列變得平穩，即數據值被數據值和以前的值之間的差異替換。
MA - 移動平均線
該模型的移動平均方面，是將觀測值與應用于滯后觀測值的移動平均模型的殘差之間的相關性合并。
ARIMA用于使模型盡可能地符合時間序列數據的特殊形式。
ARIMA模型建立
一般步驟
① 首先需要對觀測值序列進行平穩性檢測，如果不平穩，則對其進行差分運算直到差分后的數據平穩；
② 在數據平穩后則對其進行白噪聲檢驗，白噪聲是指零均值常方差的隨機平穩序列；
③ 如果是平穩非白噪聲序列就計算ACF（自相關系數）、PACF（偏自相關系數），進行ARMA等模型識別；
④ 對已識別好的模型，確定模型參數，最后應用預測并進行誤差分析。
一般地，對于給定的時間序列?，平穩序列的建模過程可以用下圖中的流程圖表示。
ARIMA實戰剖析
導入必要的庫
導入statmodelsPython庫已使用ARIMA模型。
import?os import?warnings import?matplotlib.pyplot?as?plt import?numpy?as?np import?pandas?as?pd import?datetime?as?dt import?mathfrom?pandas.plotting?import?autocorrelation_plot from?statsmodels.tsa.statespace.sarimax?import?SARIMAX from?sklearn.preprocessing?import?MinMaxScaler from?common.utils?import?load_data,?mape from?IPython.display?import?Imagefrom?statsmodels.graphics.tsaplots?import?plot_acf,?plot_pacf?? from?statsmodels.tsa.stattools?import?adfuller??#?adf檢驗庫 from?statsmodels.stats.diagnostic?import?acorr_ljungbox??#?隨機性檢驗庫 from?statsmodels.tsa.arima_model?import?ARMA?%matplotlib?inline plt.rcParams['figure.figsize']?=?(12,6) pd.options.display.float_format?=?'{:,.2f}'.format np.set_printoptions(precision=2) warnings.filterwarnings("ignore")?#?specify?to?ignore?warning?messages導入數據
energy?=?pd.read_csv('./data/energy.csv') energy.head(10)繪制從2012年1月到2014年12月的所有可用能源數據。看到這些數據，并不陌生，因為在之前的文章中已經展示了部分數據。
energy.plot(y='load',?subplots=True,?figsize=(15,?8),?fontsize=12) plt.xlabel('timestamp',?fontsize=12) plt.ylabel('load',?fontsize=12) plt.show()劃分訓練和測試數據集
現在已經加載了數據，可以將其劃分為訓練集和測試集。要在訓練集上訓練模型。通常，在模型完成訓練后，將使用測試集評估它的準確性。需要確保測試集涵蓋了來自訓練集的較晚時間段，以確保模型不會從未來時間段獲取信息。
從2014年9月1日到10月31日，分配兩個月的時間給訓練集。測試集將包括2014年11月1日至12月31日兩個月的時間段:
train_start_dt?=?'2014-11-01?00:00:00' test_start_dt?=?'2014-12-30?00:00:00'由于這一數據反映的是每日能源消費，因此存在強烈的季節性模式，但當前消費與最近幾天的消費規律最為相似。
可視化差異
為了更加直觀地看出訓練集和測試集的差異，我們在同一張圖中用不同顏色區分兩個測試集，藍色為訓練集、橙色為測試集。
energy[(energy.index?<?test_start_dt)?&?(energy.index?>=?train_start_dt)][['load']].rename(columns={'load':'train'})?\.join(energy[test_start_dt:][['load']].rename(columns={'load':'test'}),?how='outer')?\.plot(y=['train',?'test'],?figsize=(15,?8),?fontsize=12) plt.xlabel('timestamp',?fontsize=12) plt.ylabel('load',?fontsize=12) plt.show()使用一個相對較小的時間窗口來訓練數據就足夠了。
準備訓練數據
現在需要通過對數據進行篩選和歸一化來為模型訓練準備數據。篩選需要的時間段和列的數據，并且對其進行歸一化，其作用的是將數據投影在0-1之間。
① 過濾原始數據集，只包括前面提到的每個set的時間段，只包括所需的列'load'加上日期索引。
train?=?energy.copy()[(energy.index?>=?train_start_dt)?&?(energy.index?<?test_start_dt)][['load']] test?=?energy.copy()[energy.index?>=?test_start_dt][['load']]print('Training?data?shape:?',?train.shape) print('Test?data?shape:?',?test.shape)Training data shape: (1416, 1) Test data shape: (48, 1)② 使用MinMaxScaler()對訓練數據進行 (0, 1) 標準化。
scaler?=?MinMaxScaler() train['load']?=?scaler.fit_transform(train) train.head(10)③ 原始數據和標準化數據進行可視化比較。
energy[(energy.index?>=?train_start_dt)?&?(energy.index?<?test_start_dt)][['load']].rename(columns={'load':'original?load'}).plot.hist(bins=100,?fontsize=12) train.rename(columns={'load':'scaled?load'}).plot.hist(bins=100,?fontsize=12) plt.show()④ 根據訓練好的歸一化模型，對測試集數據歸一化。
test['load']?=?scaler.transform(test) test.head()穩定性檢驗
adfuller(Augmented Dickey-Fuller)測試可用于在存在串行相關的情況下在單變量過程中測試單位根。
statsmodels.tsa.stattools.adfuller(x， maxlag?=?None，regression?='c'，autolag?='AIC'， store?=?False，regresults?=?False?)adfuller中可進行adf校驗，一般傳入一個data就行，包括?list, numpy array 和 pandas series都可以作為輸入，其他參數可以保留默認。
返回值：
adf（float）
測試統計
pvalue（float）
MacKinnon基于MacKinnon的近似p值（1994年，2010年）
usedlag（int）
使用的滯后數量
nobs（int）
用于ADF回歸的觀察數和臨界值的計算
critical values（dict）
測試統計數據的臨界值為1％，5％和10％。基于MacKinnon（2010）
icbest（float）
如果autolag不是None，則最大化信息標準。
resstore?（ResultStore，可選）
一個虛擬類，其結果作為屬性附加
如何確定該序列能否平穩呢？主要看：
1%、%5、%10不同程度拒絕原假設的統計值和ADF Test result的比較，ADF Test result同時小于1%、5%、10%即說明非常好地拒絕該假設。另外，P-value是否非常接近0，接近0，則是平穩的，否則，不平穩。
若不平穩，則需要進行差分，差分后再進行檢測。
def?adf_val(ts,?ts_title):'''ts:?時間序列數據，Series類型ts_title:?時間序列圖的標題名稱，字符串'''?????#?穩定性（ADF）檢驗adf,?pvalue,?usedlag,?nobs,?critical_values,?icbest?=?adfuller(ts)??name?=?['adf',?'pvalue',?'usedlag','nobs',?'critical_values',?'icbest']????values?=?[adf,?pvalue,?usedlag,?nobs,?critical_values,?icbest]??print(list(zip(name,?values)))return?adf,?pvalue,?critical_values,??#?返回adf值、adf的p值、三種狀態的檢驗值用上面定義的函數進行平穩性檢驗。
#?讀取數據 ts_data?=?df['load'].astype('float32')?? adf,?pvalue1,?critical_values?=?adf_val(ts_data,?'raw?time?series')[('adf', -10.404080285485218), ('pvalue', 1.876514522339643e-18), ('usedlag', 49), ('nobs', 26254), ('critical_values', {'1%': -3.430599102593299, '5%': -2.8616500960359854, '10%': -2.5668286008605627}), ('icbest', 265656.2951464001)]adf結果為-10.4， 小于三個level的統計值。pvalue也是接近于0 的，所以是平穩的。
白噪聲檢測
白噪聲檢驗也稱為純隨機性檢驗，當數據是純隨機數據時，再對數據進行分析就沒有任何意義了，所以拿到數據后最好對數據進行一個純隨機性檢驗。
#?數據的純隨機性檢驗函數 acorr_ljungbox(x,?lags=None,?boxpierce=False,?model_df=0,?period=None,?return_df=True,?auto_lag=False)主要參數
lags為延遲期數，如果為整數，則是包含在內的延遲期數，如果是一個列表或數組，那么所有時滯都包含在列表中最大的時滯中。
boxpierce為True時表示除開返回LB統計量還會返回Box和Pierce的Q統計量
返回值
lbvalue:?(float or array)
測試的統計量
pvalue:?(float or array)
基于卡方分布的p統計量
bpvalue:?((optionsal), float or array)
基于 Box-Pierce 的檢驗的p統計量
bppvalue:?((optional), float or array)
基于卡方分布下的Box-Pierce檢驗的p統計量
若p值遠小于0.01，因此我們拒絕原假設，認為該時間序列是平穩的。（這里原假設是存在單位根，即時間序列為非平穩的。)
def?acorr_val(ts):'''#?白噪聲（隨機性）檢驗ts:?時間序列數據，Series類型返回白噪聲檢驗的P值'''lbvalue,?pvalue?=?acorr_ljungbox(ts,?lags=1)??#?白噪聲檢驗結果return?lbvalue,?pvalueacorr_val(ts_data)24056.19, 0.ARIMA模型實現
可以使用statsmodels?庫創建 ARIMA 模型。并遵循以下幾個步驟。
通過調用SARIMAX()并傳入模型參數：?p, d, q參數，以及?P, D, Q參數定義模型。
通過調用fit()函數為訓練數據準備模型。
通過調用forecast()函數進行預測，并指定要預測的步驟數(horizon)。
在ARIMA模型中，有3個參數用于幫助對時間序列的主要方面進行建模：季節性、趨勢和噪聲。
p：與模型的自回歸方面相關的參數，模型中包含的滯后觀測數，也稱為滯后階數。
d：與模型集成部分相關的參數，原始觀測值差異的次數，也稱為差異度。它影響到應用于時間序列的差分的數量。
q：與模型的移動平均部分相關的參數。移動平均窗口的大小，也稱為移動平均的階數。
值 0 可用于參數，表示不使用模型的該元素。這樣，ARIMA 模型可以配置為執行 ARMA 模型的功能，甚至是簡單的 AR、I 或 MA 模型。
Note:?如果數據具有季節性——我們使用季節性ARIMA模型(SARIMA)。在這種情況下，您需要使用另一組參數:' P '， ' D '和' Q '，它們描述了與' p '， ' d '和' q '相同的關聯，不同的是對應于模型的季節性成分。
確定時間序列的差分?
ARIMA 模型對時間序列的要求是平穩型。因此，當你得到一個非平穩的時間序列時，首先要做的即是做時間序列的差分，直到得到一個平穩時間序列。如果你對時間序列做??次差分才能得到一個平穩序列，那么可以使用?模型，其中??是差分次數。
fig?=?plt.figure(figsize=(20,16)) ax1=?fig.add_subplot(211) diff1?=?train.diff(1) diff1.plot(ax=ax1) ax2=?fig.add_subplot(212) diff2?=?train.diff(2) diff2.plot(ax=ax2)可以看出一階差分的時間序列的均值和方差已經基本平穩，二階差分后的時間序列與一階差分相差不大，并且二者隨著時間推移，時間序列的均值和方差保持不變。因此可以將差分次數??設置為1。
確定合適的?
現在我們已經得到一個平穩的時間序列，接來下就是選擇合適的ARIMA模型，即ARIMA模型中合適的?。
模式識別
可通過下面的代碼，計算自相關系數（Autocorrelation Function, SAF）和偏自相關系數（Partial Autocorrelation Function, PACF）。繪制并檢查平穩時間序列的自相關圖和偏自相關圖。
自相關（Autocorrelation）：?對一個時間序列，現在值與其過去值的相關性。如果相關性為正，則說明現有趨勢將繼續保持。
偏自相關（Partial Autocorrelation）：?可以度量現在值與過去值更純正的相關性。
比如，當我們計算??與??的相關性時，?可能會受到??的影響，同時??也會受到??的影響。而偏自相關就是用來計算剔除??影響后，?與??的相關性。
偏自相關的通俗計算過程：
有三個自變量?、、，一個因變量?：
線性回歸建模：通過??和??預測?，取殘差
線性回歸建模：通過??和??預測?，取殘差
由于以上兩個殘差都剔除了??和??的影響，因此對兩個殘差取相關性就是??與??的偏自相關
如果一個時間序列滿足以下兩個條件:
?具有拖尾性，即??不會在??大于某個常數之后就恒等于 0。
?具有截尾性，即??在??時變為 0。
第 2 個條件還可以用來確定階數?。考慮到存在隨機誤差的存在，因此??在??階延遲后未必嚴格為 0 ，而是在 0 附近的小范圍內波動。具體來說，設??階偏自相關系數為?，若階數大于??大部分的偏自相關系數滿足下式，則模型的階數取?。
其中??表示樣本序列長度。
ACF 和 PACF 圖：?通過差分對時間序列進行平穩化后，擬合 ARIMA 模型的下一步是確定是否需要 AR 或 MA 項來校正差分序列中剩余的任何自相關。結合自相關圖和偏自相關圖共同進行判斷時間序列模型。
關于ARMA通用判斷標準說明如下表格：
| AR(p) | 拖尾 | p階截尾 |
| MA(q) | q階截尾 | 拖尾 |
| ARMA(p,q) | 拖尾 | 拖尾 |
| 模型不適合 | 截尾 | 截尾 |
拖尾和截尾說明如下：
拖尾：?始終有非零取值，不會在大于某階后就快速趨近于0（而是在0附近波動），可簡單理解為無論如何都不會為0，而是在某階之后在0附近隨機變化。
截尾：?在大于某階(k)后快速趨于0為k階截尾，可簡單理解為從某階之后直接就變為0。
通常情況下:
如果說自相關圖拖尾，并且偏自相關圖在p階截尾時，此模型應該為AR(p)。
如果說自相關圖在q階截尾并且偏自相關圖拖尾時，此模型應該為MA(q)。
如果說自相關圖和偏自相關圖均顯示為拖尾，那么可結合ACF圖中最顯著的階數作為q值，選擇PACF中最顯著的階數作為p值，最終建立ARMA(p,q)模型。
如果說自相關圖和偏自相關圖均顯示為截尾，那么說明不適合建立ARMA模型。
然后根據如下常用準則選擇模型：
?赤池信息量 akaike information criterion
?貝葉斯信息量 bayesian information criterion
?hannan-quinn criterion
具體方法可以參考：ACF 和 PACF?[6]
手動選擇超參數
為ARIMA模型的參數選擇最佳值可能是一個挑戰，因為這有點主觀，也有點耗時。可以考慮使用'pyramid'庫[7]?中的?'auto_arima()'?函數。文末提供一種網格搜索方法來自動選擇超參數。
本文通過手動選擇參數的方式，也許模型效果不是很理想，目的是進行快速演示建模過程。
① 首先設置horizon值。先試試3個小時:
#?指定要提前預測的步驟數 HORIZON?=?3 print('Forecasting?horizon:',?HORIZON,?'hours')Forecasting horizon: 3 hours② 現在嘗試一些手動選擇參數來找到一個相對好的模型。
order?=?(4,?1,?0) seasonal_order?=?(1,?1,?0,?24)model?=?SARIMAX(endog=train,?order=order,?seasonal_order=seasonal_order) results?=?model.fit()print(results.summary())打印結果。
現在已經建立了一個時序模型，現在我們需要找到一種方法來計算它。
模型評估
為了評估模型，可以使用walk forward驗證。在實踐中，每次有新的數據可用時，時間序列模型都要重新訓練。這使得模型可以在每個時間步驟中做出最好的預測。
從使用該技術的時間序列的開始，在訓練數據集上訓練模型。然后對下一個時間步驟進行預測。根據已知值對預測進行評估。然后將訓練集擴展到包含已知值，并重復該過程。
Note:?為了更有效的訓練，應該保持訓練集窗口固定，以便每次向訓練集添加新的觀察值時，并將該觀察值從集合的開始處刪除。
這個過程為模型在實踐執行提供了更可靠的估計。然而，這是以創建眾多模型的計算成本為代價的。如果數據較小或模型簡單的話，這是可以接受的，但如果數據量大，或者模型規模大可能是一個問題。
Walk-forward validation 是時間序列模型評估的黃金標準，可以考慮用于你自己的項目。
① 首先，為每個HORIZON步驟創建一個測試數據點。
test_shifted?=?test.copy()for?t?in?range(1,?HORIZON+1):test_shifted['load+'+str(t)]?=?test_shifted['load'].shift(-t,?freq='H')test_shifted?=?test_shifted.dropna(how='any') test_shifted.head(5)數據根據它的地平線上點水平移動。
② 在循環中使用滑動窗口方法預測測試數據的長度
%%time training_window?=?720?#?投入30天(720小時)進行訓練train_ts?=?train['load'] test_ts?=?test_shiftedhistory?=?[x?for?x?in?train_ts] history?=?history[(-training_window):]predictions?=?list()order?=?(2,?1,?0) seasonal_order?=?(1,?1,?0,?24)for?t?in?range(test_ts.shape[0]):model?=?SARIMAX(endog=history,?order=order,?seasonal_order=seasonal_order)model_fit?=?model.fit()yhat?=?model_fit.forecast(steps?=?HORIZON)predictions.append(yhat)obs?=?list(test_ts.iloc[t])#?move?the?training?windowhistory.append(obs[0])history.pop(0)print(test_ts.index[t])print(t+1,?':?predicted?=',?yhat,?'expected?=',?obs)我們可以看出訓練的過程
2014-12-30 00:00:00 1 : predicted = [0.32 0.29 0.28] expected = [0.32945389435989236, 0.2900626678603402, 0.2739480752014323]2014-12-30 01:00:00 2 : predicted = [0.3 0.29 0.3 ] expected = [0.2900626678603402, 0.2739480752014323, 0.26812891674127126]2014-12-30 02:00:00 3 : predicted = [0.27 0.28 0.32] expected = [0.2739480752014323, 0.26812891674127126, 0.3025962399283795]③ 將預測結果與實際負荷進行比較:
eval_df?=?pd.DataFrame(predictions,?columns=['t+'+str(t)?for?t?in?range(1,?HORIZON+1)]) eval_df['timestamp']?=?test.index[0:len(test.index)-HORIZON+1] eval_df?=?pd.melt(eval_df,?id_vars='timestamp',?value_name='prediction',?var_name='h') eval_df['actual']?=?np.array(np.transpose(test_ts)).ravel() eval_df[['prediction',?'actual']]?=?scaler.inverse_transform(eval_df[['prediction',?'actual']]) eval_df.head()觀察每小時數據的預測，并與實際負載進行比較。這有多準確?評估模型的準確性
通過測試所有預測的平均絕對百分比誤差(MAPE)來評估模型的準確性。
MAPE是在一個預測方法的預測精度的測量統計。由上述公式定義。實際和預測的差除以實際。“這個計算的絕對值是對每個預測時間點求和，然后除以擬合點的數目n。”?wikipedia[8]
① 用代碼表示方程:
if(HORIZON?>?1):eval_df['APE']?=?(eval_df['prediction']?-?eval_df['actual']).abs()?/?eval_df['actual']print(eval_df.groupby('h')['APE'].mean())② 計算一步的MAPE:
print('One?step?forecast?MAPE:?',?(mape(eval_df[eval_df['h']?==?'t+1']['prediction'],?eval_df[eval_df['h']?==?'t+1']['actual']))*100,?'%')One step forecast MAPE: 0.5570581332313952 %③ 打印多步預測MAPE:
print('Multi-step?forecast?MAPE:?',?mape(eval_df['prediction'],?eval_df['actual'])*100,?'%')Multi-step forecast MAPE: 1.1460048657704118 %結果值較低是很好的：考慮到一個MAPE為10的預測會下降10%。
④ 為更加容易直觀地看到這種精度測量，把他們可視化出來。
if(HORIZON?==?1):##?Plotting?single?step?forecasteval_df.plot(x='timestamp',?y=['actual',?'prediction'],?style=['r',?'b'],?figsize=(15,?8))else:##?Plotting?multi?step?forecastplot_df?=?eval_df[(eval_df.h=='t+1')][['timestamp',?'actual']]for?t?in?range(1,?HORIZON+1):plot_df['t+'+str(t)]?=?eval_df[(eval_df.h=='t+'+str(t))]['prediction'].valuesfig?=?plt.figure(figsize=(15,?8))ax?=?plt.plot(plot_df['timestamp'],?plot_df['actual'],?color='red',?linewidth=4.0)ax?=?fig.add_subplot(111)for?t?in?range(1,?HORIZON+1):x?=?plot_df['timestamp'][(t-1):]y?=?plot_df['t+'+str(t)][0:len(x)]ax.plot(x,?y,?color='blue',?linewidth=4*math.pow(.9,t),?alpha=math.pow(0.8,t))ax.legend(loc='best')plt.xlabel('timestamp',?fontsize=12) plt.ylabel('load',?fontsize=12) plt.show()綜上所述，這個過程的步驟如下：
模型識別。使用繪圖和匯總統計來識別趨勢、季節性和自回歸元素，以了解所需的差異量和滯后大小。
參數估計。使用擬合程序找到回歸模型的系數。
模型檢查。使用殘差的繪圖和統計檢驗來確定模型未捕獲的時間結構的數量和類型。
重復該過程，直到在樣本內或樣本外觀察（例如訓練或測試數據集）上達到理想的擬合水平。
網格搜索選擇超參數
將網格搜索定義為一個函數evaluate_arima_model()，該函數以時間序列數據集作為輸入，以及元組(p,d,q)作為參數用于評估模型。
數據集分為兩部分：初始訓練數據集為 66％，測試數據集為剩余的 34％。
迭代測試集的每個時間步。一次迭代就可以訓練一個模型，然后使用該模型對新數據進行預測。每次迭代都進行預測并存儲在列表中。最后用測試集將所有預測值與預期值列表進行比較，并計算并返回均方誤差分數。
#?evaluate?an?ARIMA?model?for?a?given?order?(p,d,q) def?evaluate_arima_model(X,?arima_order):#?prepare?training?datasettrain_size?=?int(len(X)?*?0.66)train,?test?=?X[0:train_size],?X[train_size:]history?=?[x?for?x?in?train]#?make?predictionspredictions?=?list()for?t?in?range(len(test)):model?=?ARIMA(history,?order=arima_order)model_fit?=?model.fit(disp=0)yhat?=?model_fit.forecast()[0]predictions.append(yhat)history.append(test[t])#?calculate?out?of?sample?errorerror?=?mean_squared_error(test,?predictions)return?error繼續定義一個evaluate_models()的函數，該函數為ARIMA指定(p,d,q)參數，并以網格循環到方式進行迭代。
#?evaluate?combinations?of?p,?d?and?q?values?for?an?ARIMA?model def?evaluate_models(dataset,?p_values,?d_values,?q_values):#?確保輸入數據是浮點值（而不是整數或字符串）dataset?=?dataset.astype('float32')best_score,?best_cfg?=?float("inf"),?Nonefor?p?in?p_values:for?d?in?d_values:for?q?in?q_values:order?=?(p,d,q)try:mse?=?evaluate_arima_model(dataset,?order)if?mse?<?best_score:best_score,?best_cfg?=?mse,?orderprint('ARIMA%s?MSE=%.3f'?%?(order,mse))except:continueprint('Best?ARIMA%s?MSE=%.3f'?%?(best_cfg,?best_score))參考資料
[1]?
ARIMA 模型:?https://people.duke.edu/~rnau/411arim.htm
[2]?autoRegressive I integrated Moving a average:?https://wikipedia.org/wiki/Autoregressive_integrated_moving_average
[3]?AR自回歸:?https://wikipedia.org/wiki/Autoregressive_integrated_moving_average
[4]?綜合:?https://wikipedia.org/wiki/Order_of_integration
[5]?移動平均:?https://wikipedia.org/wiki/Moving-average_model
[6]?ACF 和 PACF :?https://people.duke.edu/~rnau/411arim3.htm
[7]?'pyramid'庫:?https://alkaline-ml.com/pmdarima/0.9.0/modules/generated/pyramid.arima.auto_arima.html
[8]?wikipedia:?https://wikipedia.org/wiki/Mean_absolute_percentage_error
往期精彩回顧適合初學者入門人工智能的路線及資料下載機器學習及深度學習筆記等資料打印機器學習在線手冊深度學習筆記專輯《統計學習方法》的代碼復現專輯 AI基礎下載黃海廣老師《機器學習課程》視頻課黃海廣老師《機器學習課程》711頁完整版課件本站qq群554839127，加入微信群請掃碼：
總結
以上是生活随笔為你收集整理的理论加实践，终于把时间序列预测ARIMA模型讲明白了的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 如何在RTSP/RTMP直播过程中加入S
- 下一篇： 【机器学习】集成学习代码练习