# 時間序列預測案例研究與 Python：波士頓每月武裝搶劫案 > 原文： [https://machinelearningmastery.com/time-series-forecast-case-study-python-monthly-armed-robberies-boston/](https://machinelearningmastery.com/time-series-forecast-case-study-python-monthly-armed-robberies-boston/) 時間序列預測是一個過程，獲得良好預測的唯一方法是實施此過程。 在本教程中，您將了解如何使用 Python 預測波士頓每月持械搶劫的數量。 完成本教程將為您提供一個框架，用于處理您自己的時間序列預測問題的步驟和工具。 完成本教程后，您將了解： * 如何檢查 Python 環境并仔細定義時間序列預測問題。 * 如何創建測試工具來評估模型，開發基線預測，并使用時間序列分析工具更好地理解您的問題。 * 如何開發自回歸集成移動平均模型，將其保存到文件中，然后加載它以對新時間步驟進行預測。 讓我們開始吧。 * **更新于 2018 年 7 月**：修正了 ACF / PACF 情節編制中的拼寫錯誤（感謝 Patrick Wolf）。  時間序列預測用 Python 進行案例研究 - 波士頓每月武裝搶劫 攝影： [Tim Sackton](https://www.flickr.com/photos/sackton/8301753905/) ，保留一些權利。 ## 概觀 在本教程中，我們將完成從端到端的時間序列預測項目，從下載數據集并定義問題到訓練最終模型和進行預測。 該項目并非詳盡無遺，但通過系統地處理時間序列預測問題，展示了如何快速獲得良好結果。 我們將要完成的這個項目的步驟如下： 1. 環境。 2. 問題描述。 3. 測試線束。 4. 持久性。 5. 數據分析。 6. ARIMA 模型 7. 模型驗證 這將提供一個模板，用于處理您可以在自己的數據集上使用的時間序列預測問題。 ## 1.環境 本教程假定已安裝且正在運行的 SciPy 環境和依賴項，包括： * SciPy 的 * NumPy 的 * Matplotlib * 熊貓 * scikit 學習 * statsmodels 我使用的是 Python 2.7。你是 Python 3 嗎？讓我知道你如何評論。 此腳本將幫助您檢查這些庫的已安裝版本。 ```py # scipy import scipy print('scipy: {}'.format(scipy.__version__)) # numpy import numpy print('numpy: {}'.format(numpy.__version__)) # matplotlib import matplotlib print('matplotlib: {}'.format(matplotlib.__version__)) # pandas import pandas print('pandas: {}'.format(pandas.__version__)) # scikit-learn import sklearn print('sklearn: {}'.format(sklearn.__version__)) # statsmodels import statsmodels print('statsmodels: {}'.format(statsmodels.__version__)) ``` 用于編寫本教程的工作站上的結果如下： ```py scipy: 0.18.1 numpy: 1.11.2 matplotlib: 1.5.3 pandas: 0.19.1 sklearn: 0.18.1 statsmodels: 0.6.1 ``` ## 2.問題描述 問題是要預測美國波士頓每月持械搶劫的數量。 該數據集提供了 1966 年 1 月至 1975 年 10 月在波士頓每月持械搶劫的數量，或者不到 10 年的數據。 這些值是一個計數，有 118 個觀察值。 該數據集歸功于 McCleary＆amp; Hay（1980）。 [您可以了解有關此數據集的更多信息，并直接從 DataMarket](https://datamarket.com/data/set/22ob/monthly-boston-armed-robberies-jan1966-oct1975-deutsch-and-alt-1977) 下載。 將數據集下載為 CSV 文件，并將其放在當前工作目錄中，文件名為“ _robberies.csv_ ”。 ## 3.測試線束 我們必須開發一個測試工具來研究數據并評估候選模型。 這涉及兩個步驟： 1. 定義驗證數據集。 2. 開發模型評估方法。 ### 3.1 驗證數據集 數據集不是最新的。這意味著我們無法輕松收集更新的數據來驗證模型。 因此，我們將假裝它是 1974 年 10 月，并保留最后一年的數據來自分析和模型選擇。 最后一年的數據將用于驗證最終模型。 下面的代碼將數據集作為 Pandas 系列加載并分成兩部分，一部分用于模型開發（ _dataset.csv_ ），另一部分用于驗證（ _validation.csv_ ）。 ```py from pandas import Series series = Series.from_csv('robberies.csv', header=0) split_point = len(series) - 12 dataset, validation = series[0:split_point], series[split_point:] print('Dataset %d, Validation %d' % (len(dataset), len(validation))) dataset.to_csv('dataset.csv') validation.to_csv('validation.csv') ``` 運行該示例將創建兩個文件并打印每個文件中的觀察數。 ```py Dataset 106, Validation 12 ``` 這些文件的具體內容是： * _dataset.csv_ ：1966 年 1 月至 1974 年 10 月的觀察結果（106 次觀察） * _validation.csv_ ：1974 年 11 月至 1975 年 10 月的觀察結果（12 次觀察） 驗證數據集是原始數據集的 10％。 請注意，保存的數據集沒有標題行，因此我們稍后在處理這些文件時無需滿足此要求。 ### 3.2。模型評估 模型評估僅對上一節中準備的 _dataset.csv_ 中的數據進行。 模型評估涉及兩個要素： 1. 表現指標。 2. 測試策略。 #### 3.2.1 績效衡量 這些觀察結果是一系列搶劫案。 我們將使用均方根誤差（RMSE）來評估預測的表現。這將更加重視嚴重錯誤的預測，并且與原始數據具有相同的單位。 在計算和報告 RMSE 之前，必須反轉對數據的任何變換，以使不同方法之間的表現直接相當。 我們可以使用 scikit-learn 庫 _mean_squared_error（）_ 中的輔助函數計算 RMSE，它計算預期值列表（測試集）和預測列表之間的均方誤差。然后我們可以取這個值的平方根來給我們一個 RMSE 分數。 例如： ```py from sklearn.metrics import mean_squared_error from math import sqrt ... test = ... predictions = ... mse = mean_squared_error(test, predictions) rmse = sqrt(mse) print('RMSE: %.3f' % rmse) ``` #### 3.2.2 測試策略 候選模型將使用前向驗證進行評估。 這是因為問題定義需要滾動預測類型模型。這是在給定所有可用數據的情況下需要一步預測的地方。 前瞻性驗證將如下工作： * 數據集的前 50％將被阻止以訓練模型。 * 剩下的 50％的數據集將被迭代并測試模型。 * 對于測試數據集中的每個步驟： * 將訓練模型。 * 進行一步預測并存儲預測以供以后評估。 * 來自測試數據集的實際觀察將被添加到訓練數據集中以用于下一次迭代。 * 將評估在測試數據集的迭代期間做出的預測并報告 RMSE 分數。 鑒于數據的小尺寸，我們將允許在每次預測之前根據所有可用數據重新訓練模型。 我們可以使用簡單的 NumPy 和 Python 代碼編寫測試工具的代碼。 首先，我們可以直接將數據集拆分為訓練集和測試集。我們小心總是將加載的數據集轉換為 float32，以防加載的數據仍然有一些 String 或 Integer 數據類型。 ```py # prepare data X = series.values X = X.astype('float32') train_size = int(len(X) * 0.50) train, test = X[0:train_size], X[train_size:] ``` 接下來，我們可以迭代測試數據集中的時間步長。訓練數據集存儲在 Python 列表中，因為我們需要在每次迭代時輕松附加新的觀察結果，并且 Numpy 數組連接感覺有點矯枉過正。 由于結果或觀察被稱為 _y_ 和 _，_（a' _y []，所以該模型所做的預測被稱為 _yhat_ 。帶有上述標記的 HTG7]是用于預測 _y_ 變量的數學符號。_ 如果模型存在問題，則在每個觀察中打印預測和觀察以進行健全性檢查預測。 ```py # walk-forward validation history = [x for x in train] predictions = list() for i in range(len(test)): # predict yhat = ... predictions.append(yhat) # observation obs = test[i] history.append(obs) print('>Predicted=%.3f, Expected=%3.f' % (yhat, obs)) ``` ## 堅持不懈 在陷入數據分析和建模之前的第一步是建立表現基準。 這將提供用于使用所提出的測試工具評估模型的模板和用于比較所有更精細的預測模型的表現測量。 時間序列預測的基線預測稱為樸素預測或持久性。 這是來自前一時間步驟的觀察被用作下一時間步驟的觀察預測的地方。 我們可以將其直接插入上一節中定義的測試工具中。 完整的代碼清單如下。 ```py from pandas import Series from sklearn.metrics import mean_squared_error from math import sqrt # load data series = Series.from_csv('dataset.csv') # prepare data X = series.values X = X.astype('float32') train_size = int(len(X) * 0.50) train, test = X[0:train_size], X[train_size:] # walk-forward validation history = [x for x in train] predictions = list() for i in range(len(test)): # predict yhat = history[-1] predictions.append(yhat) # observation obs = test[i] history.append(obs) print('>Predicted=%.3f, Expected=%3.f' % (yhat, obs)) # report performance mse = mean_squared_error(test, predictions) rmse = sqrt(mse) print('RMSE: %.3f' % rmse) ``` 運行測試工具會為測試數據集的每次迭代打印預測和觀察。 該示例以打印模型的 RMSE 結束。 在這種情況下，我們可以看到持久性模型實現了 51.844 的 RMSE。這意味著平均而言，每次預測大約有 51 次搶劫，模型錯誤。 ```py ... >Predicted=241.000, Expected=287 >Predicted=287.000, Expected=355 >Predicted=355.000, Expected=460 >Predicted=460.000, Expected=364 >Predicted=364.000, Expected=487 RMSE: 51.844 ``` 我們現在有一個基線預測方法和表現;現在我們可以開始挖掘我們的數據了。 ## 5.數據分析 我們可以使用匯總統計數據和數據圖來快速了解有關預測問題結構的更多信息。 在本節中，我們將從四個角度來看待數據： 1. 摘要統計。 2. 線圖。 3. 密度圖。 4. 盒子和晶須圖。 ### 5.1 摘要統計。 在文本編輯器中打開數據 _dataset.csv_ 文件和/或原始 _robberies.csv_ 文件并查看數據。 快速檢查表明沒有明顯缺失的觀察結果。如果我們試圖將系列強制為浮點值和 _NaN_ 或'_ 等值，我們可能已經注意到了這一點。_ '在數據中。 摘要統計信息可快速查看觀察值的限制。它可以幫助您快速了解我們正在使用的內容。 以下示例計算并打印時間序列的摘要統計信息。 ```py from pandas import Series series = Series.from_csv('dataset.csv') print(series.describe()) ``` 運行該示例提供了許多要查看的摘要統計信息。 這些統計數據的一些觀察包括： * 觀察數量（計數）符合我們的預期，這意味著我們正確處理數據。 * 平均值約為 173，我們可能會考慮這個系列中的水平。 * 在 112 次搶劫中，標準偏差（平均值的平均值差異）相對較大。 * 百分位數和標準偏差確實表明數據的大量傳播。 ```py count 106.000000 mean 173.103774 std 112.231133 min 29.000000 25% 74.750000 50% 144.500000 75% 271.750000 max 487.000000 ``` 如果由隨機波動（例如非系統性）引起，則該系列中的大量傳播可能使得高度準確的預測變得困難。 ### 5.2 線圖 時間序列的線圖可以提供對問題的大量洞察。 下面的示例創建并顯示數據集的線圖。 ```py from pandas import Series from matplotlib import pyplot series = Series.from_csv('dataset.csv') series.plot() pyplot.show() ``` 運行示例并查看繪圖。注意系列中任何明顯的時間結構。 該圖的一些觀察結果包括： * 隨著時間的推移，搶劫的趨勢越來越明顯。 * 似乎沒有任何明顯的異常值。 * 每年都有相對較大的波動，上下波動。 * 晚年的波動似乎比早年的波動大。 * 趨勢意味著數據集幾乎肯定是非平穩的，波動的明顯變化也可能有所貢獻。  每月波士頓搶劫線圖 這些簡單的觀察結果表明，我們可能會看到對趨勢進行建模并從時間序列中刪除它的好處。 或者，我們可以使用差分來使系列靜止以進行建模。如果晚年的波動有增長趨勢，我們甚至可能需要兩個差異水平。 ### 5.3 密度圖 回顧觀察密度圖可以進一步了解數據結構。 下面的示例創建了沒有任何時間結構的觀測的直方圖和密度圖。 ```py from pandas import Series from matplotlib import pyplot series = Series.from_csv('dataset.csv') pyplot.figure(1) pyplot.subplot(211) series.hist() pyplot.subplot(212) series.plot(kind='kde') pyplot.show() ``` 運行示例并查看繪圖。 這些情節的一些觀察包括： * 分布不是高斯分布。 * 分布左移并且可以是指數的或雙高斯的。  每月波士頓搶劫密度情節 我們可能會看到在建模之前探索數據的一些功率變換有一些好處。 ### 5.4 盒子和晶須圖 我們可以按月對月度數據進行分組，并了解每年的觀察結果的傳播情況以及這種情況可能會如何變化。 我們確實希望看到一些趨勢（增加平均值或中位數），但看看分布的其他部分可能會如何變化可能會很有趣。 下面的例子按年份對觀測值進行分組，并為每年的觀測值創建一個方框和胡須圖。去年（1974 年）僅包含 10 個月，與其他年份的其他 12 個月的觀察結果可能無差別。因此，僅繪制了 1966 年至 1973 年之間的數據。 ```py from pandas import Series from pandas import DataFrame from pandas import TimeGrouper from matplotlib import pyplot series = Series.from_csv('dataset.csv') groups = series['1966':'1973'].groupby(TimeGrouper('A')) years = DataFrame() for name, group in groups: years[name.year] = group.values years.boxplot() pyplot.show() ``` 運行該示例并排創建 8 個框和胡須圖，每 8 年選定數據一個。 審查該情節的一些觀察包括： * 每年的中值（紅線）顯示可能不是線性的趨勢。 * 傳播或中間 50％的數據（藍框）有所不同，但隨著時間的推移可能不一致。 * 前幾年，也許是前兩年，與其他數據集完全不同。  每月波士頓搶劫箱和晶須地塊 觀察結果表明，年度波動可能不是系統性的，難以建模。他們還建議，如果確實存在很大差異，那么從建模中剪下前兩年的數據可能會有一些好處。 這種年度數據視圖是一個有趣的途徑，可以通過查看逐年的摘要統計數據和每年的摘要統計數據的變化來進一步追求。 接下來，我們可以開始研究該系列的預測模型。 ## 6\. ARIMA 模型 在本節中，我們將針對該問題開發自動回歸集成移動平均線（ARIMA）模型。 我們將分四步處理： 1. 開發手動配置的 ARIMA 模型。 2. 使用 ARIMA 的網格搜索來查找優化模型。 3. 分析預測殘差以評估模型中的任何偏差。 4. 使用電源轉換探索對模型的改進。 ### 6.1 手動配置 ARIMA 非季節性 ARIMA（p，d，q）需要三個參數，并且傳統上是手動配置的。 對時間序列數據的分析假設我們正在使用固定的時間序列。 時間序列幾乎肯定是非平穩的。我們可以通過首先對系列進行差分并使用統計檢驗確認結果是靜止的來使其靜止。 下面的示例創建了該系列的固定版本并將其保存到文件 _stationary.csv_ 。 ```py from pandas import Series from statsmodels.tsa.stattools import adfuller # create a differe def difference(dataset): diff = list() for i in range(1, len(dataset)): value = dataset[i] - dataset[i - 1] diff.append(value) return Series(diff) series = Series.from_csv('dataset.csv') X = series.values # difference data stationary = difference(X) stationary.index = series.index[1:] # check if stationary result = adfuller(stationary) print('ADF Statistic: %f' % result[0]) print('p-value: %f' % result[1]) print('Critical Values:') for key, value in result[4].items(): print('\t%s: %.3f' % (key, value)) # save stationary.to_csv('stationary.csv') ``` 運行該示例輸出系列是否靜止的統計顯著性檢驗的結果。具體來說，增強 Dickey-Fuller 測試。 結果表明，檢驗統計值-3.980946 小于-2.893 的 5％的臨界值。這表明我們可以拒絕具有小于 5％的顯著性水平的零假設（即，結果是統計吸蟲的概率很低）。 拒絕原假設意味著該過程沒有單位根，反過來，時間序列是靜止的或沒有時間依賴的結構。 ```py ADF Statistic: -3.980946 p-value: 0.001514 Critical Values: 5%: -2.893 1%: -3.503 10%: -2.584 ``` 這表明至少需要一個差分水平。我們 ARIMA 模型中的 _d_ 參數至少應為 1。 下一步是分別為自回歸（AR）和移動平均（MA）參數 p 和 q 選擇滯后值。 我們可以通過查看自相關函數（ACF）和部分自相關函數（PACF）圖來做到這一點。 下面的示例為該系列創建了 ACF 和 PACF 圖。 ```py from pandas import Series from statsmodels.graphics.tsaplots import plot_acf from statsmodels.graphics.tsaplots import plot_pacf from matplotlib import pyplot series = Series.from_csv('stationary.csv') pyplot.figure() pyplot.subplot(211) plot_acf(series, ax=pyplot.gca()) pyplot.subplot(212) plot_pacf(series, ax=pyplot.gca()) pyplot.show() ``` 運行該示例并查看繪圖，以獲得有關如何為 ARIMA 模型設置 _p_ 和 _q_ 變量的見解。 以下是該圖的一些觀察結果。 * ACF 顯示 1 個月的顯著滯后。 * PACF 顯示出可能持續 2 個月的顯著滯后，顯著滯后可能持續 12 個月。 * ACF 和 PACF 都在同一點下降，可能表明 AR 和 MA 的混合。 _p_ 和 _q_ 值的良好起點是 1 或 2。  每月波士頓搶劫 ACF 和 PACF 情節 這種快速分析表明原始數據上的 ARIMA（1,1,2）可能是一個很好的起點。 實驗表明，ARIMA 的這種配置不會收斂并導致底層庫出錯。一些實驗表明，該模型看起來并不穩定，同時定義了非零 AR 和 MA 訂單。 該模型可以簡化為 ARIMA（0,1,2）。下面的示例演示了此 ARIMA 模型在測試工具上的表現。 ```py from pandas import Series from sklearn.metrics import mean_squared_error from statsmodels.tsa.arima_model import ARIMA from math import sqrt # load data series = Series.from_csv('dataset.csv') # prepare data X = series.values X = X.astype('float32') train_size = int(len(X) * 0.50) train, test = X[0:train_size], X[train_size:] # walk-forward validation history = [x for x in train] predictions = list() for i in range(len(test)): # predict model = ARIMA(history, order=(0,1,2)) model_fit = model.fit(disp=0) yhat = model_fit.forecast()[0] predictions.append(yhat) # observation obs = test[i] history.append(obs) print('>Predicted=%.3f, Expected=%3.f' % (yhat, obs)) # report performance mse = mean_squared_error(test, predictions) rmse = sqrt(mse) print('RMSE: %.3f' % rmse) ``` 運行此示例會導致 RMSE 為 49.821，低于持久性模型。 ```py ... >Predicted=280.614, Expected=287 >Predicted=302.079, Expected=355 >Predicted=340.210, Expected=460 >Predicted=405.172, Expected=364 >Predicted=333.755, Expected=487 RMSE: 49.821 ``` 這是一個良好的開端，但我們可以通過更好的配置 ARIMA 模型獲得改進的結果。 ### 6.2 網格搜索 ARIMA 超參數 許多 ARIMA 配置在此數據集上不穩定，但可能有其他超參數導致表現良好的模型。 在本節中，我們將搜索 _p_ ， _d_ 和 _q_ 的值，以查找不會導致錯誤的組合，并找到導致最棒的表演。我們將使用網格搜索來探索整數值子集中的所有組合。 具體來說，我們將搜索以下參數的所有組合： * _p_ ：0 到 12。 * _d_ ：0 到 3。 * _q_ ：0 到 12。 這是（ _13 * 4 * 13）或 676_ 運行的測試工具，需要一些時間才能執行。 下面列出了測試工具的網格搜索版本的完整工作示例。 ```py import warnings from pandas import Series from statsmodels.tsa.arima_model import ARIMA from sklearn.metrics import mean_squared_error from math import sqrt # evaluate an ARIMA model for a given order (p,d,q) and return RMSE def evaluate_arima_model(X, arima_order): # prepare training dataset X = X.astype('float32') train_size = int(len(X) * 0.50) train, test = X[0:train_size], X[train_size:] history = [x for x in train] # make predictions predictions = 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 error mse = mean_squared_error(test, predictions) rmse = sqrt(mse) return rmse # 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"), None for 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, order print('ARIMA%s MSE=%.3f' % (order,mse)) except: continue print('Best ARIMA%s MSE=%.3f' % (best_cfg, best_score)) # load dataset series = Series.from_csv('dataset.csv') # evaluate parameters p_values = range(0,13) d_values = range(0, 4) q_values = range(0, 13) warnings.filterwarnings("ignore") evaluate_models(series.values, p_values, d_values, q_values) ``` 運行該示例將遍歷所有組合，并在收斂且無錯誤的情況下報告結果。 該示例在現代硬件上運行需要不到 2 小時。 結果表明，發現的最佳配置是 ARIMA（0,1,2）;巧合的是，這在前一節中得到了證明。 ```py ... ARIMA(6, 1, 0) MSE=52.437 ARIMA(6, 2, 0) MSE=58.307 ARIMA(7, 1, 0) MSE=51.104 ARIMA(7, 1, 1) MSE=52.340 ARIMA(8, 1, 0) MSE=51.759 Best ARIMA(0, 1, 2) MSE=49.821 ``` ### 6.3 查看殘留錯誤 對模型進行良好的最終檢查是檢查殘差預測誤差。 理想情況下，殘差的分布應該是具有零均值的高斯分布。 我們可以通過用直方圖和密度圖繪制殘差來檢查這一點。 下面的示例計算測試集上預測的殘差，并創建這些密度圖。 ```py from pandas import Series from pandas import DataFrame from sklearn.metrics import mean_squared_error from statsmodels.tsa.arima_model import ARIMA from math import sqrt from matplotlib import pyplot # load data series = Series.from_csv('dataset.csv') # prepare data X = series.values X = X.astype('float32') train_size = int(len(X) * 0.50) train, test = X[0:train_size], X[train_size:] # walk-foward validation history = [x for x in train] predictions = list() for i in range(len(test)): # predict model = ARIMA(history, order=(0,1,2)) model_fit = model.fit(disp=0) yhat = model_fit.forecast()[0] predictions.append(yhat) # observation obs = test[i] history.append(obs) # errors residuals = [test[i]-predictions[i] for i in range(len(test))] residuals = DataFrame(residuals) pyplot.figure() pyplot.subplot(211) residuals.hist(ax=pyplot.gca()) pyplot.subplot(212) residuals.plot(kind='kde', ax=pyplot.gca()) pyplot.show() ``` 運行該示例將創建兩個圖。 圖表顯示具有較長右尾的高斯分布。 這可能是預測偏差的一個標志，在這種情況下，在建模之前可能對原始數據進行基于功率的轉換可能是有用的。  剩余誤差密度圖 檢查任何類型的自相關的殘差的時間序列也是一個好主意。如果存在，則表明該模型有更多機會對數據中的時間結構進行建模。 下面的示例重新計算殘差，并創建 ACF 和 PACF 圖以檢查是否存在任何顯著的自相關。 ```py from pandas import Series from pandas import DataFrame from sklearn.metrics import mean_squared_error from statsmodels.tsa.arima_model import ARIMA from statsmodels.graphics.tsaplots import plot_acf from statsmodels.graphics.tsaplots import plot_pacf from math import sqrt from matplotlib import pyplot # load data series = Series.from_csv('dataset.csv') # prepare data X = series.values X = X.astype('float32') train_size = int(len(X) * 0.50) train, test = X[0:train_size], X[train_size:] # walk-foward validation history = [x for x in train] predictions = list() for i in range(len(test)): # predict model = ARIMA(history, order=(0,1,2)) model_fit = model.fit(disp=0) yhat = model_fit.forecast()[0] predictions.append(yhat) # observation obs = test[i] history.append(obs) # errors residuals = [test[i]-predictions[i] for i in range(len(test))] residuals = DataFrame(residuals) pyplot.figure() pyplot.subplot(211) plot_acf(residuals, ax=pyplot.gca()) pyplot.subplot(212) plot_pacf(residuals, ax=pyplot.gca()) pyplot.show() ``` 結果表明，該模型捕獲了時間序列中存在的小自相關。  剩余錯誤 ACF 和 PACF 圖 ### 6.4 Box-Cox 轉換數據集 Box-Cox 變換是一種能夠評估一組功率變換的方法，包括但不限于數據的對數，平方根和互易變換。 下面的示例執行數據的對數轉換，并生成一些圖以查看對時間序列的影響。 ```py from pandas import Series from pandas import DataFrame from scipy.stats import boxcox from matplotlib import pyplot from statsmodels.graphics.gofplots import qqplot series = Series.from_csv('dataset.csv') X = series.values transformed, lam = boxcox(X) print('Lambda: %f' % lam) pyplot.figure(1) # line plot pyplot.subplot(311) pyplot.plot(transformed) # histogram pyplot.subplot(312) pyplot.hist(transformed) # q-q plot pyplot.subplot(313) qqplot(transformed, line='r', ax=pyplot.gca()) pyplot.show() ``` 運行該示例將創建三個圖形：轉換時間序列的折線圖，顯示轉換值分布的直方圖，以及顯示值與理想化高斯分布相比如何分布的 Q-Q 圖。 這些圖的一些觀察結果如下： * 從時間序列的線圖中刪除了大的波動。 * 直方圖顯示更平坦或更均勻（表現良好）的值分布。 * Q-Q 圖是合理的，但仍然不適合高斯分布。  Box Cox 變換每月波士頓搶劫案的情節 毫無疑問，Box-Cox 變換已經對時間序列做了些什么，可能會有用。 在繼續使用轉換后的數據測試 ARIMA 模型之前，我們必須有一種方法來反轉變換，以便將使用經過轉換的數據訓練的模型所做的預測轉換回原始比例。 示例中使用的 _boxcox（）_ 函數通過優化成本函數找到理想的 lambda 值。 lambda 用于以下函數中以轉換數據： ```py transform = log(x), if lambda == 0 transform = (x^lambda - 1) / lambda, if lambda != 0 ``` 此轉換函數可以直接反轉，如下所示： ```py x = exp(transform) if lambda == 0 x = exp(log(lambda * transform + 1) / lambda) ``` 這個逆 Box-Cox 變換函數可以在 Python 中實現，如下所示： ```py # invert box-cox transform from math import log from math import exp def boxcox_inverse(value, lam): if lam == 0: return exp(value) return exp(log(lam * value + 1) / lam) ``` 我們可以使用 Box-Cox 變換重新評估 ARIMA（0,1,2）模型。 這包括在擬合 ARIMA 模型之前首先變換歷史，然后在存儲之前反轉預測變換，以便稍后與預期值進行比較。 _boxcox（）_ 功能可能會失敗。在實踐中，我已經看到了這一點，并且似乎通過返回的 _lambda_ 值小于-5 來表示。按照慣例，_λ_ 值在-5 和 5 之間評估。 對于小于-5 的 lambda 值添加檢查，如果是這種情況，則假定 _lambda_ 值為 1，并且原始歷史記錄用于擬合模型。 _λ_ 值為 1 與“_ 無變換 _”相同，因此逆變換無效。 下面列出了完整的示例。 ```py from pandas import Series from sklearn.metrics import mean_squared_error from statsmodels.tsa.arima_model import ARIMA from math import sqrt from math import log from math import exp from scipy.stats import boxcox # invert box-cox transform def boxcox_inverse(value, lam): if lam == 0: return exp(value) return exp(log(lam * value + 1) / lam) # load data series = Series.from_csv('dataset.csv') # prepare data X = series.values X = X.astype('float32') train_size = int(len(X) * 0.50) train, test = X[0:train_size], X[train_size:] # walk-foward validation history = [x for x in train] predictions = list() for i in range(len(test)): # transform transformed, lam = boxcox(history) if lam < -5: transformed, lam = history, 1 # predict model = ARIMA(transformed, order=(0,1,2)) model_fit = model.fit(disp=0) yhat = model_fit.forecast()[0] # invert transformed prediction yhat = boxcox_inverse(yhat, lam) predictions.append(yhat) # observation obs = test[i] history.append(obs) print('>Predicted=%.3f, Expected=%3.f' % (yhat, obs)) # report performance mse = mean_squared_error(test, predictions) rmse = sqrt(mse) print('RMSE: %.3f' % rmse) ``` 運行該示例會在每次迭代時打印預測值和期望值。 注意，使用 box cox（）轉換函數時可能會看到警告;例如： ```py RuntimeWarning: overflow encountered in square llf -= N / 2.0 * np.log(np.sum((y - y_mean)**2\. / N, axis=0)) ``` 這些可以暫時忽略。 轉換數據模型的最終 RMSE 為 49.443。對于未轉換的數據，這是一個比 ARIMA 模型更小的錯誤，但只是略有不同，它可能有統計差異，也可能沒有統計差異。 ```py ... >Predicted=276.253, Expected=287 >Predicted=299.811, Expected=355 >Predicted=338.997, Expected=460 >Predicted=404.509, Expected=364 >Predicted=349.336, Expected=487 RMSE: 49.443 ``` 我們將使用此模型與 Box-Cox 變換作為最終模型。 ## 7.模型驗證 在開發模型并選擇最終模型后，必須對其進行驗證和最終確定。 驗證是流程的一個可選部分，但提供“最后檢查”以確保我們沒有欺騙或欺騙自己。 本節包括以下步驟： 1. **Finalize Model** ：訓練并保存最終模型。 2. **進行預測**：加載最終模型并進行預測。 3. **驗證模型**：加載并驗證最終模型。 ### 7.1 完成模型 最終確定模型涉及在整個數據集上擬合 ARIMA 模型，在這種情況下，在整個數據集的變換版本上。 一旦適合，模型可以保存到文件中供以后使用。因為 Box-Cox 變換也是對數據執行的，所以我們需要知道所選擇的 lambda，以便模型中的任何預測都可以轉換回原始的，未變換的比例。 下面的示例適合 Box-Cox 變換數據集上的 ARIMA（0,1,2）模型，并將整個擬合對象和 lambda 值保存到文件中。 ```py from pandas import Series from statsmodels.tsa.arima_model import ARIMA from scipy.stats import boxcox import numpy # load data series = Series.from_csv('dataset.csv') # prepare data X = series.values X = X.astype('float32') # transform data transformed, lam = boxcox(X) # fit model model = ARIMA(transformed, order=(0,1,2)) model_fit = model.fit(disp=0) # save model model_fit.save('model.pkl') numpy.save('model_lambda.npy', [lam]) ``` 運行該示例將創建兩個本地文件： * **model.pkl** 這是調用 _ARIMA.fit（）_ 的 ARIMAResult 對象。這包括系數和擬合模型時返回的所有其他內部數據。 * **model_lambda.npy** 這是存儲為單行，單列 NumPy 數組的 lambda 值。 這可能是矯枉過正的，操作使用真正需要的是模型中的 AR 和 MA 系數，差異數量的 _d_ 參數，可能是滯后觀察和模型殘差，以及λ值為了變革。 ### 7.2 進行預測 一個自然的例子可能是加載模型并進行單一預測。 這是相對簡單的，包括恢復保存的模型和 lambda 并調用 _forecast（）_ 方法。 嗯，它應該是，但是當前穩定版本的 statsmodels 庫（v0.6.1）中存在一個錯誤，它失敗并出現錯誤： ```py TypeError: __new__() takes at least 3 arguments (1 given) ``` 當我測試它時，這個 bug 似乎也出現在 statsmodels 的 0.8 版本候選版本 1 中。 有關詳細信息，請參閱 [Zae Myung Kim](http://zaemyung.com/) 的[討論并在 GitHub 上解決此問題](https://github.com/statsmodels/statsmodels/pull/3217)。 我們可以使用一個猴子補丁來解決這個問題，該補丁在保存之前將 ___getnewargs __（）_ 實例函數添加到 ARIMA 類。 以下示例與上一部分相同，只有這一小改動。 ```py from pandas import Series from statsmodels.tsa.arima_model import ARIMA from scipy.stats import boxcox import numpy # monkey patch around bug in ARIMA class def __getnewargs__(self): return ((self.endog),(self.k_lags, self.k_diff, self.k_ma)) ARIMA.__getnewargs__ = __getnewargs__ # load data series = Series.from_csv('dataset.csv') # prepare data X = series.values X = X.astype('float32') # transform data transformed, lam = boxcox(X) # fit model model = ARIMA(transformed, order=(0,1,2)) model_fit = model.fit(disp=0) # save model model_fit.save('model.pkl') numpy.save('model_lambda.npy', [lam]) ``` 我們現在可以加載模型并進行單一預測。 下面的示例加載模型，對下一個時間步進行預測，反轉 Box-Cox 變換，并打印預測。 ```py from pandas import Series from statsmodels.tsa.arima_model import ARIMAResults from math import exp from math import log import numpy # invert box-cox transform def boxcox_inverse(value, lam): if lam == 0: return exp(value) return exp(log(lam * value + 1) / lam) model_fit = ARIMAResults.load('model.pkl') lam = numpy.load('model_lambda.npy') yhat = model_fit.forecast()[0] yhat = boxcox_inverse(yhat, lam) print('Predicted: %.3f' % yhat) ``` 運行該示例打印約 452 的預測。 如果我們查看 validation.csv，我們可以看到下一個時間段第一行的值是 452.模型得到它 100％正確，這是非常令人印象深刻的（或幸運）。 ```py Predicted: 452.043 ``` ### 7.3 驗證模型 我們可以加載模型并以假裝操作方式使用它。 在測試工具部分中，我們將原始數據集的最后 12 個月保存在單獨的文件中以驗證最終模型。 我們現在可以加載這個 _validation.csv_ 文件并使用它看看我們的模型在“看不見的”數據上的真實程度。 我們可以通過兩種方式進行： * 加載模型并使用它來預測未來 12 個月。超過前一兩個月的預測將很快開始降低技能。 * 加載模型并以滾動預測方式使用它，更新每個時間步的變換和模型。這是首選方法，因為它將如何在實踐中使用此模型，因為它將實現最佳表現。 與前幾節中的模型評估一樣，我們將以滾動預測的方式進行預測。這意味著我們將在驗證數據集中逐步超過提前期，并將觀察結果作為歷史記錄的更新。 ```py from pandas import Series from matplotlib import pyplot from statsmodels.tsa.arima_model import ARIMA from statsmodels.tsa.arima_model import ARIMAResults from scipy.stats import boxcox from sklearn.metrics import mean_squared_error from math import sqrt from math import exp from math import log import numpy # invert box-cox transform def boxcox_inverse(value, lam): if lam == 0: return exp(value) return exp(log(lam * value + 1) / lam) # load and prepare datasets dataset = Series.from_csv('dataset.csv') X = dataset.values.astype('float32') history = [x for x in X] validation = Series.from_csv('validation.csv') y = validation.values.astype('float32') # load model model_fit = ARIMAResults.load('model.pkl') lam = numpy.load('model_lambda.npy') # make first prediction predictions = list() yhat = model_fit.forecast()[0] yhat = boxcox_inverse(yhat, lam) predictions.append(yhat) history.append(y[0]) print('>Predicted=%.3f, Expected=%3.f' % (yhat, y[0])) # rolling forecasts for i in range(1, len(y)): # transform transformed, lam = boxcox(history) if lam < -5: transformed, lam = history, 1 # predict model = ARIMA(transformed, order=(0,1,2)) model_fit = model.fit(disp=0) yhat = model_fit.forecast()[0] # invert transformed prediction yhat = boxcox_inverse(yhat, lam) predictions.append(yhat) # observation obs = y[i] history.append(obs) print('>Predicted=%.3f, Expected=%3.f' % (yhat, obs)) # report performance mse = mean_squared_error(y, predictions) rmse = sqrt(mse) print('RMSE: %.3f' % rmse) pyplot.plot(y) pyplot.plot(predictions, color='red') pyplot.show() ``` 運行該示例將打印驗證數據集中時間步長的每個預測值和預期值。 驗證期的最終 RMSE 預計為 53 次搶劫。這與 49 的預期誤差沒有太大的不同，但我希望它與簡單的持久性模型也沒有太大區別。 ```py >Predicted=452.043, Expected=452 >Predicted=423.088, Expected=391 >Predicted=408.378, Expected=500 >Predicted=482.454, Expected=451 >Predicted=445.944, Expected=375 >Predicted=413.881, Expected=372 >Predicted=413.209, Expected=302 >Predicted=355.159, Expected=316 >Predicted=363.515, Expected=398 >Predicted=406.365, Expected=394 >Predicted=394.186, Expected=431 >Predicted=428.174, Expected=431 RMSE: 53.078 ``` 還提供了與驗證數據集相比較的預測圖。 預測確實具有持久性預測的特征。這確實表明雖然這個時間序列確實有明顯的趨勢，但它仍然是一個相當困難的問題。  驗證時間步驟的預測 ## 教程擴展 本教程并非詳盡無遺;你可以做更多的事情來改善結果。 本節列出了一些想法。 * **統計顯著性檢驗**。使用統計檢驗來檢查不同模型之間的結果差異是否具有統計學意義。學生 t 檢驗將是一個很好的起點。 * **使用數據轉換進行網格搜索**。使用 Box-Cox 變換在 ARIMA 超參數中重復網格搜索，并查看是否可以實現不同且更好的參數集。 * **檢查殘留物**。使用 Box-Cox 變換調查最終模型上的殘差預測誤差，以查看是否存在可以解決的進一步偏差和/或自相關。 * **精益模型保存**。簡化模型保存僅存儲所需的系數而不是整個 ARIMAResults 對象。 * **手動處理趨勢**。使用線性或非線性模型直接對趨勢建模，并將其明確地從系列中刪除。如果趨勢是非線性的并且可以比線性情況更好地建模，則這可以導致更好的表現。 * **置信區間**。顯示驗證數據集上預測的置信區間。 * **數據選擇**。考慮在沒有前兩年數據的情況下對問題進行建模，看看這是否會對預測技能產生影響。 您是否實施了這些擴展？你能獲得更好的結果嗎？ 在下面的評論中分享您的發現。 ## 摘要 在本教程中，您使用 Python 發現了時間序列預測項目的步驟和工具。 我們在本教程中介紹了很多內容，具體來說： * 如何開發具有表現測量和評估方法的測試工具，以及如何快速開發基線預測和技能。 * 如何使用時間序列分析來提出如何最好地模擬預測問題的想法。 * 如何開發 ARIMA 模型，保存它，然后加載它以對新數據進行預測。 你是怎么做的？您對本教程有任何疑問嗎？ 在下面的評論中提出您的問題，我會盡力回答。