【Kaggle的】LightGBM とOptunaで住宅価格予測【RMSE: 0.440】

# はじめに
前々回、california_housingデータセットを使って、回帰を行った。今回はOptunaを導入して、ハイパーパラメータの最適化を行う。

# できたこと
- データセット: california_housing
- 問題: 回帰
- 機械学習モデル: LightGBM
- 前処理, 特徴量エンジニアリング: 何もせず
- 交差検証: K-Fold (5-Fold)
- ハイパーパラメータチューニング
    - ライブラリ: Optuna
    - 項目: num_leaves, learning_rate, feature_fraction, etc...
    - 試行回数: 50
- 評価指標: RMSE
- CVスコア: 0.4494 (前々回は、0.4739)
- テストスコア: 0.4402 (前々回は、0.4622)

```
=== RMSE for test data ===
0.44027
```

前々回（チューニングなし）に比べて、わずかにスコアアップした。よかった。

# できてないこと
EDAは割愛する。

# 実装
## ライブラリ
```Python
# 関連ライブラリをインポート
import lightgbm as lgb
import numpy as np
import optuna
import pandas as pd
import random

# scikit-learn関連をインポート
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split, KFold, StratifiedKFold
from sklearn.metrics import mean_squared_error

# warningを非表示
import warnings
warnings.filterwarnings('ignore')
```

出力を一定にするため、乱数を固定する。
```Python
SEED = 42

# 参考: 乱数固定
def seed_everything(seed=42):
    random.seed(seed)
    np.random.seed(seed)
    #os.environ['PYTHONHASHSEED'] = str(seed)
    # torch.manual_seed(seed)
    # torch.cuda.manual_seed(seed)
    # torch.backends.cudnn.daterministic = True
seed_everything(SEED)
```

## データ読み込み
```Python
data = fetch_california_housing()
X = pd.DataFrame(data.data, columns=data.feature_names)
y = pd.Series(data.target)

print(X.shape)
print(y.shape)
```
```
(20640, 8)
(20640,)
```

データを訓練用（検証用含む）と、テスト用に分割する。インデックス番号を振りなおしているのは、K-Foldを効率的に行うため。
```Python
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=SEED)

X_train.reset_index(drop=True, inplace=True)
y_train.reset_index(drop=True, inplace=True)
X_test.reset_index(drop=True, inplace=True)
y_test.reset_index(drop=True, inplace=True)

print(X_train.shape) 
display(X_train.head(3))
print(y_train.shape) 
display(y_train.head(3))
print(X_test.shape) 
display(X_test.head(3))
print(y_test.shape) 
display(y_test.head(3))
```

訓練データが1.6万件、テストデータが0.4万件になった。（データフレームの表示は割愛）
```
(16512, 8)
(16512,)
(4128, 8)
(4128,)
```

## 交差検証の準備
今回は5-Foldとする。
```Python
kf = KFold(n_splits=5, shuffle=True, random_state=SEED)
```

## ハイパーパラメータの最適化
Optunaを使って、ハイパーパラメータを最適化する。最適化する項目と探索範囲はChat-GPTに聞いた…笑。

```Python
categorical_features = [] # カテゴリ変数なし
```

