はじめに

今回私は最近はやりのchatGPTに興味を持ち、深層学習について学んでみたいと思い立ちました！
深層学習といえばPythonということなので、最終的にはPythonを使って深層学習ができるとこまでコツコツと学習していくことにしました。
ただ、勉強するだけではなく少しでもアウトプットをしようということで、備忘録として学習した内容をまとめていこうと思います。
この記事が少しでも誰かの糧になることを願っております！
※投稿主の環境はWindowsなのでMacの方は多少違う部分が出てくると思いますが、ご了承ください。
最初の記事：Python初心者の備忘録 #01
前の記事：Python初心者の備忘録 #23 ~機械学習入門編09~
次の記事：まだ

今回は特徴量エンジニアリング、特徴量選択、ハイパーパラメータチューニングについてまとめております。

■学習に使用している資料

Udemy：【本番編】米国データサイエンティストがやさしく教える機械学習超入門【Pythonで実践】

■特徴量エンジニアリング

▶特徴量エンジニアリングとは

すでにある特徴量から予測に有益な特徴量を作り出すことで非常に重要な処理となる
（他の前処理や特徴量選択などを行う時間を削ってでも注力した方が良いほど）
※なるべく多くの特徴量を作成し、試すのがいい
- 日付情報から年、月、日、曜日...etc
- 多項式特徴量($X^2_1, X_1X_2$...etc)
- binning(数値 -> カテゴリカル)
- 四則演算(例：利益 = 売上 - 経費)
- 多テーブルからの積み上げ値
- ドメイン的に意味のある特徴量
  ...etc

▶日付情報

最もよく使われる特徴量エンジニアリングの1つで時系列データに限らず有効
データの性質によって注意が必要なものがある
- 例①：テストデータが学習データの後続データである場合はYear情報はそのまま使えない
- 例②：2年未満のデータの場合、月の周期的な季節性をとらえられない
  ※時、分、秒などはデータの粒度に応じて分析/処理する

Pythonで日付情報による特徴量エンジニアリング

次で使用するデータは、下記から取得してください。
https://github.com/Yushin-Tati/Learning_machine_learning/tree/main
./bike_share.csv

カラムの説明
datetime - hourly date + timestamp
season - 1 = spring, 2 = summer, 3 = fall, 4 = winter holiday - whether the day is considered a holiday
workingday - whether the day is neither a weekend nor holiday
weather -
1: Clear, Few clouds, Partly cloudy, Partly cloudy
2: Mist + Cloudy, Mist + Broken clouds, Mist + Few clouds, Mist
3: Light Snow, Light Rain + Thunderstorm + Scattered clouds, Light Rain + Scattered clouds
4: Heavy Rain + Ice Pallets + Thunderstorm + Mist, Snow + Fog
temp - temperature in Celsius
atemp - "feels like" temperature in Celsius
humidity - relative humidity
windspeed - wind speed
casual - number of non-registered user rentals initiated
registered - number of registered user rentals initiated
count - number of total rentals

日付カラムの確認
df.info()
日付形式に変換
pandas.to_datetime()
特徴量生成
- month：.dt.month
- year：.dt.year
- day of year：.dt.dayofyear
- day of week：.dt.dayofweek
- week of year：.dt.isocalendar().week
- quarter：.dt.quarter
- leap year：.dt.is_leap_year

date feature engineering

import pandas as pd
import seaborn as sns
from sklearn.preprocessing import PolynomialFeatures

# データ準備
df = pd.read_csv('bike_share.csv')

# 日付カラムの確認
df.info()

# datetimeカラムを日付のデータタイプに変換
df['datetime'] = pd.to_datetime(df['datetime'])

# datetimeカラムのDtypeがdatetimeになっていることを確認
df.info()

# 日付データからそれぞれの値を取得
df['year'] = df['datetime'].dt.year # 年
df['month'] = df['datetime'].dt.month # 月
df['dayofyear'] = df['datetime'].dt.dayofyear # その年の何日目か
df['dayofweek'] = df['datetime'].dt.dayofweek # 曜日(月曜スタート)
df['weekofyear'] = df['datetime'].dt.isocalendar().week # その年の何週目か
df['quarter'] = df['datetime'].dt.quarter # その日時がどの四半期に該当するかのリスト
df['isleap'] = df['datetime'].dt.is_leap_year # うるう年かどうか

▶多項式特徴量(Polynomial Features)

[復習]
下記記事参照
▶多項式特徴量

補足：決定木と交互作用項

決定木は分岐を繰り返すことで交互作用項を表現できる
ただ、決定木(GBDT含む)を使い場合は無理にいろいろな種類の交互作用項を作る必要はない

