MLエンジニアへの道 #12

このブログ記事では、PythonにおけるXGBoostの実装方法を紹介します。

注意: XGBoostを本格的な用途で実装しようとしている場合は、公式のXGBoostのドキュメンテーションをよく理解してから行うようにしてください。本記事はXGBoostの Pythonパッケージのすべての側面をカバーするものではなく、初めてXGBoostを実装しようとしているカジュアルな読者向けに書かれています。

ステップ1 & 2. データの探索 & 前処理

この記事では、新しいデータセット「California Housing Price回帰データセット」を使用します。このデータセットは、カリフォルニア州のブロックグループ（米国国勢調査局がサンプルデータを公開する最小の地理的単位）の20,640サンプルを含み、それぞれに8つの特徴量（中央値収入、中央値家屋年齢、平均部屋数など）があります。Kerasからload_dataを使用してデータセットをダウンロードしましょう。

import keras
 
(X_train, y_train), (X_test, y_test) = keras.datasets.california_housing.load_data(
    version="large", path="california_housing.npz", test_split=0.2, seed=113
)

データセットはすでに訓練データとテストデータに分割されており、numpy配列として提供されていますので、さらに前処理を行う必要はありません。次に、ペアプロットを使用して特徴量間の関係を視覚化してみましょう。

columns = [
    "Longitude", "Latitude", "AveOccup", "Population", 
    "AveBedrms", "AveRooms", "HouseAge", "MedInc"
]
 
train_df = pd.DataFrame(X_train, columns=columns)
train_df['MedHouseVal'] = y_train
 
sns.pairplot(train_df)
plt.show()

このプロットは解釈が難しいですが、いくつかの洞察を得ようと試みます。平均部屋数、人口、および平均世帯人数に関して、いくつかの外れ値があるように見え、それがいくつかのプロットの解釈を難しくしているようです。これらの外れ値は、 Kerasのドキュメントに部分的に説明されており、バケーションリゾートのような空き家を含むいくつかの家庭が特定の列で驚くほど高い数値を持つ傾向があるとされています。また、これらの変数の間には強い正の相関関係があり、すべての特徴量を予測に含める必要がない(予測に重要でない)ことが示唆されています。

私たちが予測しようとしている住宅の中央値価値に関しては、ほとんどの特徴量とは明確な関係が見られないようです（実世界で予想される通り）。ただし、中央値収入にはわずかに正の相関関係が見られます。収入が高い家庭が価値の高い住宅に住む傾向があるのは当然でしょう。緯度と経度のパターンをより詳しくみるため、プロットしてみましょう。

plt.figure(figsize=(10, 6))
scatter = plt.scatter(train_df['Longitude'], train_df['Latitude'], 
                      c=train_df['MedHouseVal'], cmap='viridis', s=10)
plt.colorbar(scatter, label='Median House Value')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.title('Latitude and Longitude with Median House Value')
plt.show()

ご覧の通り、海岸線沿いの住宅はより高い価値を持つ傾向があるようです。（この仮説を厳密にテストするためには統計分析が必要ですが、ここでは最適なXGBoost予測モデルを構築することが目的なので、その分析は行いません。）したがって、XGBoostがこのパターンを捉えることが期待できます。

ステップ3. モデル

交差検証を行ってハイパーパラメータのチューニングをする前に、基本的なXGBoostモデルにおけるハイパーパラメータを理解する必要があります。以下は、XGBoostパッケージで提供されているハイパーパラメータの一部を示しています。

XGBoostモデルにおけるハイパーパラメータ
ハイパーパラメータ	引数名	説明
学習率	`eta`	各木の出力値に掛け算される、徐々に学習するためのパラメータです ( $\nu$ )。デフォルトは0.3です。
ガンマ ( $\gamma$ )	`gamma`	分割による最小のゲインを設定するもので、剪定に使用されます。デフォルトは0です。
ラムダ ( $\lambda$ )	`lambda`	出力値と類似度スコアの分母に加えられるL2正則化率です。デフォルトは1です。
アルファ ( $\alpha$ )	`alpha`	L1正則化率です。L1正則化項を追加することもできますが、過学習が深刻でない限り、通常L2正則化を使用します。デフォルトは0です。
最大深度	`max_depth`	許可される木の最大深度です。デフォルトは6です。
最小の重み	`min_child_weight`	リーフノードに対して許可される最小の重みです。回帰の場合、重みは残差の数であり、分類の場合は $\sum p(1-p)$ です。デフォルトは1です。
サブサンプル	`subsample`	木の分岐を見つけるために使用されるサンプル（行）の比率です。デフォルトは1で、すべてのサンプルが使用されるように設定されています。
列サブサンプリング	`colsample_by*`	木の分岐を見つけるために使用される特徴量（列）の比率です。ツリー、レベル、ノードのいずれか毎にサンプリングできます。デフォルトは、`bytree`、`bylevel`、`bynode` は全て1で、すべての分割で全ての特徴量が利用されるように設定されています。
木の構築法	`tree_method`	すべての可能な分岐をテストする（`exact`）か、重み付き分位数スケッチ（`approx` または `hist`）を使用するかを指定します。デフォルトは `auto` に設定されています。

