多くの人にとって馴染みがあるのは、時系列データ系の数理モデル（アルゴリズム）よりも、テーブルデータ系の数理モデル（アルゴリズム）の方です。
例えば、以下の数理モデル（アルゴリズム）はテーブルデータ系のものです。
- 線形回帰モデル（単回帰、重回帰、など）
- 正則化回帰モデル（Ridge回帰、Lasso回帰、など）
- 一般化線形モデル（GLMM）
- 一般化加法モデル（GAM）
- 階層線形モデル、マルチレベルモデル、一般化混合モデル
- 決定木（ディシジョンツリー）
- ランダムフォレスト
- ブースティングモデル（AdaBoost、XGBoost、LightGBMなど）
- ニューラルネットワークモデル
……などなど。
前回は、時系列特徴量付きデータセットを使い、ランダムフォレストで時系列予測モデルを構築しました。
今回は、前回と同じ時系列特徴量付きデータセットを使い、XGBoostで時系列予測モデルを構築します。
XGBoostと言っても、以下の2パターンの方法で作ります。
- XGBoost（デフォルト設定のまま）
- XGBoost（Optunaでハイパーパラメータチューニング）
Optunaに関しては、以下の記事を参考にして下さい。
XGBoostとは
ランダムフォレストと似たような数理モデルに、ブースティングモデルというものがあります。
ブースティングモデルの中に、XGBoostという数理モデルがあります。どれも決定木（ディシジョンツリー）をたくさん作る数理モデルです。作り方が異なります。
- 並列アンサンブル：ランダムフォレスト
- 直列アンサンブル：XGBoost
アンサンブル学習とは、複数の学習器（数理モデル）を構築し融合させ1つの数理モデルを生成する学習方法です。アンサンブルを直訳すると合奏や合唱などになります。
決定木（ディシジョンツリー）を活用した並列アンサンブル学習で構築する数理モデルの代表がランダムフォレストで、直列アンサンブル学習で構築する数理モデルの代表がブースティングモデルです。
要は、XGBoostは決定木（ディシジョンツリー）を活用した直列アンサンブル学習で構築する数理モデルの1つということです。
前回の復習になります。
並列アンサンブルであるランダムフォレストは、元のデータセットから複数の独立したデータセットをたくさん作り、そのデータセットごとに決定木（ディシジョンツリー）を構築します。それぞれの決定木の予測値を集約し最終的な予測値にします。
並列アンサンブルに対し、直列アンサンブルはどうちがうのでしょうか？
直列アンサンブルであるXGBoostなどは、決定木（ディシジョンツリー）を順番に作っていきます。先ず、とりあえず決定木を作ります。次に、その決定木の誤りを修正するような決定木を作ります。その次に、その決定木の誤りを修正するような決定木を作ります。このように直列的に決定木をどんどん作っていきます。
必要なライブラリーの読み込み
先ず、必要なライブラリーなどを読み込みます。
以下、コードです。
import numpy as np
import pandas as pd
import optuna
import xgboost as xgb
from sklearn import model_selection
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_absolute_percentage_error
from sklearn.inspection import PartialDependenceDisplay
from matplotlib import 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()
以下、実行結果です。
学習データでXGBoostモデルを構築し、構築したモデルをテストデータで精度検証します。
予測精度の評価指標
今回の予測精度の評価指標は、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}}|
XGBoost（デフォルト設定のまま）
学習データを使って、XGBoost（デフォルト設定のまま）を学習します。
以下、コードです。
# 学習 regressor = xgb.XGBRegressor() regressor.fit(X_train.values, y_train.values)
特徴量重要度（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.values,
features=[0,1,2])
以下、実行結果です。
テストデータで精度検証します。
以下、コードです。
# 予測
train_pred = regressor.predict(X_train.values)
test_pred = regressor.predict(X_test.values)
# 精度指標（テストデータ）
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="XGBoost", linestyle="dotted", lw=2, color="m") plt.legend()
以下、実行結果です。
XGBoost（Optunaでハイパーパラメータチューニング）
学習データを使って、XGBoost（Optunaでハイパーパラメータチューニング）を学習します。
先ずは、Optunaで最適なハイパーパラメータを探索の準備をします。
以下、コードです。
# XGBoostモデル
regressor = xgb.XGBRegressor()
# 目的関数の設定
def objective(trial):
#ハイパーパラメータの探索の設定
params = {
'max_depth': trial.suggest_int('max_depth', 1, 9),
'n_estimators': trial.suggest_int('n_estimators', 10, 1000),
'learning_rate': trial.suggest_loguniform('learning_rate', 1e-8, 1.0)
}
regressor.set_params(**params)
#CV（クロスバリデーション）の設定
score = model_selection.cross_val_score(regressor,
X_train,
y_train,
n_jobs=-1,
cv=10,
scoring = 'neg_mean_squared_error'
)
val = - score.mean()
return val
では、Optunaで最適なハイパーパラメータを探索します。
以下、コードです。非常に時間がかかりますので、ご注意下さい。
# ハイパーパラメータの探索の実施 study = optuna.create_study(direction="minimize") study.optimize(objective, n_trials=100)
最適なハイパーパラメータを出力します。
以下、コードです。
# 最適解の出力
print(f"The best parameters are : \n {study.best_params}")
以下、実行結果です。
この最適なハイパーパラメータを使い、XGBoostモデルを学習します。
以下、コードです。
# 最適パラメータでモデル学習 regressor = xgb.XGBRegressor(**study.best_params) regressor.fit(X_train.values, y_train.values)
特徴量重要度（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.values,
features=[0,1,2])
以下、実行結果です。
テストデータで精度検証します。
以下、コードです。
# 予測
train_pred = regressor.predict(X_train.values)
test_pred = regressor.predict(X_test.values)
# 精度指標（テストデータ）
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()
以下、実行結果です。
次回
今回も含め何回かに分け、テーブルデータ系モデルで時系列予測モデルを構築する方法について説明してきました。
テーブルデータ系モデルで予測し精度検証してきました。
もしかしたら、ある予測タイミング（ある1時点）で複数先予測（Multi-Step ahead prediction）をしているかのように見えていたかもしれません。実は予測するタイミングを1時点ずらしながら1期先予測（1-Step ahead prediction）を複数回実施していただけです。
複数先予測（Multi-Step ahead prediction）とは、例えば、2022年5月18日（ある1時点の予測タイミング）に、5月19日から5月25日の各日の売上などを予測する、という感じです。
予測するタイミングを1時点ずらしながら1期先予測（1-Step ahead prediction）とは、例えば、2022年5月18日（予測タイミング）に次の日の5月19日の売上を予測、2022年5月19日（予測タイミング）に次の日の5月20日の売上を予測、2022年5月20日（予測タイミング）に次の日の5月21日の売上を予測、・・・みたいな感じのことを言っています。
実務的には、ある予測タイミング（ある1時点）で複数先予測（Multi-Step ahead prediction）を実施することが多いです。そのため、時系列系の数理モデルのパッケージの場合、もとから複数先予測（Multi-Step ahead prediction）できるようなっています。
次回は、今回と同じ時系列特徴量付きデータセットを使い、テーブルデータ系モデルで複数先予測（Multi-Step ahead prediction）方法について説明し、よく利用される「1期先予測モデルを1つ作り再帰的に利用する方法」について説明します。
