Pythonで時系列解析・超入門(その9)

テーブルデータ系モデルで構築する時系列予測モデル
(ランダムフォレスト)

Pythonで時系列解析・超入門(その9)テーブルデータ系モデルで構築する時系列予測モデル(ランダムフォレスト)

多くの人にとって馴染みがあるのは、時系列データ系の数理モデル(アルゴリズム)よりも、テーブルデータ系の数理モデル(アルゴリズム)の方です。

例えば、以下の数理モデル(アルゴリズム)テーブルデータ系のものです。

  • 線形回帰モデル(単回帰、重回帰、など)
  • 正則化回帰モデル(Ridge回帰、Lasso回帰、など)
  • 一般化線形モデル(GLMM
  • 一般化加法モデル(GAM
  • 階層線形モデル、マルチレベルモデル、一般化混合モデル
  • 決定木(ディシジョンツリー)
  • ランダムフォレスト
  • ブースティングモデル(AdaBoostXGBoostLightGBMなど)
  • ニューラルネットワークモデル

……などなど。

前回は、時系列特徴量付きデータセットを使い、ディシジョンツリー(決定木)で時系列予測モデルを構築しました。

Pythonで時系列解析・超入門(その8)テーブルデータ系モデルで構築する時系列予測モデル(ディシジョンツリー)

今回は、前回と同じ時系列特徴量付きデータセットを使い、ランダムフォレストで時系列予測モデルを構築します。

ランダムフォレストと言っても、以下の2パターンの方法で作ります。

  • ランダムフォレスト(デフォルト設定のまま)
  • ランダムフォレスト(Optunaでハイパーパラメータチューニング)

Optunaに関しては、以下の記事を参考にして下さい。

ランダムフォレストとは

ランダムフォレストは、複数のディシジョンツリー(決定木)の結果を束ね予測結果とします。結果を束ねるとは、平均値を計算したり、多数決を取ったリ、することです。

複数のディシジョンツリー(決定木)は、同じデータセットで作るわけではありません。それぞれ別のデータセットで作ります。

元のデータセットから複数のデータセットをランダムに作ります。そのデータセットごとにディシジョンツリー(決定木)を作ります。ランダムなデータセットで、複数のディシジョンツリー(決定木)を作るため、ランダムフォレストと呼ばれます。

元のデータセットから複数のデータセットの作り方は意外とシンプルで、元のデータセットの行(ケース)と列(変数)をランダムに選びます。

  • 行(ケース):元のデータセットから復元抽出(ブートストラップ法)でランダムにデータを抽出しデータセットを複数作ります
  • 列(変数):そのデータセットの説明変数(特徴量)をランダムに選択し、新たなデータセットとします

 

必要なライブラリーの読み込み

先ず、必要なライブラリーなどを読み込みます。

以下、コードです。

import numpy as np
import pandas as pd
import optuna

from sklearn.ensemble import RandomForestRegressor
from sklearn.inspection import PartialDependenceDisplay
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_absolute_percentage_error

import matplotlib.pyplot as plt

plt.style.use('ggplot') #グラフのスタイル
plt.rcParams['figure.figsize'] = [12, 9] # グラフサイズ設定

 

利用するデータ

今回利用するデータは、前回準備した時系列特徴量付きデータセットです。

以下からダウンロードできます。

dataset.csv
https://www.salesanalytics.co.jp/6ro8

 

このURLから直接データセットを読み込めます。

以下、コードです。

# データセットの読み込み
url='dataset.csv'
df=pd.read_csv(url,                         #読み込むデータのURL
               index_col='Month',           #変数「Month」をインデックスに設定
               parse_dates=True)            #インデックスを日付型に設定

df.head() #確認

 

以下、実行結果です。

 

グラフ化し確認します。

以下、コードです。

# プロット
df.plot()

plt.title('Passengers')                            #グラフタイトル
plt.ylabel('Monthly Number of Airline Passengers') #タテ軸のラベル
plt.xlabel('Month')                                #ヨコ軸のラベル
plt.show()

 

以下、実行結果です。

 

次に、読み込んだデータセットを、学習データとテストデータに分割します。

以下、コードです。

# 学習データ
train = df.iloc[:-12]

y_train = train['y']              #目的変数y
X_train = train.drop('y', axis=1) #説明変数X

# テストデータ
test = df.iloc[-12:]  #テストデータ

y_test = test['y']              #目的変数y
X_test = test.drop('y', axis=1) #説明変数X

 

グラフ化します。

以下、コードです。

# グラフ化
fig, ax = plt.subplots()

ax.plot(y_train.index, y_train.values, label="actual(train dataset)")
ax.plot(y_test.index, y_test.values, label="actual(test dataset)")

plt.legend()

 

以下、実行結果です。

 

学習データディシジョンツリー(決定木)モデルを構築し、構築したモデルをテストデータで精度検証します。

 

予測精度の評価指標

今回の予測精度の評価指標は、RMSE二乗平均平方根誤差、Root Mean Squared Error)とMAE平均絶対誤差、Mean Absolute Error)、MAPE平均絶対パーセント誤差、Mean absolute percentage error)を使います。 

以下の記号を使い精度指標の説明をします。

  • y_i^{actual} ・・・i番目の実測値
  • y_i^{pred} ・・・i番目の予測値
  • n ・・・実測値・予測値の数

■ 二乗平均平方根誤差RMSE、Root Mean Squared Error)

\sqrt{\frac{1}{n}\sum_{i=1}^n(y_i^{actual}-{y_i^{pred}})^2}

■ 平均絶対誤差MAE、Mean Absolute Error)

\frac{1}{n}\sum_{i=1}^n|y_i^{actual}-{y_i^{pred}}|