```Python
def objective(trial):
    params = {
        'objective': 'regression',
        'metric': 'rmse',
        'boosting_type': 'gbdt',
        'verbosity': -1,
        'num_leaves': trial.suggest_int('num_leaves', 20, 150),
        'learning_rate': trial.suggest_float('learning_rate', 1e-3, 1e-1, log=True),
        'feature_fraction': trial.suggest_float('feature_fraction', 0.4, 1.0),
        'bagging_fraction': trial.suggest_float('bagging_fraction', 0.4, 1.0),
        'bagging_freq': trial.suggest_int('bagging_freq', 1, 7),
        'min_child_samples': trial.suggest_int('min_child_samples', 5, 100),
        'seed': SEED 
    }

    rmse_scores = []

    # 交差検証を実行
    for tr_index, val_index in kf.split(X_train):
        # 訓練データと検証データに分割
        X_tr = X_train.loc[tr_index]
        X_val = X_train.loc[val_index]
        y_tr = y_train.loc[tr_index]
        y_val = y_train.loc[val_index]

        lgb_tr = lgb.Dataset(X_tr, y_tr, categorical_feature=categorical_features)
        lgb_val = lgb.Dataset(X_val, y_val, reference=lgb_tr, categorical_feature=categorical_features)

        # 訓練
        model = lgb.train(params, 
                          lgb_tr, 
                          valid_sets=[lgb_val], 
                          #num_boost_round=10000,
                          callbacks=[
                              lgb.early_stopping(
                                  stopping_rounds=100,
                                  verbose=False),
                            lgb.log_evaluation(period=100)
                          ]
                          )

        # 予測
        y_pred_val = model.predict(X_val, num_iteration=model.best_iteration)
        rmse = np.sqrt(mean_squared_error(y_val, y_pred_val))
        rmse_scores.append(rmse)
    
    # 交差検証の平均スコアを返す
    return np.mean(rmse_scores)
```

最適化を実行する。

```Python
# 最適化を実行
study = optuna.create_study(direction='minimize')
study.optimize(objective, n_trials=50)
```

最適化したパラメータを表示する。
```Python
print('Best parameters:', study.best_params)
```
```
Best parameters: {'num_leaves': 90, 'learning_rate': 0.07995616486006746, 'feature_fraction': 0.6460671201267157, 'bagging_fraction': 0.7290593706997489, 'bagging_freq': 3, 'min_child_samples': 25}
```

## 最適化したハイパーパラメータで再訓練

最適化したハイパーパラメータで再度訓練する。まずはパラメータを設定する。

```Python
params_tuned = {
    'objective': 'regression',
    'metric': 'rmse',
    'boosting_type': 'gbdt',
    'verbosity': -1,
    'num_leaves': study.best_params['num_leaves'], #trial.suggest_int('num_leaves', 20, 150),
    'learning_rate': study.best_params['learning_rate'], #trial.suggest_float('learning_rate', 1e-3, 1e-1, log=True),
    'feature_fraction': study.best_params['feature_fraction'], #trial.suggest_float('feature_fraction', 0.4, 1.0),
    'bagging_fraction': study.best_params['bagging_fraction'], #trial.suggest_float('bagging_fraction', 0.4, 1.0),
    'bagging_freq': study.best_params['bagging_freq'], #trial.suggest_int('bagging_freq', 1, 7),
    'min_child_samples': study.best_params['min_child_samples'], #trial.suggest_int('min_child_samples', 5, 100),
    'seed': SEED 
}
```

次に、訓練結果を格納するための変数を用意する。
```Python
# 結果格納用変数
y_pred_oof = np.zeros(X_train.shape[0])
y_preds_test = []
rmse_scores = []
models = []
```

訓練実施。
```Python
# 再訓練
for tr_index, val_index in kf.split(X_train, y_train):
    # 訓練データと検証データに分割
    X_tr = X_train.loc[tr_index]
    X_val = X_train.loc[val_index]
    y_tr = y_train.loc[tr_index]
    y_val = y_train.loc[val_index]

    lgb_tr = lgb.Dataset(X_tr, y_tr, categorical_feature=categorical_features)
    lgb_val = lgb.Dataset(X_val, y_val, reference=lgb_tr, categorical_feature=categorical_features)

    # 訓練
    model = lgb.train(params_tuned, 
                      lgb_tr, 
                      valid_sets=[lgb_val], 
                      #num_boost_round=10000,
                      callbacks=[
                        lgb.early_stopping(
                            stopping_rounds=100,
                            verbose=False),
                        lgb.log_evaluation(period=100)
                        ]
                      )

    # 予測
    y_pred_val = model.predict(X_val, num_iteration=model.best_iteration)
    y_pred_test = model.predict(X_test, num_iteration=model.best_iteration)

    # rmseを計算
    rmse = np.sqrt(mean_squared_error(y_val, y_pred_val))

    # 予測値やスコア、学習済みモデルを格納
    y_pred_oof[val_index] = y_pred_val
    y_preds_test.append(y_pred_test)
    rmse_scores.append(rmse)
    models.append(model)
```