他にも max_leaves や min_bin などの、各木で許可される最大の葉の数や、重み付き分位数スケッチで使用されるビンの数を指定するハイパーパラメータもありますが、上記が主なものです。データセットが大きくないため、主に eta、gamma、および lambda をチューニングします。（max_depth と min_child_weight は剪定のためのもので、gamma と lambda で管理できます。subsample と colsample_by* は大規模なデータセットを扱うためのものです。）

次に、交差検証のための関数を定義します。今回は、sklearn の RandomSearchCV を使用します。これはハイパーパラメータの確率分布を取り、ランダムにハイパーパラメータの組み合わせを選んでスコアを比較するものです。GridSearchCV を使用することもできますが、時間がかかりすぎると判断しました。

from sklearn.model_selection import StratifiedKFold, RandomizedSearchCV
from xgboost import XGBRegressor, plot_importance
from sklearn.metrics import make_scorer, mean_squared_error
 
def randomsearchCV_XGBR(X_train, y_train, parameters, n_iter=20, n_splits=10, shuffle=True):
  estimator = XGBRegressor(objective='reg:squarederror',
                           n_estimators=100,
                           tree_method='hist',
                           seed=42)
 
  cv = StratifiedKFold(n_splits=n_splits, shuffle=shuffle)
 
  rand_search = RandomizedSearchCV(
      estimator=estimator,
      param_distributions=parameters,
      scoring = make_scorer(mean_squared_error),
      n_iter=n_iter, ## number of samples to test
      ## n_jobs=10, ## if multithreads are available
      cv = cv,
      verbose=0
  )
 
  rand_search.fit(X_train, y_train)
 
  return rand_search.best_params_, rand_search.best_score_

ハイパーパラメータの確率分布を渡すために、scipy.stats の uniform 関数(一様分布)を使用します。一様分布を使用するのはどのハイパーパラメータの値が最適な結果をもたらすのかについて何もわかっていないためです。

from scipy.stats import uniform
 
parameters = {
      'eta': uniform(loc=0.001, scale=1),
      'gamma': uniform(loc=0, scale=1),
      'lambda': uniform(loc=0, scale=1)
}
 
params, score = randomsearchCV_XGBR(X_train, y_train, parameters)
 
print("best parameters: ", params)
print("best score: ", score)
 
# =>
# best parameters:  {'eta': 0.013192469939280582, 'gamma': 0.4946327995662789, 'lambda': 0.1587494731756135}
# best score:  4679027302.4

上記のハイパーパラメータから、XGBoost Regressor を構築し、訓練データをモデルに学習させます。

xgb = XGBRegressor(objective='reg:squarederror',
                  eta=params['eta'],
                  gamma=params['gamma'],
                  reg_lambda=params['lambda'],
                  n_estimators=100,
                  tree_method='hist',
                  seed=42)
 
r = xgb.fit(X_train, y_train)

ステップ4. モデルの評価

モデルを評価するために、先ほど学習させたモデルを使ってテストデータに対して予測を行います。

pred = xgb.predict(X_test)
 
mean_squared_error(y_test, pred)
# => 4806821000.0

この場合、MSE（平均二乗誤差）は解釈が難しいため、r2_score または決定係数を使用しましょう。決定係数は、予測変数の変動のどれだけが説明変数によって説明されるかを表します。（詳しく知りたい場合は、StatQuestの R-squared, Clearly Explained!!! をチェックすることをお勧めします。）

from sklearn.metrics import r2_score
 
r2_score(y_test, pred)
# => 0.6384312017848871

同じデータに対して多くの前処理を行い、線形回帰モデルをフィットさせたところ、決定係数は約0.6でした。したがって、XGBoost Regressorはそれほど手間をかけずに、より良い予測を達成したと言えます！（より厳密に確認するためには、パフォーマンスを複数回測定する必要があります。）

XGBoostのもう一つの利点は、木や特徴の重要性を簡単に可視化できることです。ここでは、木ではなく特徴の重要性を可視化します（木は大きすぎてここでは可視化しないことにしました）。

F_score は、その特徴を使ってどれだけ多くの分割が行われたかを示す指標です。上記の図は、緯度、経度、中央値収入を使った分割が頻繁に行われたことを示しており、これはステップ1で観察した分離性と関係性を考えると理解できます。他の特徴の重要性は比較的低く、これもステップ1での予想と一致しています。

結論

XGBoostを使った機械学習パイプラインがどのようなものか、そしてそれをPythonでどのように実装するかについて理解していただけたでしょうか。理解できたという方は、知識とスキルが向上のためXGBoostを使って分類問題に取り組んでみることをお勧めします。

追記：Cプログラマへの道シリーズをフォローしている方は、公式のXGBoostのドキュメンテーションを見ながら、C言語でXGBoostを実装することに挑戦してみてください。頑張ってください！

リソース

Keras. n.d. California Housing price regression dataset. Keras.
n.d. 2022. XGBoost Documentation. XGBoost.