■ 平均絶対パーセント誤差MAPE、Mean absolute percentage error)

\frac{1}{n}\sum_{i=1}^n|\frac{y_i^{actual}-{y_i^{pred}}}{y_i^{actual}}|

 

ランダムフォレスト
(デフォルト設定のまま)

学習データを使って、ランダムフォレスト(デフォルト設定のまま)を学習します。

以下、コードです。

# 学習
regressor = RandomForestRegressor(random_state=123)
regressor.fit(X_train, y_train)

 

特徴量重要度(Feature Importances)を見てみます。

以下、コードです。

# 特徴量重要度(Feature Importances)
df_importance = pd.DataFrame(zip(X_train.columns, regressor.
                                 feature_importances_),
                             columns=["Features","Importance"])

df_importance = df_importance.sort_values("Importance",
                                          ascending=False)

df_importance #確認

 

以下、実行結果です。

 

特徴量重要度グラフ化します。

以下、コードです。

# グラフ化
df_importance.plot.bar(x='Features',y='Importance', rot=0)

 

以下、実行結果です。

 

部分従属プロット(Partial Dependence Plot)を見てみます。

以下、コードです。

 

テストデータで精度検証します。

以下、コードです。

# 予測
train_pred = regressor.predict(X_train)
test_pred = regressor.predict(X_test)

# 精度指標(テストデータ)
print('RMSE:')
print(np.sqrt(mean_squared_error(y_test, test_pred)))
print('MAE:')
print(mean_absolute_error(y_test, test_pred))
print('MAPE:')
print(mean_absolute_percentage_error(y_test, test_pred))

 

以下、実行結果です。

 

グラフ化します。

以下、コードです。

# グラフ化
fig, ax = plt.subplots()

ax.plot(y_train.index, y_train.values, label="actual(train dataset)")
ax.plot(y_test.index, y_test.values, label="actual(test dataset)")
ax.plot(y_train.index, train_pred, linestyle="dotted", lw=2,color="m")
ax.plot(y_test.index, test_pred, label="CART", linestyle="dotted", lw=2, color="m") 

plt.legend()

 

以下、実行結果です。

 

ランダムフォレスト
(Optunaでハイパーパラメータチューニング)

学習データを使って、ランダムフォレスト(Optunaでハイパーパラメータチューニング)を学習します。

先ずは、Optunaで最適なハイパーパラメータを探索の準備をします。

以下、コードです。

# ランダムフォレストモデル
regressor = RandomForestRegressor(random_state=123)

# ハイパーパラメータの探索の設定
params = {
    "n_estimators":optuna.distributions.IntUniformDistribution(10,10000,1),
    "max_depth":optuna.distributions.IntLogUniformDistribution(2,10,1),
}

# CV(クロスバリデーション)の設定
optuna_search = optuna.integration.OptunaSearchCV(regressor,
                                                  params,
                                                  cv=10,
                                                  n_jobs=-1,
                                                  n_trials=100,
                                                 )

 

では、Optunaで最適なハイパーパラメータを探索します。

以下、コードです。非常に時間がかかりますので、ご注意下さい。

# ハイパーパラメータの探索の実施
optuna_search.fit(X_train, y_train)

 

最適なハイパーパラメータを出力します。

以下、コードです。

# 最適解の出力
print(f"The best parameters are : \n {optuna_search.best_params_}")

 

以下、実行結果です。

 

この最適なハイパーパラメータを使い、ランダムフォレスト学習します。

以下、コードです。

# 学習
regressor = RandomForestRegressor(**optuna_search.best_params_)
regressor.fit(X_train, y_train)

 

特徴量重要度(Feature Importances)を見てみます。

以下、コードです。

# 特徴量重要度(Feature Importances)
df_importance = pd.DataFrame(zip(X_train.columns, regressor.feature_importances_),
                             columns=["Features","Importance"])

df_importance = df_importance.sort_values("Importance",
                                          ascending=False)

df_importance #確認

 

以下、実行結果です。

 

特徴量重要度グラフ化します。

以下、コードです。

# グラフ化
df_importance.plot.bar(x='Features',y='Importance', rot=0)

 

以下、実行結果です。

 

部分従属プロット(Partial Dependence Plot)を見てみます。

以下、コードです。

# 部分従属プロット(Partial Dependence Plot) 
PartialDependenceDisplay.from_estimator(regressor, X_train, [0,1,2])

 

以下、実行結果です。

 

テストデータで精度検証します。

以下、コードです。

# 予測
train_pred = regressor.predict(X_train)
test_pred = regressor.predict(X_test)

# 精度指標(テストデータ)
print('RMSE:')
print(np.sqrt(mean_squared_error(y_test, test_pred)))
print('MAE:')
print(mean_absolute_error(y_test, test_pred))
print('MAPE:')
print(mean_absolute_percentage_error(y_test, test_pred))

 

以下、実行結果です。

 

グラフ化します。

以下、コードです。

# グラフ化
fig, ax = plt.subplots()

ax.plot(y_train.index, y_train.values, label="actual(train dataset)")
ax.plot(y_test.index, y_test.values, label="actual(test dataset)")
ax.plot(y_train.index, train_pred, linestyle="dotted", lw=2,color="m")
ax.plot(y_test.index, test_pred, label="CART", linestyle="dotted", lw=2, color="m") 

plt.legend()

 

以下、実行結果です。

 

次回

今回は、前回と同じ時系列特徴量付きデータセットを使い、ランダムフォレストで時系列予測モデルを構築しました。

次回は、今回と同じ時系列特徴量付きデータセットを使い、XGBoostで時系列予測モデルを構築します。