Data Driven Knowledgebase

!pip install pmdarima

!pip install optuna

!pip install japanize_matplotlib



import japanize_matplotlib as plt

import pandas as pd

import numpy as np

from scipy import signal

from statsmodels.tsa.seasonal import seasonal_decompose

from statsmodels.tsa.seasonal import STL

from statsmodels.tsa.stattools import adfuller

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

import pmdarima as pm

from pmdarima import utils

from pmdarima import arima

from pmdarima import model_selection

from statsmodels.tsa.statespace.sarimax import SARIMAX

import optuna

from prophet import Prophet

from prophet.diagnostics import cross_validation

from prophet.diagnostics import performance_metrics



url='https://www.salesanalytics.co.jp/591h' #データセットのあるURL

df=pd.read_csv(url, #読み込むデータのURL

index_col='Month', #変数「Month」をインデックスに設定

parse_dates=True) #インデックスを日付型に設定



# データを読み込む

url='https://www.salesanalytics.co.jp/591h'

df=pd.read_csv(url, index_col='Month', parse_dates=True)



# matplotlibでグラフのスタイルを設定

plt.style.use('ggplot')



# グラフの描画

df.plot(figsize=(10, 5)) # figsizeでグラフのサイズを指定



# グラフのタイトルと軸ラベルを日本語で設定

plt.title('飛行機の乗客数推移') # グラフのタイトルを設定

plt.xlabel('日付') # 横軸（X軸）のラベルを「日付」とする

plt.ylabel('客数') # 縦軸（Y軸）のラベルを「値」とする



# グラフを表示

plt.show()

今回はオープンソースである乗客数推移で練習します。

import pandas as pd

from statsmodels.tsa.seasonal import STL



# STL分解を行う

stl = STL(df['Passengers'], seasonal=13)

result = stl.fit()



# 分解された成分をプロットする

seasonal, trend, resid = result.seasonal, result.trend, result.resid





plt.figure(figsize=(8, 6))



# トレンド成分

plt.subplot(311)

plt.plot(trend)

plt.title('トレンド成分')



# 季節成分

plt.subplot(312)

plt.plot(seasonal)

plt.title('季節成分')



# 残差成分

plt.subplot(313)

plt.plot(resid)

plt.title('残差成分')



plt.tight_layout()

plt.show()

このコードでSTL分解を行い、以下のように考察しました。

plot_acf(df.Passengers, lags=20)

plot_pacf(df.Passengers, lags=20)

from sklearn.metrics import mean_squared_error,mean_absolute_error,mean_absolute_percentage_error



df_train, df_test = model_selection.train_test_split(df, test_size=12)



arima_model = pm.auto_arima(df_train,

seasonal=False,

m=12,

trace=True,

n_jobs=-1,

maxiter=10)



# 予測

##学習データの期間の予測値

train_pred = arima_model.predict_in_sample()

##テストデータの期間の予測値

test_pred, test_pred_ci = arima_model.predict(

n_periods=df_test.shape[0],

return_conf_int=True

)

# テストデータで精度検証

print('RMSE:')

print(np.sqrt(mean_squared_error(df_test, test_pred)))

print('MAE:')

print(mean_absolute_error(df_test, test_pred))

print('MAPE:')

print(mean_absolute_percentage_error(df_test, test_pred))



# グラフ化

fig, ax = plt.subplots()

ax.plot(df_train[24:].index, df_train[24:].values, label="actual(train dataset)")

ax.plot(df_test.index, df_test.values, label="actual(test dataset)", color="gray")

ax.plot(df_train[24:].index, train_pred[24:], color="c")

ax.plot(df_test.index, test_pred, label="auto ARIMA", color="c")

ax.fill_between(

df_test.index,

test_pred_ci[:, 0],

test_pred_ci[:, 1],

color='c',

alpha=.2)

ax.legend()

これを実行すればこのような実行結果になります。

最後に残差をプロットして予測精度とモデルの適合具合を確認しましょう。

# ここから残差の分布を描画します。

residuals = df_test.values.flatten() - test_pred

plt.figure()

plt.hist(residuals, bins=6, alpha=0.7, color='blue')

plt.xlabel('残差')

plt.ylabel('頻度')

plt.show()

# モデルのサマリーを表示する

print(arima_model.summary())