Pythonで多項式特徴量による特徴量エンジニアリング

Polynomial Features

# データロード
df = pd.read_csv('penguins_size.csv')
df.dropna(inplace=True)

# 多項式特徴量作成
poly = PolynomialFeatures(degree=2, include_bias=False)
df[['culmen_length_mm', 'culmen_depth_mm','culmen_length_mm^2', 'culmen_length_mm * culmen_depth_mm', 'culmen_depth_mm^2' ]] = poly.fit_transform(df[['culmen_length_mm', 'culmen_depth_mm']])
sns.pairplot(df, hue='species')

▶binning

数値の特徴量からカテゴリカルの特徴量を生成する方法で、元のカラムと組み合わせて使用してもいい
※実際には全データの範囲に対してパーセンタイルによってbinnnigをすることが多い

Pythonでbinnnigによる特徴量エンジニアリング

pandas.cut(df['カラム名'], bins, labels)で実行可能
[引数]
- bins：binの数を指定
- labels：変換後のカテゴリーのリスト(Falseの場合はbin番号を返す)

binnnig

# binnnig
df['body_mass_g_bin'] = pd.cut(df['body_mass_g'], bins=10, labels=False)

# カテゴリカル変数として扱う場合dtypeをobjectに変換。(数値変数として扱ってもよい)
df['body_mass_g_bin'] = df['body_mass_g_bin'].astype('object')
df.info()

▶四則演算

特徴量同士を四則演算して新たな特徴量を生成する方法
網羅的にやるのは大変なので、ドメイン的に意味のある特徴量を生成するように考える
集約値との四則演算も有効
カテゴリカル変数同士の結合にも有効
- 例：'good' + 'east' = 'good_east'
- target encodinでより目的変数とカテゴリ地の組み合わせを細分化できるが、過学習には注意する
  (※カテゴリに対してのレコード数が減るため)

Pythonで四則演算による特徴量エンジニアリング

Series同士の演算を行えばよい
- df['new_feature'] = df['feature1'] + df['feature2']

四則演算でfeature engineering

# 今まで加工してきたペンギンデータをそのまま使用
df['culmen_diff'] = df['culmen_length_mm'] - df['culmen_depth_mm']
df['culmen_ratio'] = df['culmen_length_mm'] / df['culmen_depth_mm']
%matplotlib inline
sns.pairplot(df, hue='species')

▶集約値

1:Nの関係にある2つのテーブルがある場合に使用できる方法で、テーブル結合の際に集約値を使用する
(※min, max, sum, mean, median, std...etc)

Pythonで集約値による特徴量エンジニアリング

次にで使用しているデータは下記から取得してください。
https://github.com/Yushin-Tati/Learning_machine_learning
./feature engineering with aggregated valus
・homecredit_application.csv
・homecredit_previous_application.csv
(・homecredit_columns_description.csv) ← 上記データのカラム説明

SK_ID_CURRをキーとして過去の累積借用額(credit)の合計を新たな特徴量とする

aggregated values

# データロード
df_app = pd.read_csv('homecredit_application.csv')
df_prev_app = pd.read_csv('homecredit_previous_application.csv')

# 見やすいようにソート
df_prev_app.sort_values('SK_ID_CURR')

# 集約値計算
credit_prev_sum = df_prev_app.groupby('SK_ID_CURR').sum()[['AMT_CREDIT']]

# SK_ID_CURRをキーに表を結合
df_app.join(credit_prev_sum, on='SK_ID_CURR', rsuffix='_PREV_SUM')

▶submodelを使った予測値

本来の目的変数以外の特徴量をtargetにしてモデルを学習する方法
※欠損値補完以外に特徴量エンジニアリングにも使うことができる

■特徴量選択

▶特徴量選択とは

特徴量は沢山あればいいというものではなく(次元の呪い)、本当に必要なものを選別する必要がある

特徴量選択の手法

greedy feature selection
recursive feature elimination
モデルの特徴量の重要度で判断
- 線形回帰等のおあらメトリックモデルの係数を見る(t検定量など)
- 決定木などの特徴量重要度を見る
L1正則化項による特徴量選択

▶greedy feature selection

greedy(貪欲)に特徴量の追加と評価を繰り返す手法
1. それぞれの特徴量単体でモデルを評価し、最高精度の特徴量を選択する
2. 残りの特徴量を1つ追加してモデルを評価し、残りの特徴量全てに実施して最高精度の特徴量を選択する
3. 2.を最高精度を更新しなくなるまで繰り返す
計算量が高くなりやすいので注意が必要だが、高精度なモデルが必要な際に有効

