Pythonでデータ分析：XGboost - データサイエンティスト(仮)

導入

前回、アンサンブル学習の方法の一つであるランダムフォレストについて紹介しました。
tekenuko.hatenablog.com
今回は、XGboostと呼ばれる、別の方法がベースになっているモデルを紹介します。

XGboostとは

XGboostは、アンサンブル学習がベースになっている手法です。
アンサンブル学習は、大きく2通りの方法があります。一つはバギングと呼ばれ、復元抽出によってたくさんモデルを作成し、それらを平均化する手法です。各モデルが無相関であれば、平均化によりモデルの平均誤差が1 / （モデルの数）になります。ランダムフォレストは、このバギングがベースになっています。ただし、各決定木ができるだけ独立になるよう、使う変数群もランダム抽出するなどの工夫がなられています。前回は、過学習気味である結果はあまり解消できませんでしたが、高い精度を持つモデルが構築できたのでした。

もう一つは、ブースティングと呼ばれ、モデルを逐次更新していく手法です。更新は、例えばAdaBoostと呼ばれる方法では、データにかける「重み」を前回作成したモデルをもとに変えることで行います。その際、作成したモデルの「正確さ」に対応するような重み係数も同時に算出し、最後に各モデルの予測値と重み係数をかけ合わせたものの和を、全体の予測値とします。

このブースティングという枠組みは、損失関数の最小化問題でも議論できることが（例えば）以下の論文で知られています*1。
projecteuclid.org
損失関数の最小化問題は、数理最適化問題であり、目的に応じて最小値を求めるための様々な手法があります。よく知られた（使われる）方法の一つに、損失関数の勾配情報（微分）を使ったアルゴリズムがあります*2。ブースティングで勾配情報を用いたアルゴリズムを用いているものを勾配ブースティングと呼び、それをC++で高速実装したものがXGboostです。

このXGboostですが、KaggleのHiggs Boson Machine Learning Challengeの優勝チームが使用したことで、有名になったようです。ちなみにちょっとだけ補足しておくと、Higgs Bosonというのは素粒子の標準模型に含まれる素粒子の一つ（だと現在は理解されています）です。標準模型とは、現在確立している、もっとも基本的な素粒子模型です。ほとんどの素粒子に関する実験結果を説明できる優れた理論ですが、理論的・実験的に問題があることも知られています。そのため、標準模型を拡張した模型の構築、およびそのような理論の実験での探索が素粒子物理学の研究内容の一つです*3。この標準模型の持つ性質の一つに、電弱対称性の自発的破れというものがあり、その際の名残りとして発生する粒子がHiggs Bosonです*4。現在の我々の世界でHiggs Bosonを発生させる方法の一つに、高いエネルギーを持った素粒子を衝突させる、という方法があり、発生したHiggs Bosonは様々な素粒子へと崩壊します（崩壊の仕方は、考える素粒子模型によって異なります。）。その中でも、2つのタウ粒子という、電子の親戚みたいな素粒子に崩壊する現象（ $h \rightarrow \tau \tau$ ）に関するデータがKaggleに提供されたのでした。
www.kaggle.com

このとき提供されたデータは、伝統的な素粒子の解析だとHiggs Bosonが存在するとは言い切れない程度のデータ量でした。Higgs Bosonの存在が初めて検証されたのは、別の崩壊現象です。タウ粒子への崩壊現象でHiggs Bosonの発見がされなかったのは、タウ粒子では高感度の検出が難しいためです。タウ粒子は質量が大きく、ハドロンと呼ばれる強い相互作用をする粒子への崩壊を許してしまいます。一般にハドロンは複雑な崩壊をするため、タウ粒子も複雑な崩壊パターンになります。しかし、それはこの崩壊の解析が無駄であることを意味しているわけではありません。標準模型の無矛盾性、もしくは新たな理論の可能性の検証のために、この崩壊の解析を進めることは一定の重要性を持っています。

