モモノキ＆ナノネと一緒にPythonで時系列分析を覚えよう（10）
時系列分析の実践練習（予測の95%信頼区間をグラフに表示する）
前回の時系列分析の続きで、予測の95%信頼区間をグラフに追加してみよう。
予測範囲みたいなやつ？
うん、予測の95%信頼区間を帯状に塗り潰す予定だよ。コレログラムのグラフ表示みたいに。
ナノネ、前回のコード呼び出してみて。予測結果のグラフが表示できるように必要箇所だけピックアップしてね。
# 飛行機乗客数の時系列データ（前回使用コード抜粋 -> 次回SARIAモデル実践に利用）
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use('ggplot') # グラフのデザイン変更（お好みで利用）
plt.xkcd() # グラフのデザイン変更（お好みで利用、日本語は未対応）
import statsmodels.api as sm # version 0.8.0以上
# CSVファイル読み込み（事前にダウンロードしたCSVファイルを利用）
df = pd.read_csv('AirPassengers.csv')
# Pandas.Seriesにデータを格納（データに乗客数、インデックスは日付）
passengers = pd.Series(df['#Passengers'], dtype='float')
passengers.index = pd.to_datetime(df['Month'])
# オリジナル passengers 1949-01 ～ 1960-12
passengers_train2 = passengers['1949-01':'1959-12'] # モデル作成用データ（訓練）1年テスト用残し
passengers_test2 = passengers['1960-01':'1960-12'] # テスト用データ1年分
# SRIMAモデル（テストデータ1年を除いてモデル作成）
sarimax_train = sm.tsa.SARIMAX(passengers_train2,
order=(3, 1, 3),
seasonal_order=(0, 1, 1, 12),
enforce_stationarity = False,
enforce_invertibility = False
).fit()
sarimax_train2_pred = sarimax_train.predict('1959-12', '1960-12') # テストデータ1年分予測
plt.figure(figsize=(8, 4))
plt.plot(passengers, label="actual")
plt.plot(sarimax_train2_pred, c="b", label="model-pred", alpha=0.7)
plt.legend(loc='best')
plt.figure(figsize=(8, 4))
plt.plot(passengers_test2, label="actual") # 正解
plt.plot(sarimax_train2_pred, c="b", label="predict", alpha=0.7) # 予測
plt.legend(loc='best')
OK、必要なとこだけピックアップした。予測まで実行しているけど必要なコードは以外に少ない。大半はデータ準備とグラフ描画。
時系列データの事前分析やモデルのパラメータチェックなどを除くと、予測に必要な処理の実装部分のコードは比較的短く書けるね。
信頼区間の範囲をグラフに表示するには、どうやればいいの？
区間の上下限値が取得できれば、なんとかなりそう。値の取得方法が分からないからグーグル先生に聞いてみる。
ポチポチ？・・・・・・・ポチポチ？ ．．．．．．．．．．．。
ポチポチ？.......ポチポチ？.......ポチっと、これでいけるかも。
ずいぶん苦戦したね。
get_prediction()とconf_int()を使えば信頼区間の値が取れそうだよ。
StatsModels公式ドキュメント（conf_int）
http://www.statsmodels.org/dev/generated/statsmodels.regression.linear_model.OLSResults.conf_int.html 外部サイト
①モデル変数名.get_prediction(予測開始, 予測終了)として、結果を変数で受ける。結果はstatsmodels.tsa.statespace.mlemodel.PredictionResultsWrapperという形式みたいだね。
predict_dy = sarimax_train.get_prediction(start ='1959-12', end='1960-12')
type(predict_dy)
②受けった結果変数.conf_int(alpha=0.05)とすれば、95%信頼区間の上下値がDataFrameの形式で取得できるみたい。alphaはデフォルトで0.05だから省略可、任意の数値も指定できるよ。
predict_dy_ci = predict_dy.conf_int(alpha=0.05) # defalut alpah=0.05 :returns a 95% confidence interval
type(predict_dy_ci)
predict_dy_ci # lower, upper取得
DataFrameの中身は上の通りで、指定した期間の信頼区間が取得できてるよ。データの1列目がlower値、2列目がupper値だね。
あと予測の平均値は.predicted_meanで取得できるよ。これまで使ってきた.predict()の結果と同じだよ。
predict_dy.predicted_mean # mean取得
lowerとupper値を使ってグラフに帯状に塗り潰すには？
Matplotlibのfill_between()を使えば簡単にできるよ。引数にx値（インデックス）、y-lower値、 y-upper値を渡せば、いい感じに塗り潰してくれるよ。色とか透明度はお好みでしてね。
plt.figure(figsize=(8, 4))
plt.plot(passengers_test2, label="Actual")
plt.plot(predict_dy.predicted_mean, c="b", label="model-pred", alpha=0.7)
plt.fill_between(predict_dy_ci.index, predict_dy_ci.iloc[:, 0], predict_dy_ci.iloc[:, 1], color='g', alpha=0.2)
plt.legend(loc='upper left')
できた。緑で薄く塗り潰したところが、95%信頼区間だね。
predict2_dy = sarimax_train.get_prediction(start ='1959-12', end='1961-12') # 2年分の未来予測
predict2_dy_ci = predict2_dy.conf_int()
plt.figure(figsize=(12, 4))
plt.plot(passengers['1959-01':], label="actual")
plt.plot(predict2_dy.predicted_mean, c="b", linestyle='--', label="model-pred", alpha=0.7)
plt.fill_between(predict2_dy_ci.index, predict2_dy_ci.iloc[:, 0], predict2_dy_ci.iloc[:, 1], color='g', alpha=0.2)
plt.legend(loc='upper left')
予測を1年増やして、表示期間も少しか変えてみた。
初めてにしては、それっぽい予測ができたね。
うまく予測できるか分からないけど、他のデータ使ったモデル作成にも挑戦してみてね。
飛行機乗客数を使った時系列分析の練習は以上だよ。
理解不足で説明が不足してたり正確ではないとこもあるので、もっと詳しい書籍やサイト利用して理解を深めてね。あと、利用しているパッケージの仕様が大きくかわるとコードが動作しなくなる可能性もあるので、その時は代替方法を検索で調べてみてね。
ラジャー、お疲れさまでした。
お疲れさまでした。またね。
# 実行環境
# import pkg_resources
# for dist in pkg_resources.working_set:
# print(dist.project_name, dist.version)
・Python 3.6.4 |Anaconda ・notebook 5.4.0 ・matplotlib 2.1.2 ・pandas 0.22.0 ・numpy 1.14.0 ・statsmodels 0.8.0
# モモノキ＆ナノネと一緒にPythonで時系列分析を覚えよう
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use('ggplot')
plt.xkcd()
import statsmodels.api as sm # version 0.8.0以上
# CSVファイル読み込み（事前にダウンロードしたCSVファイルを利用）
df = pd.read_csv('AirPassengers.csv')
# Pandas.Seriesにデータを格納
passengers = pd.Series(df['#Passengers'], dtype='float') # 乗客数
passengers.index = pd.to_datetime(df['Month']) # 日付インデックス
# 元データ分割：passengers 1949-01 ～ 1960-12
passengers_train2 = passengers['1949-01':'1959-12'] # 訓練データ（モデル作成用）
passengers_test2 = passengers['1960-01':'1960-12'] # テストデータ1年分（検証用）
# SRIMAモデル
sarimax_train = sm.tsa.SARIMAX(passengers_train2,
order=(3, 1, 3),
seasonal_order=(0, 1, 1, 12),
enforce_stationarity = False,
enforce_invertibility = False
).fit()
# モデルによる予測
predict2_dy = sarimax_train.get_prediction(start ='1959-12', end='1961-12') # 2年分の未来予測
predict2_dy_ci = predict2_dy.conf_int() # 信頼区間取得
# グラフ表示
plt.figure(figsize=(12, 4))
plt.plot(passengers['1959-01':], label="actual") # 実データプロット
plt.plot(predict2_dy.predicted_mean, c="b", linestyle='--', label="model-pred", alpha=0.7) # 予測プロット
# 予測の95%信頼区間プロット（帯状）
plt.fill_between(predict2_dy_ci.index, predict2_dy_ci.iloc[:, 0], predict2_dy_ci.iloc[:, 1], color='g', alpha=0.2)
plt.legend(loc='upper left')
0 件のコメント :
コメントを投稿