## 評価
訓練結果を確認する。まずは各foldでの検証データ（validation)に対するRMSEを表示する。0,435から0.468くらい。前々回のOptuna無しの時は、0.468から0.478くらいだったので、良化したっぽい。
```Python
# rmseを表示
for i in range(len(rmse_scores)):
    print(f'fold{i+1} rmse: {rmse_scores[i]:.4f}')
```
```
fold1 rmse: 0.4681
fold2 rmse: 0.4448
fold3 rmse: 0.4398
fold4 rmse: 0.4590
fold5 rmse: 0.4352
```

平均を取り、CVスコアを算出すると、0.4494となった。前々回は0.4739だったので、平均値も良化。
```Python
cv_score = np.mean(rmse_scores)
print('=== CV score (rmse) ===')
print(f'{cv_score:.4f}')
```
```
=== CV score (rmse) ===
0.4494
```

ハイパーパラメータチューニングの効果が見られたか…。

## テスト
一般的なKaggleコンペなら、テストデータに対する予測を行い、形式を整理し、Submitする。今回はScikit-learn上でデータセットをダウンロードし、最初に訓練データとテストデータに分割したので、テストデータの正解ラベルも手元にあり、RMSEスコアを算出可能。まずは各foldでのテストデータに対する予測値を表示する。

```Python
# 各foldのテストデータに対する予測値を表示
for i in range(len(y_preds_test)):
    print(f'fold{i+1}: {y_preds_test[i][:5]} ...')
```
```
fold1: [0.48554674 1.24362198 5.08410152 2.3051161  2.3623417 ] ...
fold2: [0.52486522 0.91217026 5.1110846  2.39300856 2.4203427 ] ...
fold3: [0.47965315 0.83645024 4.94657118 2.32503299 2.34725133] ...
fold4: [0.50302413 0.94841206 5.1756785  2.45468694 2.33222264] ...
fold5: [0.4919892  0.8870616  5.06584356 2.35460799 2.26655994] ...
```

最初の5個しか表示していないが、変な値は出ていなさそう。平均を取って、最終的な予測値とする。この辺りはいろいろやり方があると思うが、書籍等を参考にして、このやり方を採用する。軽くアンサンブル学習になっているがよいところか？

```Python
y_sub = sum(y_preds_test) / len(y_preds_test)
print(f'{y_sub[:5]} ...')
```
```
[0.49701569 0.96554323 5.07665587 2.36649051 2.34574366] ...
```

目視で見た感じ、計算間違いはしていなさそう。今回、提出はしないが、csvファイルを作成する。
```Python
# 提出用ファイルを作成
sub = pd.DataFrame({'target': y_sub})
sub.to_csv('submission.csv', index=False)
display(sub.head())
```

今回、テストデータの正解ラベルは手元にあるので、採点する。
```Python
# テストデータの予測値を使って採点
rmse_sub = np.sqrt(mean_squared_error(y_test, y_sub))
print('=== RMSE for test data ===')
print(f'{rmse_sub:.5f}')
```
```
=== RMSE for test data ===
0.44027
```
結果は、0.4402。CVスコアが0.4494なので、過学習はしていなさそう。むしろアンサンブルの効果が出てる？前々回（Optuna無し）と比べて、0.4622から0.4402とスコアアップしており、ハイパーパラメータチューニングの効果が見て取れた。よかった。

# おわりに
california_housingデータセットを使って回帰問題を解いた。Optunaを導入したところ、スコアアップした。効果はわずかだが、実際のコンペでもこれくらいのイメージ。上位を狙うなら、チューニングは欠かせないか…。特徴量エンジニアリングと違って、基本的にやれば効果につながるのがハイパーパラメータチューニングのいいところか（過学習は別の話…）。