時系列パート(基礎編:予測モデルの構築)#
ここでは,カフェの顧客データ(cafe_customers.csv)を使用して,時系列予測モデルの構築から評価まで実践的なスキルを学習します.
必要なライブラリのインポート#
pandas, numpy: データ処理の基本ライブラリ
matplotlib: グラフ作成ライブラリ
sklearn: 機械学習ライブラリ(Scikit-learn)
TimeSeriesSplit: 時系列データ専用のクロスバリデーションGridSearchCV: ハイパーパラメータの自動調整各種回帰アルゴリズム: LinearRegression, Ridge, RandomForest, GradientBoosting
評価指標: RMSE, MAE, R²など
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from datetime import timedelta
from tqdm import tqdm
# 警告メッセージを非表示にするライブラリ・設定
import warnings
warnings.filterwarnings('ignore')
# 機械学習ライブラリ
from sklearn.model_selection import TimeSeriesSplit, GridSearchCV
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
データの読み込み#
# データの読み込み
feature_df = pd.read_csv('../data/processed/feature_engineered_data.csv')
# Timestampをdatetime型に変換
feature_df['Timestamp'] = pd.to_datetime(feature_df['Timestamp'])
feature_df.head()
時系列データの適切な分割#
なぜ時系列データは特別な分割が必要なのか?#
通常の機械学習では,データをランダムに訓練用・テスト用に分割します.
しかし,時系列データでは以下の理由から時間順序を保った分割が必要です:
データリークの防止: 未来の情報を使って過去を予測してしまう問題を避ける
現実的な評価: 実際の運用では常に過去のデータで未来を予測する
時間依存性の考慮: 時系列データには時間的な依存関係がある
分割方法:
訓練データ: 過去のデータ(データ全体の80%)
テストデータ: 最新のデータ(データ全体の20%)
この方法により,モデルが実際の運用環境に近い条件で評価されます.
# テストデータの日数を決定(データ全体の20%)
total_days = (feature_df['Timestamp'].max() - feature_df['Timestamp'].min()).days
test_days = total_days // 5
# 時系列データの適切な分割
split_date = feature_df['Timestamp'].max() - timedelta(days=test_days)
trainval_df = feature_df[feature_df['Timestamp'] <= split_date].copy() # 訓練・検証データ
test_df = feature_df[feature_df['Timestamp'] > split_date].copy() # テストデータ
print("訓練・テストデータの期間")
print(f"訓練期間 : {trainval_df['Timestamp'].min()} - {trainval_df['Timestamp'].max()}")
print(f"テスト期間: {test_df['Timestamp'].min()} - {test_df['Timestamp'].max()}")
print("\n訓練・テストデータの行数")
print(f"訓練データ : {len(trainval_df)} 行")
print(f"テストデータ: {len(test_df)} 行")
説明変数と目的変数の分離#
機械学習では,データを以下のように分けて考えます:
説明変数(X): 予測に使用する特徴量(時間,曜日,ラグ特徴量など)
目的変数(y): 予測したい値(この場合は顧客数)
重要なポイント:
Timestampカラムは予測には直接使用しないため除外Customersカラムが予測対象なので,これも説明変数から除外これにより,モデルは時間情報と顧客数以外の全ての特徴量を学習に使用
# 説明変数Xと目的変数yを分離
feature_cols = [col for col in feature_df.columns if col not in ['Timestamp', 'Customers']]
X_train = trainval_df[feature_cols] # 訓練・検証データの説明変数
y_train = trainval_df['Customers'] # 訓練・検証データの目的変数
X_test = test_df[feature_cols] # テストデータの説明変数
y_test = test_df['Customers'] # テストデータの目的変数
モデルの構築と比較#
複数のアルゴリズムを比較する理由#
異なる機械学習アルゴリズムにはそれぞれ特徴があります:
線形回帰(Linear Regression)
最もシンプルな手法
特徴量と目的変数の線形関係を学習
解釈しやすいが,複雑なパターンは捉えられない
Ridge回帰
線形回帰に正則化を追加
過学習を防ぎ,汎化性能を向上
多重共線性の問題にも対処
ランダムフォレスト(Random Forest)
複数の決定木を組み合わせた手法
非線形関係も学習可能
特徴量重要度が分析できる
勾配ブースティング(Gradient Boosting)
弱学習器を順次改善していく手法
高い予測精度を実現
ハイパーパラメータの調整が重要
# 複数のモデルを定義
models = {
'Linear Regression': LinearRegression(),
'Ridge Regression' : Ridge(alpha=1.0),
'Random Forest' : RandomForestRegressor(n_estimators=100, random_state=42),
'Gradient Boosting': GradientBoostingRegressor(n_estimators=100, random_state=42)
}
# モデルの訓練と評価
results = {} # 評価結果を格納する辞書
predictions = {} # 予測結果を格納する辞書
for name, model in tqdm(models.items()):
# モデル訓練
model.fit(X_train, y_train)
# 予測
y_pred = model.predict(X_test)
predictions[name] = y_pred
# 評価指標計算(RMSE, MAE, R²)
rmse = np.sqrt(mean_squared_error(y_test, y_pred))
mae = mean_absolute_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)
results[name] = {
'RMSE↓': rmse,
'MAE↓' : mae,
'R²↑' : r2
}
評価指標の理解#
モデルの性能を測るため,以下の指標を使用します:
RMSE(Root Mean Squared Error)
予測誤差の二乗平均の平方根
大きな誤差により敏感
単位は目的変数と同じ(この場合は「人」)
値が小さいほど良い ↓
MAE(Mean Absolute Error)
予測誤差の絶対値の平均
外れ値に対してRMSEより頑健
解釈しやすい(平均的な誤差)
値が小さいほど良い ↓
R²(決定係数)
モデルがデータの分散をどれだけ説明できるかを示す
0〜1の値(1に近いほど良い予測)
値が大きいほど良い ↑
これらの指標を組み合わせることで,モデルの性能を多角的に評価できます.
予測結果の可視化#
# 学習結果の表示
results_df = pd.DataFrame(results).T
results_df.head().round(4)
実際の値 vs 予測値で可視化#
model_name = "Random Forest"
# 可視化(折れ線グラフ)
plt.figure(figsize=(15, 6))
plt.plot(test_df['Timestamp'], y_test, label='Actual', linewidth=1, alpha=0.8)
plt.plot(test_df['Timestamp'], predictions[model_name], label=f'{model_name} Prediction', linewidth=1, alpha=0.8)
plt.title(f'Time Series Prediction: {model_name}')
plt.xlabel('Date')
plt.ylabel('Number of Customers')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 可視化(散布図)
plt.figure(figsize=(15, 6))
plt.scatter(y_test, predictions[model_name], alpha=0.6, s=20)
plt.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], 'r--', lw=2)
plt.xlabel('Actual Values')
plt.ylabel('Predicted Values')
plt.title(f'Actual vs Predicted: {model_name}')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
散布図(Actual vs Predicted)の見方#
理想的な予測では,すべての点が赤い対角線(y=x)上に乗ります.
散布図から分かること:
点が対角線に近い: 予測精度が高い
点が対角線より上: 実際の値より低く予測している(過小評価)
点が対角線より下: 実際の値より高く予測している(過大評価)
点のばらつき: 予測の安定性(ばらつきが小さい方が良い)
残差分析#
# 残差を計算(予測値 - 実測値)
residuals = y_test - predictions[model_name]
# 可視化(散布図)
plt.figure(figsize=(15, 6))
plt.scatter(predictions[model_name], residuals, alpha=0.6, s=20)
plt.axhline(y=0, color='r', linestyle='--')
plt.xlabel('Predicted Values')
plt.ylabel('Residuals')
plt.title('Residuals vs Predicted')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
残差分析とは#
残差(Residuals) = 実際の値 - 予測値
残差分析は,モデルの予測誤差のパターンを調べる重要な手法です.
理想的な残差プロット:
残差が水平線(y=0)の周りにランダムに散らばっている
特定のパターンや傾向が見られない
分散が一定(等分散性)
問題のあるパターン:
曲線状のパターン: モデルが非線形関係を捉えきれていない
扇形のパターン: 異分散性(予測値によって誤差の大きさが変わる)
周期的なパターン: 季節性や周期性を捉えきれていない
ハイパーパラメータチューニング(グリッドサーチ)#
ハイパーパラメータとは,学習前に設定する必要があるパラメータです(学習によって更新されるパラメータとは異なります).
主要なハイパーパラメータ:
Random Forest:
n_estimators: 決定木の数(多いほど性能向上,計算時間増加)max_depth: 木の最大深度(深いほど複雑なパターンを学習,過学習リスク増)min_samples_split: 分割に必要な最小サンプル数(大きいほど過学習を防ぐ)
Gradient Boosting:
n_estimators: ブースティング段階数learning_rate: 学習率(小さいほど慎重に学習,計算時間増加)max_depth: 各弱学習器の最大深度
Ridge Regression:
alpha: 正則化の強さ(大きいほど過学習を防ぐ,小さいほど訓練データに適合)
チューニングの方法:
グリッドサーチ: 全ての組み合わせを試す(確実だが時間がかかる)
ランダムサーチ: ランダムに組み合わせを試す(効率的)
ベイズ最適化: より賢い探索方法(高度な手法)
# モデルごとのハイパーパラメータグリッド
if model_name == 'Random Forest':
param_grid = {
'n_estimators' : [50, 100, 200],
'max_depth' : [None, 10, 20],
'min_samples_split': [2, 5]
}
base_model = RandomForestRegressor(random_state=42)
elif model_name == 'Gradient Boosting':
param_grid = {
'n_estimators' : [50, 100, 200],
'learning_rate': [0.05, 0.1, 0.2],
'max_depth' : [3, 5, 7]
}
base_model = GradientBoostingRegressor(random_state=42)
elif model_name == 'Ridge Regression':
param_grid = {
'alpha': [0.1, 1.0, 10.0, 100.0]
}
base_model = Ridge()
else: # Linear Regression
param_grid = {}
base_model = LinearRegression()
時系列クロスバリデーション(TimeSeriesSplit)#
通常のクロスバリデーションは,データをランダムに分割して交差検証を行います.しかし,時系列データでは時間順序を保った分割が必要です.
TimeSeriesSplitの仕組み:
Split 1: Train [1,2,3] → Test [4]
Split 2: Train [1,2,3,4] → Test [5]
Split 3: Train [1,2,3,4,5] → Test [6]
特徴:
訓練データは常に過去,テストデータは未来
各分割で訓練データが増加(より現実的)
データリークを完全に防ぐ
パラメータ:
n_splits=3: 3回の交差検証を実行各分割で異なる時期のデータでモデルを評価
最終的には全分割の平均スコアで最適パラメータを決定
注意点:
通常のKFoldより時間がかかる
データ量が少ないと訓練データが不足する可能性
# 時系列クロスバリデーション
tscv = TimeSeriesSplit(n_splits=3)
# グリッドサーチ実行(数分かかります)
grid_search = GridSearchCV(base_model, param_grid, cv=tscv, scoring='neg_mean_squared_error', n_jobs=-1)
grid_search.fit(X_train, y_train)
# 最良パラメータで再予測
best_tuned_model = grid_search.best_estimator_
y_pred_tuned = best_tuned_model.predict(X_test)
# チューニング後の評価
rmse_tuned = np.sqrt(mean_squared_error(y_test, y_pred_tuned))
rmse_original = results[model_name]['RMSE↓']
print(f"チューニング前RMSE: {rmse_original:.4f}")
print(f"チューニング後RMSE: {rmse_tuned:.4f}")
特徴量重要度の分析(決定木ベースモデルのみ)#
特徴量重要度とは#
特徴量重要度は,各特徴量が予測にどれだけ貢献しているかを数値化したものです.
Random ForestとGradient Boostingの特徴量重要度:
各決定木における特徴量の分割回数と改善量に基づいて計算
値が大きいほど,その特徴量が予測により重要
全ての重要度の合計は1.0になる
活用方法:
特徴選択: 重要度の低い特徴量を除去してモデルを簡素化
ドメイン知識の検証: 業務的に重要だと思われる特徴量が実際に重要か確認
新しい特徴量のヒント: 重要な特徴量を参考に新しい特徴量を作成
# 特徴量重要度を取得
importances = best_tuned_model.feature_importances_
# 重要度をDataFrameに整理
feature_importance_df = pd.DataFrame({
'feature' : feature_cols,
'importance': importances
}).sort_values('importance', ascending=False)
# 可視化(全特徴量)
plt.figure(figsize=(15, 6))
plt.barh(range(len(feature_importance_df)), feature_importance_df['importance'])
plt.yticks(range(len(feature_importance_df)), feature_importance_df['feature'])
plt.xlabel('Feature Importance')
plt.title(f'Feature Importances - {model_name}')
plt.gca().invert_yaxis()
plt.tight_layout()
plt.show()
# 可視化(重要度上位15特徴量)
top_features = feature_importance_df.head(15)
plt.figure(figsize=(10, 8))
plt.barh(range(len(top_features)), top_features['importance'])
plt.yticks(range(len(top_features)), top_features['feature'])
plt.xlabel('Feature Importance')
plt.title(f'Top 15 Feature Importances - {model_name}')
plt.gca().invert_yaxis()
plt.tight_layout()
plt.show()