[例]
1. それぞれの特徴量単体でモデルを評価し、最高精度の特徴量を選択する

2. 残りの特徴量を1つ追加してモデルを評価し、残りの特徴量全てに実施して最高精度の特徴量を選択する

3. 最高精度を更新しなくなるまで繰り返す

最終的にAとCが選択される

Pythonでgreedy feature selection

前処理
- 欠損値が多いデータは削除
- カテゴリカルの欠損値にはNaNを代入して新しいカテゴリとして扱う
- one-hot encodingでダミー変数化
- 標準化
特徴量エンジニアリング
- 多項式特徴量
- 四則演算
greedy feature selectionをクラスで実装
- k-Foldで評価
- 引数にPipelineを受け取るようにする
制度の推移と、選択された特徴量を確認

greedy feature selection

import pandas as pd
import numpy as np
from sklearn.impute import SimpleImputer
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.preprocessing import StandardScaler, PolynomialFeatures
from sklearn.pipeline import Pipeline
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import KFold, cross_val_score
from sklearn.feature_selection import RFE, RFECV, SelectFromModel
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier

# データ準備
df = pd.read_csv('penguins_size.csv')
df.loc[df[df['sex']=='.'].index, 'sex'] = np.nan # sexカラムの'.'をNaNに置き換え

# 欠損値確認
df.info()

★#2～6に欠損値がありそう

# 欠損が多いデータは落とす
df.dropna(thresh=3, inplace=True)

# X, y 作成
X = df.drop('species', axis=1)
y = df['species']

# 欠損値代入
imputer = SimpleImputer(strategy='constant', fill_value='NaN')
imputer.set_output(transform='pandas')
# X全体にfit_transformをすると，全てのカラムのDtypeがobjectになってしまい，後続のダミー変数作成の際に全カラムに対してOHEしてしまうため，
# ここではsexカラムを指定して，sexカラムにのみ実施する
X['sex'] = imputer.fit_transform(X[['sex']])

# OHE
class GetDummies(BaseEstimator, TransformerMixin):
    
    def __init__(self):
        self.columns = None
        
    def fit(self, X, y=None):
        self.columns = pd.get_dummies(X).columns
        return self
    
    def transform(self, X):
        X_new = pd.get_dummies(X)
        return X_new.reindex(columns=self.columns, fill_value=0)

# 特徴量エンジニアリング
# 多項式特徴量
poly = PolynomialFeatures(degree=2, include_bias=False)
X[['culmen_length_mm', 'culmen_depth_mm', 'culmen_length_mm^2', 'culmen_length_mm * culmen_depth_mm', 'culmen_depth_mm^2']] = poly.fit_transform(X[['culmen_length_mm', 'culmen_depth_mm']])

# 四則演算
X['culmen_diff'] = X['culmen_length_mm'] - X['culmen_depth_mm']
X['culmen_ratio'] = X['culmen_length_mm'] / X['culmen_depth_mm']

# Pipeline
pipeline = Pipeline([('dummy', GetDummies()),
          ('scaler', StandardScaler()),
          ('model', LogisticRegression())])

# cv
cv = KFold(n_splits=3, random_state=0, shuffle=True)

# クラス作成
class GreedyFeatureSelection():
    
    def __init__(self, pipeline, cv):
        self.pipeline = pipeline
        self.cv = cv
        self.selected_features = []
        self.scores = [0]
    
    def select_feature(self, X, y):
        
        all_features = X.columns
        
        while True:
            # print('greedy selection started')
            best_score = self.scores[-1]
            candidate_feature = None
            for feature in all_features:
                if feature in self.selected_features:
                    continue
                # print(f'{feature} started')
                features = self.selected_features + [feature]
                X_train = X[features]
                # 評価
                score = cross_val_score(self.pipeline, X_train, y, cv=self.cv).mean()
                # print(f'{features} score: {score}')
                if score > best_score:
                    # print(f'best score updated {best_score} -> {score}')
                    best_score = score
                    candidate_feature = feature
            
            if candidate_feature is not None:
                # print(f'========{candidate_feature} is selected=============')
                self.scores.append(best_score)
                self.selected_features.append(candidate_feature)
            else:
                break


# Greedy feature selection
gfs = GreedyFeatureSelection(pipeline=pipeline, cv=cv)
gfs.select_feature(X, y)

# スコアの結果と選択された特徴量を確認
print(gfs.scores)  # -> [0, 0.9619883040935672, 0.9912280701754387, 0.9970760233918128, 1.0]
print(gfs.selected_features)  # -> ['culmen_ratio', 'island', 'culmen_diff', 'body_mass_g']