脱線しすぎたので話を戻します（笑）。
このXGboost（というかブースティング全般に関して）ですが、なぜか汎化性能がよいそうです。もう少し詳しくいうと、ブースティングでは統計的学習理論で定められた（緩い）誤差の上限が求まるのですが、実際にはそれよりも性能が良くなることが多い、とのことです。なぜこのようなことが起きるかは、自分はあまり理解できていません。ひとまず、今回はXGboostという手法が非常に良い精度を出すらしいので、ちょっとモデル作ってみよう、という立場で話を進めます。

インストール方法（ざっくり）

PythonでXGboostと使うためには、以下のサイトを参考にインストールします。
xgboost/python-package at master · dmlc/xgboost · GitHub

私はMacユーザなので、そこまで問題はありませんでしたが、Window（特に32bit）に入れようとすると闇が深そうです。インストール方法に関しては、色々なサイトでご紹介されているので、ここではあまり深く触れないことにします。

データセット

毎度おなじみ、ボストン近郊の住宅情報のデータを使用します。

# パッケージのインポート
import numpy as np
import pandas as pd
from pandas import DataFrame
import xgboost as xgb
from sklearn.datasets import load_boston
# データマート作成
boston = load_boston()
df = DataFrame(boston.data, columns = boston.feature_names)
df['MEDV'] = np.array(boston.target)

モデル構築

デフォルトのパラメータを用いて構築

XGboostは非常に多くのパラメータを持っています。まずはデフォルトの状態でどこまで性能のよいモデルができるかを見てみます。

# 説明変数、目的変数
X = df.iloc[:, :-1].values
y = df.loc[:, 'MEDV'].values
# 学習用、検証用データ作成
from sklearn.cross_validation import train_test_split
(X_train, X_test, y_train, y_test) = train_test_split(X, y, test_size = 0.3, random_state = 666)

XGboostによる回帰は、xgboostのXGBRegressorで行います。

# XGboostのライブラリをインポート
import xgboost as xgb
# モデルのインスタンス作成
mod = xgb.XGBRegressor()
mod.fit(X_train, y_train)

平均二乗誤差（MSE）と $R^2$ がどうなるかを見てみます。

y_train_pred = mod.predict(X_train)
y_test_pred = mod.predict(X_test)
# MSE
from sklearn.metrics import mean_squared_error
print('MSE train : %.3f, test : %.3f' % (mean_squared_error(y_train, y_train_pred), mean_squared_error(y_test, y_test_pred)) )
# R^2
from sklearn.metrics import r2_score
print('R^2 train : %.3f, test : %.3f' % (r2_score(y_train, y_train_pred), r2_score(y_test, y_test_pred)) )

# 出力
>>>print('MSE train : %.3f, test : %.3f' % (mean_squared_error(y_train, y_train_pred), mean_squared_error(y_test, y_test_pred)) )
MSE train : 1.687, test : 10.921
>>>print('R^2 train : %.3f, test : %.3f' % (r2_score(y_train, y_train_pred), r2_score(y_test, y_test_pred)) )
R^2 train : 0.981, test : 0.847

今までで最も良いモデルになりました。学習データの方に適合しがちな性質は残っていますが、ランダムフォレストよりもMSEが小さくて $R^2$ が大きいです。究極的にはなぜこうなるかわかりませんが（汗）、すごいです。

残差プロット

モデルの性能を大雑把に見るのに、残差を可視化してみます。

# matplotlibのインポートとおまじない
import matplotlib.pyplot as plt
%matplotlib inline

# プロット
plt.figure(figsize = (10, 7))
plt.scatter(y_train_pred, y_train_pred - y_train, c = 'black', marker = 'o', s = 35, alpha = 0.5, label = 'Training data')
plt.scatter(y_test_pred, y_test_pred - y_test, c = 'lightgreen', marker = 's', s = 35, alpha = 0.7, label = 'Test data')
plt.xlabel('Predicted values')
plt.ylabel('Residuals')
plt.legend(loc = 'upper left')
plt.hlines(y = 0, xmin = -10, xmax = 50, lw = 2, color = 'red')
plt.xlim([-10, 50])
plt.show()

f:id:tekenuko:20160922115049p:plain

ランダムフォレストの場合は、学習データに関して右下の方にランダムっぽくないものがあったのですが、XGboostではそれが原点に近づいています。そのため、精度は上がったと思われます。しかしながら、ランダムっぽくはない（原点に寄っただけ）のでXGboostで以前拾えていなかった情報を拾えるようにはなっていないのではないかというのが個人的な意見です。

