# Python 的時間序列預測案例研究：巴爾的摩的年度用水量 > 原文： [https://machinelearningmastery.com/time-series-forecast-study-python-annual-water-usage-baltimore/](https://machinelearningmastery.com/time-series-forecast-study-python-annual-water-usage-baltimore/) 時間序列預測是一個過程，獲得良好預測的唯一方法是實施此過程。 在本教程中，您將了解如何使用 Python 預測巴爾的摩的年用水量。 完成本教程將為您提供一個框架，用于處理您自己的時間序列預測問題的步驟和工具。 完成本教程后，您將了解： * 如何確認您的 Python 環境并仔細定義時間序列預測問題。 * 如何創建測試工具來評估模型，開發基線預測，并使用時間序列分析工具更好地理解您的問題。 * 如何開發自回歸集成移動平均模型，將其保存到文件中，然后加載它以對新時間步驟進行預測。 讓我們開始吧。  使用 Python 進行時間序列預測研究 - 巴爾的摩年度用水量 照片由 [Andy Mitchell](https://www.flickr.com/photos/monstermunch/3085401491/) 拍攝，保留一些權利。 ## 概觀 在本教程中，我們將完成從端到端的時間序列預測項目，從下載數據集并定義問題到訓練最終模型和進行預測。 該項目并非詳盡無遺，但通過系統地處理時間序列預測問題，展示了如何快速獲得良好結果。 我們將要完成的這個項目的步驟如下。 1. 環境。 2. 問題描述。 3. 測試線束。 4. 持久性。 5. 數據分析。 6. ARIMA 模型。 7. 模型驗證。 這將提供一個模板，用于處理您可以在自己的數據集上使用的時間序列預測問題。 ## 1.環境 本教程假定已安裝且正在運行的 SciPy 環境和依賴項，包括： * SciPy 的 * NumPy 的 * Matplotlib * 熊貓 * scikit 學習 * statsmodels 如果您需要在工作站上安裝 Python 和 SciPy 環境的幫助，請考慮為您管理大部分內容的 [Anaconda 發行版](https://www.continuum.io/downloads)。 此腳本將幫助您檢查這些庫的已安裝版本。 ```py # scipy import scipy print('scipy: %s' % scipy.__version__) # numpy import numpy print('numpy: %s' % numpy.__version__) # matplotlib import matplotlib print('matplotlib: %s' % matplotlib.__version__) # pandas import pandas print('pandas: %s' % pandas.__version__) # scikit-learn import sklearn print('sklearn: %s' % sklearn.__version__) # statsmodels import statsmodels print('statsmodels: %s' % 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.問題描述 問題是預測年用水量。 該數據集提供了 1885 年至 1963 年巴爾的摩的年用水量，或 79 年的數據。 這些值以每人每天的升數為單位，有 79 個觀測值。 該數據集歸功于 Hipel 和 McLeod，1994。 [您可以了解有關此數據集的更多信息，并直接從 DataMarket](https://datamarket.com/data/set/22sl/baltmore-city-annual-water-use-liters-per-capita-per-day-1885-1968) 下載。 將數據集下載為 CSV 文件，并將其放在當前工作目錄中，文件名為“ _water.csv_ ”。 ## 3.測試線束 我們必須開發一個測試工具來研究數據并評估候選模型。 這涉及兩個步驟： 1. 定義驗證數據集。 2. 開發模型評估方法。 ### 3.1 驗證數據集 數據集不是最新的。這意味著我們無法輕松收集更新的數據來驗證模型。 因此，我們將假裝它是 1953 年并且從分析和模型選擇中扣留最近 10 年的數據。 最后十年的數據將用于驗證最終模型。 下面的代碼將數據集作為 Pandas 系列加載并分成兩部分，一部分用于模型開發（ _dataset.csv_ ），另一部分用于驗證（ _validation.csv_ ）。 ```py from pandas import Series series = Series.from_csv('water.csv', header=0) split_point = len(series) - 10 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 69, Validation 10 ``` 這些文件的具體內容是： * _dataset.csv_ ：1885 年至 1953 年的觀察結果（69 次觀察）。 * _validation.csv_ ：1954 年至 1963 年的觀察結果（10 次觀察）。 驗證數據集約占原始數據集的 12％。 請注意，保存的數據集沒有標題行，因此我們稍后在處理這些文件時無需滿足此要求。 ### 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 代碼編寫測試工具的代碼。 首先，我們可以直接將數據集拆分為訓練集和測試集。如果加載的數據仍然有一些 _String_ 或 _Integer_ 數據類型，我們小心地始終將加載的數據集轉換為 _float32_ 。 ```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 結束。 在這種情況下，我們可以看到持久性模型實現了 21.975 的 RMSE。這意味著平均而言，每個預測的模型錯誤人均每天約 22 升。 ```py ... >Predicted=613.000, Expected=598 >Predicted=598.000, Expected=575 >Predicted=575.000, Expected=564 >Predicted=564.000, Expected=549 >Predicted=549.000, Expected=538 RMSE: 21.975 ``` 我們現在有一個基線預測方法和表現;現在我們可以開始挖掘我們的數據了。 ## 5.數據分析 我們可以使用匯總統計數據和數據圖來快速了解有關預測問題結構的更多信息。 在本節中，我們將從四個角度來看待數據： 1. 摘要統計。 2. 線圖。 3. 密度圖。 4. 盒子和晶須圖。 ### 5.1。摘要統計 摘要統計信息可快速查看觀察值的限制。它可以幫助您快速了解我們正在使用的內容。 以下示例計算并打印時間序列的摘要統計信息。 ```py from pandas import Series series = Series.from_csv('dataset.csv') print(series.describe()) ``` 運行該示例提供了許多要查看的摘要統計信息。 這些統計數據的一些觀察包括： * 觀察數量（計數）符合我們的預期，這意味著我們正確處理數據。 * 平均值約為 500，我們可能會考慮這個系列中的水平。 * 標準偏差和百分位數表明平均值的合理緊密擴散。 ```py count 69.000000 mean 500.478261 std 73.901685 min 344.000000 25% 458.000000 50% 492.000000 75% 538.000000 max 662.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。盒子和晶須圖 我們可以將年度數據分組十年，并了解每十年觀測的傳播以及這可能會如何變化。 我們確實希望看到一些趨勢（增加平均值或中位數），但看看分布的其他部分可能會如何變化可能會很有趣。 下面的例子將觀察結果按十年分組，并為每個十年的觀測創建一個盒子和須狀圖。過去十年只有 9 年，可能與其他幾十年沒有用。因此，僅繪制了 1885 年至 1944 年之間的數據。 ```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['1885':'1944'].groupby(TimeGrouper('10AS')) decades = DataFrame() for name, group in groups: decades[name.year] = group.values decades.boxplot() pyplot.show() ``` 運行該示例并排創建 6 個框和胡須圖，一個用于 60 年的選定數據。 審查該情節的一些觀察包括： * 每年的中值（紅線）可能顯示可能不是線性的增加趨勢。 * 傳播或中間 50％的數據（藍色框）確實顯示出一些變化。 * 在幾十年內可能存在異常值（在盒子外面交叉和胡須）。 * 第二個到過去十年似乎平均消費較低，可能與第一次世界大戰有關。  年度用水箱和晶須地塊 這種年度數據視圖是一個有趣的途徑，可以通過查看十年到十年的匯總統計數據和匯總統計數據的變化來進一步追求。 ## 6\. ARIMA 模型 在本節中，我們將針對該問題開發自回歸集成移動平均線或 ARIMA 模型。 我們將通過手動和自動配置 ARIMA 模型來進行建模。接下來是調查所選模型的殘差的第三步。 因此，本節分為 3 個步驟： 1. 手動配置 ARIMA。 2. 自動配置 ARIMA。 3. 查看殘留錯誤。 ### 6.1 手動配置 ARIMA ARIMA（ _p，d，q_ ）模型需要三個參數，并且傳統上是手動配置的。 對時間序列數據的分析假設我們正在使用固定的時間序列。 時間序列可能是非平穩的。我們可以通過首先對系列進行差分并使用統計檢驗確認結果是靜止的來使其靜止。 下面的示例創建了該系列的固定版本并將其保存到文件 _stationary.csv_ 。 ```py from pandas import Series from statsmodels.tsa.stattools import adfuller from matplotlib import pyplot # 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 X = X.astype('float32') # 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)) # plot differenced data stationary.plot() pyplot.show() # save stationary.to_csv('stationary.csv') ``` 運行該示例輸出差異系列是否靜止的統計顯著性檢驗的結果。具體來說，增強 Dickey-Fuller 測試。 結果表明，檢驗統計值-6.126719 小于-3.534 的 1％的臨界值。這表明我們可以拒絕具有小于 1％的顯著性水平的零假設（即，結果是統計僥幸的低概率）。 拒絕原假設意味著該過程沒有單位根，反過來，時間序列是靜止的或沒有時間依賴的結構。 ```py ADF Statistic: -6.126719 p-value: 0.000000 Critical Values: 5%: -2.906 1%: -3.534 10%: -2.591 ``` 這表明至少需要一個差分水平。我們 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('dataset.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 沒有顯著滯后。 * PACF 也沒有顯著滯后。 _p_ 和 _q_ 值的良好起點也是 0。  固定年用水量數據集的 ACF 和 PACF 圖 這種快速分析表明原始數據上的 ARIMA（0,1,0）可能是一個很好的起點。 這實際上是一種持久性模型。下面列出了完整的示例。 ```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,0)) 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 為 22.311，略高于上面的持久性模型。 這可能是因為 ARIMA 實施的細節，例如計算和添加的自動趨勢常數。 ```py ... >Predicted=617.079, Expected=598 >Predicted=601.781, Expected=575 >Predicted=578.369, Expected=564 >Predicted=567.152, Expected=549 >Predicted=551.881, Expected=538 RMSE: 22.311 ``` ### 6.2 網格搜索 ARIMA 超參數 ACF 和 PACF 圖顯示我們不能比此數據集上的持久性模型做得更好。 為了確認這一分析，我們可以對一套 ARIMA 超參數進行網格搜索，并檢查沒有模型可以獲得更好的樣本外 RMSE 表現。 在本節中，我們將搜索組合的 _p_ ， _d_ 和 _q_ 的值（跳過那些未收斂的組合），并找到結果的組合在最好的表現。我們將使用網格搜索來探索整數值子集中的所有組合。 具體來說，我們將搜索以下參數的所有組合： * _p_ ：0 到 4。 * _d_ ：0 到 2。 * _q_ ：0 到 4。 這是（ _5 * 3 * 5_ ），或測試線束的 300 次潛在運行，并且將花費一些時間來執行。 我們還將禁用從模型中自動添加趨勢常數，方法是將'_ 趨勢 _'參數設置為' _nc_ '，當調用 [fit（）[HTG5 時]為無常數]。](http://statsmodels.sourceforge.net/stable/generated/statsmodels.tsa.arima_model.ARIMA.fit.html) 下面列出了測試工具的網格搜索版本的完整工作示例。 ```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) model_fit = model.fit(trend='nc', 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 RMSE=%.3f' % (order,mse)) except: continue print('Best ARIMA%s RMSE=%.3f' % (best_cfg, best_score)) # load dataset series = Series.from_csv('dataset.csv') # evaluate parameters p_values = range(0, 5) d_values = range(0, 3) q_values = range(0, 5) warnings.filterwarnings("ignore") evaluate_models(series.values, p_values, d_values, q_values) ``` 運行該示例將遍歷所有組合，并在收斂且無錯誤的情況下報告結果。該示例需要 2 分多鐘才能在現代硬件上運行。 結果顯示，發現的最佳配置是 ARIMA（2,1,0），RMSE 為 21.733，略低于之前測試的手動持久性模型，但可能會有所不同或可能沒有顯著差異。 ```py ... ARIMA(4, 1, 0) RMSE=24.802 ARIMA(4, 1, 1) RMSE=25.103 ARIMA(4, 2, 0) RMSE=27.089 ARIMA(4, 2, 1) RMSE=25.932 ARIMA(4, 2, 2) RMSE=25.418 Best ARIMA(2, 1, 0) RMSE=21.733 ``` 我們將繼續選擇這個 ARIMA（2,1,0）模型。 ### 6.3 查看殘留錯誤 對模型進行良好的最終檢查是檢查殘差預測誤差。 理想情況下，殘差的分布應該是具有零均值的高斯分布。 我們可以使用匯總統計和圖來檢查這一點，以研究 ARIMA（2,1,0）模型的殘差。以下示例計算并總結了殘差預測誤差。 ```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-forward validation history = [x for x in train] predictions = list() for i in range(len(test)): # predict model = ARIMA(history, order=(2,1,0)) model_fit = model.fit(trend='nc', 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) print(residuals.describe()) pyplot.figure() pyplot.subplot(211) residuals.hist(ax=pyplot.gca()) pyplot.subplot(212) residuals.plot(kind='kde', ax=pyplot.gca()) pyplot.show() ``` 首先運行該示例描述了殘差的分布。 我們可以看到分布具有正確的偏移，并且均值在 1.081624 處非零。 這或許表明預測存在偏見。 ```py count 35.000000 mean 1.081624 std 22.022566 min -52.103811 25% -16.202283 50% -0.459801 75% 12.085091 max 51.284336 ``` 還繪制了殘差的分布。 圖表顯示具有較長右尾的類高斯分布，提供了可能值得探索的功率變換的進一步證據。  剩余預測誤差密度圖 我們可以通過將每個預測的平均殘差誤差 1.081624 添加到偏差校正預測來使用此信息。 以下示例執行此偏差校正。 ```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-forward validation history = [x for x in train] predictions = list() bias = 1.081624 for i in range(len(test)): # predict model = ARIMA(history, order=(2,1,0)) model_fit = model.fit(trend='nc', disp=0) yhat = bias + float(model_fit.forecast()[0]) predictions.append(yhat) # observation obs = test[i] history.append(obs) # report performance mse = mean_squared_error(test, predictions) rmse = sqrt(mse) print('RMSE: %.3f' % rmse) # summarise residual errors residuals = [test[i]-predictions[i] for i in range(len(test))] residuals = DataFrame(residuals) print(residuals.describe()) # plot residual errors pyplot.figure() pyplot.subplot(211) residuals.hist(ax=pyplot.gca()) pyplot.subplot(212) residuals.plot(kind='kde', ax=pyplot.gca()) pyplot.show() ``` 預測的表現從 21.733 略微改善至 21.706，這可能是也可能不是很重要。 預測殘差的總結表明，平均值確實被移動到非常接近零的值。 ```py RMSE: 21.706 0 count 3.500000e+01 mean -3.537544e-07 std 2.202257e+01 min -5.318543e+01 25% -1.728391e+01 50% -1.541425e+00 75% 1.100347e+01 max 5.020271e+01 ``` 最后，剩余誤差的密度圖確實顯示向零的小偏移。  偏差校正殘差預測誤差密度圖 這種偏差校正是否值得，值得商榷，但我們現在將使用它。 ## 7.模型驗證 在開發模型并選擇最終模型后，必須對其進行驗證和最終確定。 驗證是流程的一個可選部分，但提供“最后檢查”以確保我們沒有被欺騙或誤導自己。 本節包括以下步驟： 1. **Finalize Model** ：訓練并保存最終模型。 2. **進行預測**：加載最終模型并進行預測。 3. **驗證模型**：加載并驗證最終模型。 ### 7.1 完成模型 最終確定模型涉及在整個數據集上擬合 ARIMA 模型，在這種情況下，在整個數據集的變換版本上。 一旦適合，模型可以保存到文件中供以后使用。 下面的示例在數據集上訓練 ARIMA（2,1,0）模型，并將整個擬合對象和偏差保存到文件中。 當前穩定版本的 statsmodels 庫（v0.6.1）中存在一個錯誤，當您嘗試從文件加載已保存的 ARIMA 模型時會導致錯誤。報告的錯誤是： ```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') # fit model model = ARIMA(X, order=(2,1,0)) model_fit = model.fit(trend='nc', disp=0) # bias constant, could be calculated from in-sample mean residual bias = 1.081624 # save model model_fit.save('model.pkl') numpy.save('model_bias.npy', [bias]) ``` 運行該示例將創建兩個本地文件： * _model.pkl_ 這是調用 _ARIMA.fit（）_ 的 ARIMAResult 對象。這包括系數和擬合模型時返回的所有其他內部數據。 * _model_bias.npy_ 這是存儲為單行，單列 NumPy 數組的偏差值。 ### 7.2 進行預測 一個自然的例子可能是加載模型并進行單一預測。 這是相對簡單的，包括恢復保存的模型和偏差并調用 _forecast（）_ 函數。 下面的示例加載模型，對下一個時間步進行預測，并打印預測。 ```py from pandas import Series from statsmodels.tsa.arima_model import ARIMAResults import numpy model_fit = ARIMAResults.load('model.pkl') bias = numpy.load('model_bias.npy') yhat = bias + float(model_fit.forecast()[0]) print('Predicted: %.3f' % yhat) ``` 運行該示例打印約 540 的預測。 ```py Predicted: 540.013 ``` 如果我們查看 _validation.csv_ ，我們可以看到下一個時間段第一行的值是 568.預測是在正確的球場。 ### 7.3 驗證模型 我們可以加載模型并以假裝操作方式使用它。 在測試工具部分中，我們將原始數據集的最后 10 年保存在單獨的文件中以驗證最終模型。 我們現在可以加載這個 _validation.csv_ 文件并使用它來查看我們的模型在“看不見的”數據上的真實程度。 我們可以通過兩種方式進行： * 加載模型并使用它來預測未來 10 年。超過前一年或兩年的預測將很快開始降低技能。 * 加載模型并以滾動預測方式使用它，更新每個時間步的變換和模型。這是首選方法，因為它是如何在實踐中使用此模型，因為它將實現最佳表現。 與前幾節中的模型評估一樣，我們將以滾動預測的方式進行預測。這意味著我們將在驗證數據集中逐步超過提前期，并將觀察結果作為歷史記錄的更新。 ```py from pandas import Series from matplotlib import pyplot from statsmodels.tsa.arima_model import ARIMA from statsmodels.tsa.arima_model import ARIMAResults from sklearn.metrics import mean_squared_error from math import sqrt import numpy # 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') bias = numpy.load('model_bias.npy') # make first prediction predictions = list() yhat = bias + float(model_fit.forecast()[0]) 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)): # predict model = ARIMA(history, order=(2,1,0)) model_fit = model.fit(trend='nc', disp=0) yhat = bias + float(model_fit.forecast()[0]) 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 預計為每人每天 16 升。這與 21 的預期誤差沒有太大的不同，但我希望它與簡單的持久性模型也沒有太大區別。 ```py >Predicted=540.013, Expected=568 >Predicted=571.589, Expected=575 >Predicted=573.289, Expected=579 >Predicted=579.561, Expected=587 >Predicted=588.063, Expected=602 >Predicted=603.022, Expected=594 >Predicted=593.178, Expected=587 >Predicted=588.558, Expected=587 >Predicted=588.797, Expected=625 >Predicted=627.941, Expected=613 RMSE: 16.532 ``` 還提供了與驗證數據集相比較的預測圖。 預測確實具有持久性預測的特征。這表明雖然這個時間序列確實有明顯的趨勢，但它仍然是一個相當困難的問題。  驗證數據集預測圖 ## 摘要 在本教程中，您使用 Python 發現了時間序列預測項目的步驟和工具。 我們在本教程中介紹了很多內容;特別： * 如何開發具有表現測量和評估方法的測試工具，以及如何快速開發基線預測和技能。 * 如何使用時間序列分析來提出如何最好地模擬預測問題的想法。 * 如何開發 ARIMA 模型，保存它，然后加載它以對新數據進行預測。 你是怎么做的？您對本教程有任何疑問嗎？ 在下面的評論中提出您的問題，我會盡力回答。