▶recursive feature elimination(RFE)

greedy feature selectionの逆で、すべての特徴量でモデルを構築してこら特徴量を1つずつ減らしていく手法
1. 全ての特徴量でモデルを学習
2. モデルへの影響が最も低い特徴量を落とし、残りの特徴量でモデルを学習する
3. 2.を指定した特徴量の数になるまで実施する
  → あらかじめ最終的に残す特徴量の数を決めておく必要がある
モデルへの影響
- 線形回帰の係数(p値)や決定木の重要度を使用する
- 重要度などの指標を持ったモデルを使用する

PythonでRFE

sklearn.feature_selection.RFEクラスを使用
1. インスタンス生成
  - estimator：sklearnのモデルインスタンス
  - n_feature_to_select：最終的な特徴量数
2. .fit(X, y)でRFEを実施し、特徴量選択
3. .transform(X)で削減された特徴量X'を返す
  ※.ranking_で特徴量の重要度順を確認できる
  ※.support_で特徴量の選択結果のマスクを取得できる
sklearn.feature_selection.RFECVクラスでRFEをCVで実行可能

RFE

# データはgreedy feature selectionの時のものをそのまま流用
# モデルは決定木、残す特徴量数は6個
rfe = RFE(DecisionTreeClassifier(), n_features_to_select=6)
X = pd.get_dummies(X, drop_first=True)
rfe.fit(X, y)
rfe.transform(X)
"""結果
array([[20.4       ,  2.09090909,  0.        ,  1.        ,  1.        ,
         0.        ],
       [22.1       ,  2.27011494,  0.        ,  1.        ,  0.        ,
         0.        ],
       [22.3       ,  2.23888889,  0.        ,  1.        ,  0.        ,
         0.        ],
       …,
       [34.7       ,  3.21019108,  0.        ,  0.        ,  1.        ,
         0.        ],
       [30.4       ,  3.05405405,  0.        ,  0.        ,  0.        ,
         0.        ],
       [33.8       ,  3.09937888,  0.        ,  0.        ,  1.        ,
         0.        ]])
"""

# 特徴量の重要度ランキング
rfe.ranking_  # -> array([8, 7, 6, 5, 4, 3, 2, 1, 1, 1, 1, 1, 1])

# 選択された特徴量のリストを取得
# rfe.feature_names_in_[rfe.ranking_==1]
rfe.feature_names_in_[rfe.support_]  # -> array(['culmen_diff', 'culmen_ratio', 'island_Dream', 'island_Torgersen', 'sex_MALE', 'sex_NaN'], dtype=object)

# 選択した特徴量をPipelineに組み込む
pipeline = Pipeline([('rfe', rfe), ('scaler', StandardScaler()), ('model', LogisticRegression())])
scores = cross_val_score(pipeline, X, y, cv=cv)
print(scores)  # -> [1.         1.         0.99122807]

# RFECV
rfecv = RFECV(DecisionTreeClassifier(), cv=cv)
rfecv.fit(X, y)
rfecv.transform(X)
"""結果
array([[20.4       ,  2.09090909,  0.        ],
       [22.1       ,  2.27011494,  0.        ],
       [22.3       ,  2.23888889,  0.        ],
       …,
       [34.7       ,  3.21019108,  0.        ],
       [30.4       ,  3.05405405,  0.        ],
       [33.8       ,  3.09937888,  0.        ]])
"""

rfecv.feature_names_in_[rfecv.support_]  # -> array(['culmen_diff', 'culmen_ratio', 'island_Dream'], dtype=object)

# 特徴量の数と精度の推移を描画
import matplotlib.pyplot as plt
plt.errorbar(range(1, 1+len(rfecv.cv_results_['mean_test_score'])),
            rfecv.cv_results_['mean_test_score'],
            yerr=rfecv.cv_results_['std_test_score'])

★1～3までは精度が向上しているが、それ以上は変化がないので、重要度3位までの特徴量が選択されている妥当性が見てとれる

■モデルの係数や重要度での特徴量選択

モデルが出力した特徴量の重要度を使って閾値を定めて選択する
- 線形回帰やロジスティック回帰の係数
- 決定木の重要度...etc
異なるモデルで重量度を確認することも可能
- 例：ロジスティック回帰で係数をもとに特徴量選択をし、決定木で予測モデル構築
気軽に利用できるが、閾値を決定するのが難しく予算内で将来取得する特徴量を決めたいときなどに有効
※精度が極端に落ちる可能性もあるので注意

復習

線形回帰の係数の統計的仮設検定
■線形回帰の解釈
決定木の解釈
■決定木