グリッドサーチでより性能のよいモデルを探す

XGboostのパラメータに関しては、以下に載っています。
xgboost/parameter.md at master · dmlc/xgboost · GitHub

今回は、幾つかのパラメータに関して、グリッドサーチを試してみます。前回同様、scikit-learnのGridSearchCVで行います。

# グリッドサーチに必要なクラスのインポート
from sklearn.grid_search import GridSearchCV
# サーチするパラメータは範囲を指定
params = {'max_depth': [3, 5, 10], 'learning_rate': [0.05, 0.1], 'max_depth': [3, 5, 10, 100], 'subsample': [0.8, 0.85, 0.9, 0.95], 'colsample_bytree': [0.5, 1.0]}
# モデルのインスタンス作成
mod = xgb.XGBRegressor()
# 10-fold Cross Validationでパラメータ選定
cv = GridSearchCV(mod, params, cv = 10, scoring= 'mean_squared_error', n_jobs =1)
cv.fit(X_train, y_train)

MSEと $R^2$ は以下のようになります。

y_train_pred = cv.predict(X_train)
y_test_pred = cv.predict(X_test)
# MSE
from sklearn.metrics import mean_squared_error
print('MSE train : %.3f, test : %.3f' % (mean_squared_error(y_train, y_train_pred), mean_squared_error(y_test, y_test_pred)) )
# R^2
from sklearn.metrics import r2_score
print('R^2 train : %.3f, test : %.3f' % (r2_score(y_train, y_train_pred), r2_score(y_test, y_test_pred)) )

# 出力
>>>print('MSE train : %.3f, test : %.3f' % (mean_squared_error(y_train, y_train_pred), mean_squared_error(y_test, y_test_pred)) )
MSE train : 1.742, test : 10.375
>>>print('R^2 train : %.3f, test : %.3f' % (r2_score(y_train, y_train_pred), r2_score(y_test, y_test_pred)) )
R^2 train : 0.981, test : 0.855

学習用データの場合は、デフォルトの場合と比べて少しMSEが悪くなり、検証用データの方は少し良くなりました。 $R^2$ の方を見ても、検証用データに関して若干上昇していることから、パラメータ調整によって少し汎化性能が良くなったようです。

残差もプロットしておきましょう。

plt.figure(figsize = (10, 7))
plt.scatter(y_train_pred, y_train_pred - y_train, c = 'black', marker = 'o', s = 35, alpha = 0.5, label = 'Training data')
plt.scatter(y_test_pred, y_test_pred - y_test, c = 'lightgreen', marker = 's', s = 35, alpha = 0.7, label = 'Test data')
plt.xlabel('Predicted values')
plt.ylabel('Residuals')
plt.legend(loc = 'upper left')
plt.hlines(y = 0, xmin = -10, xmax = 50, lw = 2, color = 'red')
plt.xlim([-10, 50])
plt.show()

f:id:tekenuko:20160922120440p:plain

見た目だと改善があまりわかりません（笑）。

おわりに

今回、XGboostというモデルを用い、実際に精度の良いモデルを構築することができました。しかしながら、

データの構造を捉えきれていない（気がする）
学習データに適合しがち

という傾向はいぜんとして残っています。こちらは、原因はそれぞれ独立なのか、それとも共通の原因があるのか、などまだまだ問題点としても検討しきれていません。これらの問題を解決できるかを考えてみて、できる/緩和できる方法があったら今後紹介してみたいと思います。

*1:PRMLの14.3.1でも少し紹介されています。

*2:最急降下法やニュートン法などです。これらは、目的関数が微分可能かつ凸であれば高速に最適値を求めることができる手法です。

*3:私は標準模型を「超対称化」した模型の現象論の研究をしていました。ようするに理論的立場から実験への橋渡しをするような研究です。

*4:電弱対称性の破れを引き起こすのはHiggs場というもので、Higgs Bosonとは意味が異なります。この辺は、説明するととても長くなるのは、このへんでお茶をにごしておきます（汗）。