Pythonでモデルの係数や重要度で特徴量選択

sklearn.feature_selection.SelectFormModelクラスを使用する
1. インスタンス生成
  - estimator：sklearnのモデルインスタンス
  - threshold：重要度の閾値、meanやmedian、1.25*mean等も可能
2. .fit(X, y)で学習
3. .transform(X)で削減された特徴量X'を返す
  ※.estimator_でモデルインスタンスにアクセス可能
  ※.get_support()で特徴量の選択結果のマスクを取得できる

# RandomForestを使って特徴量の重要度が低い特徴量を落とす
sfm = SelectFromModel(RandomForestClassifier(random_state=0))
X_selected = sfm.fit_transform(X, y)
"""結果
array([[ 39.1       ,  18.7       , 181.        , ..., 349.69      ,
         20.4       ,   2.09090909],
       [ 39.5       ,  17.4       , 186.        , ..., 302.76      ,
         22.1       ,   2.27011494],
       [ 40.3       ,  18.        , 195.        , ..., 324.        ,
         22.3       ,   2.23888889],
       …,
       [ 50.4       ,  15.7       , 222.        , ..., 246.49      ,
         34.7       ,   3.21019108],
       [ 45.2       ,  14.8       , 212.        , ..., 219.04      ,
         30.4       ,   3.05405405],
       [ 49.9       ,  16.1       , 213.        , ..., 259.21      ,
         33.8       ,   3.09937888]])
"""

# .estimator_で特徴量選択に使用したモデルにアクセス
sfm.estimator_.feature_importances_  # -> array([0.08924232, 0.07961904, 0.09933245, 0.0525517 , 0.08072858, 0.01467092, 0.10925035, 0.14543881, 0.24562145, 0.07424349, 0.00630753, 0.00299336, 0.        ])

# 選択された特徴量のリストを取得
sfm.feature_names_in_[sfm.get_support()]  # -> array(['culmen_length_mm', 'culmen_depth_mm', 'flipper_length_mm', 'culmen_length_mm^2', 'culmen_depth_mm^2', 'culmen_diff', 'culmen_ratio'], dtype=object)

▶正則化項(L1ノルム：Lasso)による特徴量選択

復習

正則化項
▶Lasso(L1ノルム)

PythonでL1ノルムによる特徴量選択

sklearn.feature_selecDon.SelectFromModelクラスを使用
1. インスタンス生成
  - esDmator：L1正則化項のモデルインスタンスを渡す
2. .fit(X, y)で学習
3. .transform(X)で削減された特徴量X’を返す
  ※.esDmator_でモデルインスタンスにアクセス可能
  ※.get_support()で特徴量の選択結果のマスクを取得

Lasso feature engineering

# estimtatorにl1正則を使用するモデルを入れると，l1正則を使って特徴量選択が可能
sfm = SelectFromModel(LogisticRegression(penalty="l1", solver='liblinear'))
pipeline = Pipeline([('scaler', StandardScaler()),
                     ('feature selection', sfm)])
pipeline.set_output(transform='pandas')
pipeline.fit_transform(X, y)

# 使用したロジスティック回帰の係数(3クラス分)を確認
sfm.estimator_.coef_
"""結果
array([[ 0.        ,  0.        , -0.17837796,  0.        ,  0.        ,
         0.        ,  0.08218938, -6.97592113,  0.        , -0.34900286,
         0.88982087,  1.17944793,  0.        ],
       [ 4.21255667,  0.        ,  0.        , -2.36976069,  0.        ,
         0.14737444,  0.        ,  0.        ,  0.        ,  2.66136579,
         0.        , -0.75593737,  0.        ],
       [ 0.        , -2.75559072,  1.69766816,  1.80788218,  0.        ,
         0.        ,  0.        ,  0.        ,  0.5335975 , -0.25339517,
         0.        ,  0.        ,  0.        ]])
"""

# 選択された特徴量のリスト
sfm.feature_names_in_[sfm.get_support()]  # -> array(['culmen_length_mm', 'culmen_depth_mm', 'flipper_length_mm', 'body_mass_g', 'culmen_length_mm * culmen_depth_mm', 'culmen_depth_mm^2', 'culmen_diff', 'culmen_ratio', 'island_Dream', 'island_Torgersen', 'sex_MALE'], dtype=object)

■ハイパーパラメータチューニング

▶ハイパーパラメータチューニングとは

それぞれのモデルには様々なハイパーパラメータがあり、それを最適化することを指す
例：学習率$\alpha$、バッチサイズ、エポック数...etc
チューニング方法
- Grid Search：総当たりで最適な組み合わせを探す
- Random Search：ランダムな組み合わせで探す
- Optimization：関数を定義してそれを最小化して探す

CatBoost、LightBoost、XGBoostのハイパーパラメータを確認したい場合は下記を見てみてください。

▶Grid Search

あらかじめ指定した範囲でハイパーパラメータの組み合わせを総当たり
で試す方法で、事前にある程度範囲を指定する必要がある
→ 公式ドキュメントや他の事例を参考にする
計算量が非常に高くなることに注意する、また必ず乱数シード(random seed)を固定する

モデルのハイパーパラメータの意味を理解することで，無駄な組み合わ
せを減らすことができる

PythonでGrid Search

sklearn.model_selec7on.GridSearchCVクラスを使用
1. インスタンス生成
  - estimator：sklearnのモデルインスタンス
  - param_grid：探索するパラメータと値のdictionary {””: [value, ..], …}
  - scoring：評価指標
  - cv：sklearnのcvオブジェクト
2. .fit(X, y)で学習
  - **fit_paramsにestimatorのfitメソッドの引数を渡すことが可能
3. .best_params_や.best_score_で結果を確認

Grid Search

import pandas as pd
import numpy as np
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split, GridSearchCV, KFold, RandomizedSearchCV, cross_val_score
from sklearn.impute import SimpleImputer
from sklearn.compose import ColumnTransformer
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.pipeline import Pipeline
from scipy.stats import uniform 
from scipy.stats import randint as sp_randint

import lightgbm as lgb
from hyperopt import hp, tpe
from hyperopt.fmin import fmin
import plotly.express as px

# データロード
dataset = fetch_california_housing()

# データの説明を確認
print(dataset['DESCR'])

# データ準備
X = pd.DataFrame(dataset['data'], columns=dataset['feature_names'])
y = dataset['target']

# hold-out
X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=0)

# モデル
lgbmr = lgb.LGBMRegressor(random_state=0)
param_grid = {'num_leaves': [10, 20, 30, 40, 50, 60],
              'max_depth': [5, 10, 15], 
              'reg_alpha': [0, 0.01, 0.03]}

cv = KFold(n_splits=3, random_state=0, shuffle=True)
gs = GridSearchCV(lgbmr, param_grid=param_grid, cv=cv)

# early_stopping
callbacks = [lgb.early_stopping(stopping_rounds=10)]
eval_set = [(X_val, y_val)]
fit_params = {'callbacks': callbacks, 'eval_set': eval_set}
gs.fit(X_train, y_train, **fit_params)

# grid searchの結果確認
cv_results_df = pd.DataFrame(gs.cv_results_)
cv_results_df.head(3)

※「rank_test_score」が全ての組み合わせの内、何番目に結果が良かったかを表しているので、これが1の組み合わせが最適なチューニングだと分かる

# top5の結果を表示
for index, row in cv_results_df[cv_results_df['rank_test_score']<=5].iterrows():
    print(f'{row["rank_test_score"]}: {row["params"]}')
    print(f'{row["mean_test_score"]}')

"""結果
4: {'max_depth': 10, 'num_leaves': 50, 'reg_alpha': 0.01}
0.8329725903236045
3: {'max_depth': 15, 'num_leaves': 50, 'reg_alpha': 0}
0.8331058990226472
2: {'max_depth': 15, 'num_leaves': 50, 'reg_alpha': 0.03}
0.8332059773155042
5: {'max_depth': 15, 'num_leaves': 60, 'reg_alpha': 0}
0.8327661059213919
1: {'max_depth': 15, 'num_leaves': 60, 'reg_alpha': 0.01}
0.8338051795879502
"""

# 最も良いハイパーパラメータの組み合わせで再度学習し評価
lgbmr = lgb.LGBMRegressor(**gs.best_params_, random_state=0, learning_rate=0.01, n_estimators=1000)
callbacks = [lgb.early_stopping(stopping_rounds=10)]
eval_set = [(X_val, y_val)]
fit_params = {'callbacks': callbacks, 'eval_set': eval_set}
lgbmr.fit(X_train, y_train, **fit_params)

# チューニング後の学習の評価結果
lgbmr.score(X_val, y_val)  # -> 0.8499232425233268

PipelineにGrid Searchを組み込む

# データロード
df = pd.read_csv('penguins_size.csv')

# データクリーニング
df.loc[df[df['sex']=='.'].index[0], 'sex'] = np.nan
df.dropna(thresh=3, inplace=True)

# X, y作成
X = df.drop('species', axis=1)
y = df['species']

# 欠損値代入
cat_cols = X.select_dtypes(exclude=np.number).columns.to_list()
imputer = SimpleImputer(strategy='most_frequent')
ct = ColumnTransformer([('cat_imputer', imputer, cat_cols)], remainder='passthrough')
ct.set_output(transform='pandas')

# OHE
class GetDummies(BaseEstimator, TransformerMixin):
    
    def __init__(self):
        self.columns = None
        
    def fit(self, X, y=None):
        self.columns = pd.get_dummies(X).columns
        return self
    
    def transform(self, X):
        X_new = pd.get_dummies(X)
        return X_new.reindex(columns=self.columns, fill_value=0)

# モデル
lgbmc = lgb.LGBMClassifier(random_state=0)
    
# pipeline
pipeline = Pipeline([('impute',ct), ('dummy', GetDummies()), ('model', lgbmc)])

# cv
cv = KFold(n_splits=3, random_state=0, shuffle=True)

# grid search
# pipelineを使う場合はキーの接頭辞にtransformerの名前を__で繋げてつけることに注意
param_grid = {'model__num_leaves': [10, 20, 30, 40, 50, 60],
              'model__max_depth': [5, 10, 15], 
              'model__reg_alpha': [0, 0.01, 0.03]}
gs = GridSearchCV(pipeline, param_grid=param_grid, cv=cv)
gs.fit(X, y)

# 結果をDFに変換
cv_results_df = pd.DataFrame(gs.cv_results_)

▶Random Search

ランダムにハイパーパラメータを組み合わせて探索する方法
ハイパーパラメータが多いモデルはGrid Searchより現実的
※必ず乱数シード(random seed)を固定する

PythonでRandom Search

sklearn.model_selec7on.RandomizedSearchCVクラスを使用する
1. インスタンス生成
  - es7mator： sklearnのモデルインスタンス
  - param_distribu7ons：探索するパラメータと値を出力するobject(または値のリスト)のdict(例: {”name”: object, …})
    ※objectは.rvs()メソッドを実装している必要があり、scipy.stats.uniformやscipy.stats.randintを使用することが多い
  - scoring：評価指標
  - cv： sklearnのcvオブジェクト
  - n_iter：イテレーションの回数
2. .fit(X, y)で学習
3. .best_params_や.best_score_で結果を確認

Random Search

# データ準備
dataset = fetch_california_housing()
X = pd.DataFrame(dataset['data'], columns=dataset['feature_names'])
y = dataset['target']
X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=0)

# モデル
lgbmr = lgb.LGBMRegressor(random_state=0)
param_dist = {'num_leaves': sp_randint(10, 60),
              'max_depth': sp_randint(5, 15), 
              'reg_alpha': uniform(0, 0.03)}

# early stopping
eval_set = [(X_val, y_val)]
callbacks = [lgb.early_stopping(stopping_rounds=10)]
fit_params = {'callbacks': callbacks, 'eval_set': eval_set}

# cv
cv = KFold(n_splits=3, shuffle=True, random_state=0)

# random search
rs = RandomizedSearchCV(lgbmr, param_distributions=param_dist, cv=cv, n_iter=36)
rs.fit(X_train, y_train, **fit_params)

# random searchの結果確認
cv_results_df = pd.DataFrame(rs.cv_results_)
cv_results_df.head(3)

# top5を表示
for _, row in cv_results_df[cv_results_df['rank_test_score']<=5].iterrows():
    print(f"{row['rank_test_score']}: {row['params']}")
    print(f"{row['mean_test_score']}")

"""結果
5: {'max_depth': 12, 'num_leaves': 45, 'reg_alpha': 0.024576442334105927}
0.8310624277355725
4: {'max_depth': 12, 'num_leaves': 45, 'reg_alpha': 0.006443170680683113}
0.8314624085008285
3: {'max_depth': 13, 'num_leaves': 49, 'reg_alpha': 0.005793462333491382}
0.8321267477795421
1: {'max_depth': 14, 'num_leaves': 59, 'reg_alpha': 0.019163009094849987}
0.8334056557875065
2: {'max_depth': 14, 'num_leaves': 58, 'reg_alpha': 0.010214275208938901}
0.8333039253544552
"""
# Grid Search -> 1: {'max_depth': 15, 'num_leaves': 60, 'reg_alpha': 0.01} 0.8338051795879502
# 今回はGrid Searchのほうが若干精度がいいが、人が決めた数値以外の検討もしてくれるので、Random Searchのほうが好ましい場合が多い

# 最適を取得
rs.best_params_  # -> {'max_depth': 14, 'num_leaves': 59, 'reg_alpha': 0.019163009094849987}

▶最適化手法を用いたチューニング

最適化アルゴリズムを用いて最適なハイパーパラメータを探索する
Grid SearchやRandom Searchでは，過去に試した結果を反映できない

ベイズ最適化

ハイパーパラメータ最適化で最も人気な手法の一つ
ベイズ最適化はブラックボックス最適化問題（真の関数が不明である場合の最適化問題）の一つ
- 今回であれば「ハイパーパラメータを入力，損失を出力とする関数は不明」という状況であり、勾配情報が使えない
- そこで、真の関数を推定した関数(獲得関数: acquisition function)を使って最適化問題を解く
  → 真の関数を最小にする入力を探す

-> Expected Improvementが高くなるように最適化を進める

Expected Improvement(EI)

評価値の改善量 (EI: Expected Improvement)
EIの値が最も高くなるxの評価値が高くなると考える

▶ベイズの定理

※先に事前確率と事後確率と条件付確率について説明する。

事前確率と事後確率と条件付き確率
$P(A)$：事象Aが起こる確率
$P(A|B)$：事象Bが起きている場合に事象Aが起こる確率

事前/事後確率と同時確率
$P(A, B)$：事象Aと事象(B)が同時に起こる確率 = $P(A)P(B)$

ベイズの定理：$P(A|B)=\frac{P(B|A)P(A)}{P(B)}$

Bが起きた時にAが起きる確率は，Aが起きた時にBが起きる確率を用いて
表すことができる
$P(A|B)$を求めるのは難しくても， $P(B|A)$は求めやすい場合が多々ある

ベイズの定理をIEの式に適用する

→ yが与えられた時のxの確率密度の考え方を次で説明する。

▶Tree-Structured Parzen Estimator (TPE)

$P(x|y)$：$y$が与えられた時の$x$の確率密度
-> 過去の結果$y$が得られた時のハイパーパラメータ$x$の確率密度
ベースライン$y^*$より上か下かで場合分けを行い、それぞれのケースで異なる分布を使用する(※上なら改善がない、下なら改善されたと判断する)

過去に得られた$x$と$y$からParzen Estimator(カーネル密度推定)を使って関数を推定する

Pythonでベイズ最適化によるハイパーパラメータチューニング

hyperoptライブラリを使用する
hyperopt.fmin.fmin()
- objective：最小化したい目的関数(pythonの関数)を渡す
  1. モデルインスタンス生成
  2. 学習および評価
  3. 精度を返す
- space：変数の範囲をdictonary形式で渡す
  -> hyperopt.hp.quniform()や hyperopt.hp.uniform()を使用
- algo： hyperopt.tpe.suggest
- max_evals：イテレーション回数

ベイズ最適化

# データの準備
dataset = fetch_california_housing()
X = pd.DataFrame(dataset['data'], columns=dataset['feature_names'])
y = dataset['target']

# 関数
def objective(params):
    
    # paramsに入っている値は全てfloatなので，intにcastする
    params = {'num_leaves': int(params['num_leaves']),
             'max_depth': int(params['max_depth']), 
             'reg_alpha': params['reg_alpha']}
    model = lgb.LGBMRegressor(**params, random_state=0)
    cv = KFold(n_splits=3, shuffle=True, random_state=0)
    scores = cross_val_score(model, X, y, cv=cv)
    
    # hyperoptにはlogを取る仕組みがないため別途logを保存させる
    log['params'].append(params)
    log['score'].append(scores.mean())
    log['score_std'].append(scores.std())
    
    # 最小化をすることを目指すので，scoreは負の形にして低ければ精度が良いという形にする
    return -scores.mean()

# 検証範囲を指定
space = {'num_leaves': hp.quniform('num_leaves', 10, 60, 2),
         'max_depth': hp.quniform('max_depth', 5, 15, 2), 
         'reg_alpha': hp.uniform('reg_alpha', 0, 0.03)}

log = {'params': [], 'score': [], 'score_std': []}

# ベイズ最適化
best = fmin(objective, space=space, algo=tpe.suggest, max_evals=100)

# logを確認
log_df = pd.DataFrame(log)
log_df = log_df.sort_values('score', ascending=False)
log_df[['num_leaves', 'max_depth', 'reg_alpha']] = log_df.apply(lambda row: pd.Series(row['params']), axis=1)
log_df.head(3)

# ベイズ最適化により選択された100個のハイパーパラメータの組み合わせを描画
# スコアにより色を変えて，スコアが高い付近に組み合わせが集中していることが確認できる
fig = px.scatter_3d(log_df, x='num_leaves', y='max_depth', z='reg_alpha', color='score')
fig.show()

次の記事

まだ