FastBDTの計算時間が速いかを確認してみる
はじめに
最近、ブースティング系のアルゴリズムでXGboostより速いものが実装されているようです。
github.com
論文を見ると、Stochastic Gradient Boosted Decision Tree (SGBDT)という手法が用いられているようです*1。用途としては、信号(例えば、ラベル1)と背景事象(1以外)を分類する場合が想定されています。このアルゴリズムは、素粒子実験の一つである、Belle II実験の解析で使用されているとのことです。ばくっと探していたらカールスルーエの学生の修論があったので、載せておきます。
https://ekp-invenio.physik.uni-karlsruhe.de/record/48602/files/EKP-2015-00001.pdf
Belle II実験
Belle II実験に関しては、修士課程のときの自分の研究内容と深く関わっているので少し補足を。
Belle II実験とは、日本のKEKにある電子-陽電子加速器を用いた実験です。この実験の目的の一つ*2は、B中間子とよばれる、ボトムクォークを含む複合粒子をたくさん生成し、その「稀な」崩壊現象の精密測定から新しい理論の探索を行うことです*3。標準模型を超える物理があった場合、直接高エネルギー加速器で直接生成されるだけでなく、間接的に様々な物理量に寄与する可能性があります*4。その寄与は、標準模型の予測値からも「ずれ」という形で表れます。その「ずれ」を実験で精密に検証することで、新しい物理の姿を間接的に見ることができます。
間接的に新しい物理の効果を見るにあたって、「稀な」崩壊現象を見ることが非常に大事です。それは、比較的起こりやすい崩壊現象よりも、新しい物理の効果を捉えやすいためです。標準模型では、例えばGlashow‐Iliopoulos‐Maiani機構(GIM機構)といった、ある種類の崩壊を抑制する機構が備わっているため、素朴な次元解析での評価よりも崩壊率が小さいものがあります。これらは標準模型では「稀な」崩壊と予測されます。しかし、新しい物理は必ずしもGIM機構のようなものが備わっているわけではないため*5、大きな補正が「稀な」崩壊現象に加わる可能性があります。そのため、比較的起こりやすい崩壊現象とくらべて変化を捉えやすくなります。
B中間子に着目する理由は、その「稀さ」が絶妙だからです。B中間子の物理現象は、抑制されているカテゴリーの中では比較的大きな寄与を持っています*6。そのため、稀ではあるが、標準模型が正しいと仮定した場合の予測値が実験で検知できる程度には寄与がある、といった現象になり、「標準模型からのずれ」の検証が他の中間子よりも容易になっています。
インストールから使用可能になるまで
つい素粒子の話になると長くなってしまうので、話を戻します。Mac OSの場合、サイトから直接ダウンロードし、makeすれば使えるようになります。cmakeは予めインストールしておきます。
$ git clone https://github.com/thomaskeck/FastBDT $ cd FastBDT/ $ cmake . $ make $ make install
上記の操作をすれば、実行ファイルが生成されるので、FastBDTのディレクトリで
$ ./FastBDF [あといろいろ必要なオプション指定]
と唱えると実行できます。
Pythonから使用したい場合
Pythonで動かす方法もあります。ただし、まだ完全に設定などを理解していないので、現状を備忘録的にまとめておきます。
修正点
まず、インポートする元は ~/FastBDT/python/FastBDT.pyにあります。ただし、このスクリプトで指定しているライブラリの名前が自分の環境では異なっていたので、少し修正が必要です。修正方法は2通りあります。
- [19行目]で指定しているライブラリをlibFastBDT_CInterface.soからlibFastBDT_CInterface.dylibに書き換える
- [14行目-16行目]のコメントアウトを外し、[19行目]をコメントアウト
こうすれば、インポートしたときに正しく読み込みができるようになりました。
インポート
以下のようにパスを上手く指定してインポートします。
import sys sys.path.append('/Users/[ユーザ名]/FastBDT/python') import FastBDT
インポートまでの一連の流れですが、不慣れなこともあり、かなり洗練されていない感じになっています(というよりセオリーを知らない。)が、ご了承ください(汗)。
モデルはClassifierで作成することができます。モデルインスタンス作成後は、scikit-learnのようにfitやpredictといったメソッドが使えます。
計算時間の比較
学習にかかった時間を計測します。比較対象として、パーセプトロン、ランダムフォレスト、XGboostを選択します。
サンプルデータ
FastBDTは'1'とそれ以外を識別する問題にのみ対応しているようです*7。そのため、簡単のため、人工的に2値データを生成して、それで比較することにしました。クラス分類のテスト用の人工データは、scikit-learn.datasetのmake_classificationで作成することが可能です。参考にしたサイトは以下です。
overlap.hatenablog.jp
説明変数が5つ(n_features = 5)、サンプル数100万(n_samples = 1000000)の2値データ(n_classes = 2)を生成します。他のオプションは適当です。
from sklearn.datasets import make_classification # データ作成 dat = make_classification( n_samples=1000000, n_features=5, n_informative=2, n_redundant=2, n_repeated=0, n_classes=2, n_clusters_per_class=2, weights=None, flip_y=0.01, class_sep=1.0, hypercube=True, shift=0.0, scale=1.0, shuffle=True, random_state=None)
データマートを作成し、さらに学習用、検証用データを作成します。
# ライブラリ import pandas as pd from pandas import DataFrame # DataFrameに格納 df = DataFrame(dat[0], columns = ['Var1', 'Var2', 'Var3', 'Var4', 'Var5']) df['Class'] = dat[1] # 説明変数、目的変数 X = df.iloc[:, :-1].values y = df.loc[:, 'Class'].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)
計算時間、誤分類数、Accuracyを算出
ばくっと計算する関数を作っておきます。
def CalcTimes(mod, X_train, X_test, y_train, y_test): import numpy as np from sklearn.metrics import accuracy_score start = time.time() # モデルインスタンス作成、学習 mod = mod mod.fit(X_train, y_train) elapsed_time = time.time() - start print('elapsed_time: %.2f' %elapsed_time + ' [sec]') # 誤分類数、Accuracy y_train_pred = np.round(mod.predict(X_train)) y_test_pred = np.round(mod.predict(X_test)) print('Misclassified samples: %d' % (y_test != np.round(y_test_pred)).sum()) print('Accuracy: %.2f' % accuracy_score(y_test, y_test_pred))
FastBDT
from FastBDT import Classifier mod = Classifier() CalcTimes(mod, X_train, X_test, y_train, y_test) # 出力 elapsed_time: 6.50 [sec] Misclassified samples: 34948 Accuracy: 0.88
パーセプトロン
from sklearn.linear_model import Perceptron mod = Perceptron() CalcTimes(mod, X_train, X_test, y_train, y_test) # 出力 elapsed_time: 0.56 [sec] Misclassified samples: 66686 Accuracy: 0.78
ランダムフォレスト
from sklearn.ensemble import RandomForestClassifier mod = RandomForestClassifier() CalcTimes(mod, X_train, X_test, y_train, y_test) # 出力 elapsed_time: 44.74 [sec] Misclassified samples: 39241 Accuracy: 0.87
XGboost
from xgboost import XGBClassifier mod = XGBClassifier() CalcTimes(mod, X_train, X_test, y_train, y_test) # 出力 elapsed_time: 45.03 [sec] Misclassified samples: 34131 Accuracy: 0.89
比較結果
今回のデータセットに関しては、Macbook Air使用で、FastBDTはランダムフォレストやXGboostと比較して6倍以上速い結果になりました。しかも、Accuracyが大きく損なわれずにです。条件を変えると倍率は変わると思われますが、概ね似たカテゴリーのアルゴリズムよりは高速なんだと思います。もっと開発が進んで使いやすくなったら、非常に強力なアルゴリズムなのではないかと期待されます。
*1:細かい部分で工夫があるようですが、まだキャッチアップできていません。
*2:他にもK中間子の稀崩壊探索なども目的です。
*3:前身として、Belle実験というものが同じKEKで行われていました。Belle実験では、B中間子の崩壊現象を用いた標準模型の検証が主目的でした。
*4:量子効果(摂動論の高次効果)で寄与します。
*5:4次元以外の「余剰な」次元を持っている理論がその例です。超対称性を持った理論でSuper GIM機構というものがある場合は、余剰次元をもつ理論よりは寄与が抑制されます。
*6:色々省略しますが、これはトップクォークという素粒子がとてつもなく重いことが起源です。
*7:使ってみた結果、こう認識しています。間違いがあればご指摘ください。
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つのタウ粒子という、電子の親戚みたいな素粒子に崩壊する現象()に関するデータが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)とがどうなるかを見てみます。
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が小さくてが大きいです。究極的にはなぜこうなるかわかりませんが(汗)、すごいです。
残差プロット
モデルの性能を大雑把に見るのに、残差を可視化してみます。
# 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()
ランダムフォレストの場合は、学習データに関して右下の方にランダムっぽくないものがあったのですが、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とは以下のようになります。
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が悪くなり、検証用データの方は少し良くなりました。の方を見ても、検証用データに関して若干上昇していることから、パラメータ調整によって少し汎化性能が良くなったようです。
残差もプロットしておきましょう。
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()
見た目だと改善があまりわかりません(笑)。
おわりに
今回、XGboostというモデルを用い、実際に精度の良いモデルを構築することができました。しかしながら、
- データの構造を捉えきれていない(気がする)
- 学習データに適合しがち
という傾向はいぜんとして残っています。こちらは、原因はそれぞれ独立なのか、それとも共通の原因があるのか、などまだまだ問題点としても検討しきれていません。これらの問題を解決できるかを考えてみて、できる/緩和できる方法があったら今後紹介してみたいと思います。
Pythonでデータ分析:ランダムフォレスト
導入
前回、非線形的な効果を表現することの一例として、決定木回帰を紹介しました。
tekenuko.hatenablog.com
決定木は、ざっくりとしたデータの特徴を捉えるのに優れています*1。しかしながら、条件がデータに依存しがちなため、過学習しやすいという欠点もあったのでした。
この欠点を緩和するための方法として、アンサンブル学習という方法があります。これは、データをサンプリングしてモデルを構築、それらを組み合わせて新たなモデルを構築する方法です。ランダムフォレストは、ざっくりいうと多数の決定木を作成し、それらを平均化する手法です。個々のモデルではデータの変化の影響が大きくても、まるっと平均化したものは影響が少なくなるため、一つの決定木でモデルを作るのに比べて過学習が緩和されやすくなります*2。
ランダムフォレストをより深く理解するためには、ある程度しっかりした機械学習の本を読んだりする必要があります。ここでは、Pythonでどのようにモデル構築ができるかについてフォーカスし、理論的は背景には深く立ち入らないことにします*3。
参考
毎度ですが、以下です。今更ですが、とても勉強になります。
book.impress.co.jp
データセット
今回も、前回までと同じく、ボストン近郊の住宅情報のデータを使用します。
# 必要なライブラリのインポート import numpy as np import pandas as pd from pandas import DataFrame from sklearn.datasets import load_boston # データのロード、マージ boston = load_boston() df = DataFrame(boston.data, columns = boston.feature_names) df['MEDV'] = np.array(boston.target)
モデル構築
デフォルトのパラメータでモデルを構築してみる
ランダムフォレストは、多くの調整できるパラメータがあります。まずは、デフォルトでどの程度の性能になるかを見てみます。
まず、説明変数と目的変数を定義し、学習用と検証用データに分けます。モデル構築には、学習用データを使用します。
# 説明変数、目的変数 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)
ランダムフォレストによる回帰は、scikit-learnのRandomForestRegressorで行います。
# 必要なライブラリのインポート from sklearn.ensemble import RandomForestRegressor # モデル構築、パラメータはデフォルト forest = RandomForestRegressor() forest.fit(X_train, y_train)
平均二乗誤差(MSE)とがどうなるかを見てみます。
# 予測値を計算 y_train_pred = forest.predict(X_train) y_test_pred = forest.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('MSE 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 : 2.546, test : 14.968 >>>print('MSE train : %.3f, test : %.3f' % (r2_score(y_train, y_train_pred), r2_score(y_test, y_test_pred)) ) MSE train : 0.972, test : 0.790
今まで考えてきた手法と比較すると、かなり性能が上がっています。残念ながら、過学習する傾向は残っています。それでも検証用データに対するは約0.8と大きい値であり、悲観するほどではないかなと思います。
残差プロット
モデルの性能を大雑把に見るのに、残差を可視化してみます。
# 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()
図を見ると、検証データよりも学習データに関するほうが散らばりは少なく、実際に学習データにより適合したモデルであることがわかります。さらに、学習データに関して、右下にランダムに見えない振る舞いが見られ、ランダムフォレストでも完全に捉えきれていない情報がある可能性があります。
グリッドサーチでより性能の良いモデルを探す
さらにモデルの性能を上げる方法として、ハイパーパラメータの調整があります。ハイパーパラメータとは、モデルに事前に手で投入するパラメータのことで、これを変化させると精度が上がる可能性があります。ハイパーパラメータの調整方法は、大雑把に以下の方法があります。
- グリッドサーチ:決められた範囲を総当り探索
- ランダムサーチ:投入するパラメータをランダムに変えて試す
今回は、グリッドサーチを試してみます。グリッドサーチは、scikit-learnのGridSearchCVで簡単かつ汎用的に行うことができます。
# 必要なライブラリのインポート from sklearn.grid_search import GridSearchCV # 動かすパラメータを明示的に表示、今回は決定木の数を変えてみる params = {'n_estimators' : [3, 10, 100, 1000, 10000], 'n_jobs': [-1]}
実際にGridSearchCVで探索してみます。検証する手法は交差検証、評価する指標としてMSEを用います。
# モデルにインスタンス生成 mod = RandomForestRegressor() # ハイパーパラメータ探索 cv = GridSearchCV(mod, params, cv = 10, scoring= 'mean_squared_error', n_jobs =1) cv.fit(X_train, y_train)
MSEとは以下のようになります。
# 予測値を計算 y_train_pred = forest.predict(X_train) y_test_pred = forest.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('MSE 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.712, test : 13.053 >>>print('MSE train : %.3f, test : %.3f' % (r2_score(y_train, y_train_pred), r2_score(y_test, y_test_pred)) ) MSE train : 0.981, test : 0.817
デフォルトの場合と比べて、若干改善しました。さらに変化させるパラメータや範囲を増やすとより性能が良いモデルができる可能性はありますが、計算時間も爆発的に増加します。総当り的に探索するときは、最初はある程度ざっくり探索して精度向上が見込めそうな領域を見つけ*4、そこをさらに細かく探す、といった方法を取るのがよいでしょう。
Pythonでデータ分析:決定木
導入
前回、線形回帰からの拡張の一つとして、非線形項をモデルに加えることを紹介しました。
tekenuko.hatenablog.com
非線形性を表現する方法は他にも幾つかあり、その一つに、決定木という手法があります。今回は、回帰に決定木を用いた方法を紹介します。
決定木とは
まず、分類問題の文脈での決定木から定義します。決定木とは、変数に対して、分割が「お得になる」ような多数決を終了条件を満たすまで繰り返す手法です。「お得かどうか」の判断は、例えば情報の乱雑さに対応するエントロピーという量が使用されます。よりエントロピーが低い(= 分類がきれい)な条件を優先的に終了条件を満たすまで選択していきます。このように、「Xという変数がY以上」という条件でデータを分ける、といった操作が複数起こるため、その条件分岐を可視化すると、木の枝葉のような構造になります。
回帰における決定木は、上記のエントロピーの代わりに、平均二乗誤差(MSE)を指標として用います。サンプルをある条件でサブセットに分割し、それぞれのサブセットの平均値を予測値とし、MSEが減少が大きくなる分割を優先的に採用していく、といったアルゴリズムになります。そのため、予測値を可視化すると、非常に「カクカクした」(微分不可能な)振る舞いを示します。
今回は、後者の決定木回帰を議論していきます。
参考
いつものことながら、以下です。
book.impress.co.jp
データセット
今回も、前回までと同じく、ボストン近郊の住宅情報のデータを使用します。
# 必要なライブラリのインポート import numpy as np import pandas as pd from pandas import DataFrame from sklearn.datasets import load_boston # データのロード、マージ boston = load_boston() df = DataFrame(boston.data, columns = boston.feature_names) df['MEDV'] = np.array(boston.target)
モデル構築
木の分岐の深さを3に固定した場合
決定木回帰は、scikit-learnで簡単に行うことが可能です。今回は、決定木の特徴を見ることにフォーカスするために、説明変数はLSTATの一つだけに絞ります。
モデル構築は、scikit-learnのDecisionTreeRegressorクラスを用いて行うことが可能です。モデリング例は以下になります。
# 必要なライブラリのインポート from sklearn.tree import DecisionTreeRegressor # 説明変数と目的変数を選択 X = df.loc[:, ['LSTAT']].values y = df.loc[:, 'MEDV'].values # モデル構築、木の深さは3に固定 tree = DecisionTreeRegressor(max_depth = 3) tree.fit(X, y)
散布図と回帰直線をプロットすると、以下になります。
# 回帰直線を図示するのに変数を並び替え sort_idx = X.flatten().argsort() # matplotlibのインポートとおまじない import matplotlib.pyplot as plt %matplotlib inline # プロット plt.figure(figsize = (10, 7)) plt.scatter(X[sort_idx], y[sort_idx], c = 'blue', label = 'Training points') plt.plot(X[sort_idx], tree.predict(X[sort_idx]), color = 'red', label = 'depth = 3') plt.xlabel('LSTAT') plt.ylabel('MEDV') plt.legend(loc = 'upper right') plt.show()
線形回帰や非線形回帰と比較して、カクカクした直線が引かれます。その「カクカクする」性質のため、ある程度データの傾向を補足できるのが決定木のメリットです。一方で、予測値が微分不可能になってしまう制約があるため、解析的には扱いづらい手法になっています。
木の深さと予測値の関係
先程は、木の深さを3に固定しました。木の深さを大きくすると、より複雑な条件分岐を許すことになるため、モデルの当てはまりが良くなります。しかし、それは過学習を引き起こすことにもつながるため、適切な値を選択する必要があります。
以下で、木の深さを幾つか変えた場合の予測値をプロットしています。
# あらかじめ色のリストを作成 color = ['red', 'green', 'yellow', 'magenta', 'cyan'] # for文で順繰り木の深さを変えたモデル結果をプロット plt.figure(figsize = (10, 7)) plt.scatter(X[sort_idx], y[sort_idx], c = 'lightgray', label = 'Training points') for t in (np.arange(5)): tree = DecisionTreeRegressor(max_depth = t + 1) tree.fit(X, y) sort_idx = X.flatten().argsort() plt.plot(X[sort_idx], tree.predict(X[sort_idx]), color = color[t], label = 'depth = %.1f' %(t + 1)) plt.xlabel('LSTAT') plt.ylabel('MEDV') plt.legend(loc = 'upper right') plt.show()
図から、木の深さが3より小さいとデータの傾向を表現しきれておらず、逆に3より大きくなると予測値が非常に複雑な振る舞いをしてしまう(過学習する)ことがわかります。今の場合は、木の深さ3がバランスがよいと判断されるでしょう。
Next Step
今回は前回までと毛色の異なった決定木の手法を紹介しました。ざっくりデータの傾向を説明できる利点がある一方、過学習の危険性を内包した手法でした。さらに、ある工夫をするより性能が高いモデルを構築できるようになるのですが、それは次回以降に紹介していこうと思います*1。
*1:というよりそれらの手法をやりたかったからあえてまず決定木を紹介したのでした。
Pythonでデータ分析:非線形効果を導入
導入
前回、ボストン近郊の住宅情報のデータを用いて線形回帰モデルを作りました。
tekenuko.hatenablog.com
今回は、モデルの性能を上げる可能性の一つとして、多項式や指数・対数などの非線形効果をモデルに投入した場合の振る舞いを見ようと思います。
使用するデータセット
前回と同様、ボストン近郊の住宅情報のデータを使用します。
# 必要なライブラリのインポート import numpy as np import pandas as pd from pandas import DataFrame from sklearn.datasets import load_boston # データのロード、マージ boston = load_boston() df = DataFrame(boston.data, columns = boston.feature_names) df['MEDV'] = np.array(boston.target)
非線形効果を入れる前に
非線形効果を入れる際に、(私個人の考えとして)注意する点が2つあります。
- 過学習
- 投入コスト
過学習に関しては、モデルの複雑さが増すことから生じます。一般に、説明変数と目的変数のデータがそれぞれn個あったとすると、n-1次の多項式を用いると残差が0で説明できてしまいます。しかしながら、このようなモデルは今のデータセットのみに強く適合するものであり、新しいデータに対する説明力が極端に低くなってしまいます。そのため、常に学習データと検証データにおける精度のバランスを見ながらモデルを評価する必要があります。
投入コストに関しては、ビジネスでのデータ分析の場面で度々考慮しなければならない問題です。実務では、納期と精度、そして費用対効果のバランスをうまくとることが求められます。非常に時間をかけて(ほんのすこし精度が向上する)非線形効果を取り入れたモデルを構築するよりも、時間をかけずに単純な(ある程度精度が担保された)モデルを作成し、早く実運用に載せて知見をためるという方が案件状況によっては良い場合があります。そのため、モデルの複雑化に関しては、要件定義や分析設計の段階でプロジェクト内での認識合わせが非常に重要になります*1。
このような注意点があるため、闇雲に変数を複雑化することはできるだけ避けたほうがよいというのが私個人の考えです。検討する場合も、事前知識を援用して判断できるもののみを考え、複雑化は最小限に留める、などします*2。
可視化
事前知識を援用する方法の一つに、データを可視化するという方法があります(あると思っています)。以下のように、変数間でどのような関係が見られるかをプロットしてみます。
# matplotlib(のなかのpyplot)をインポート import matplotlib.pyplot as plt # seaborn(可視化の補助を行う)のインポート import seaborn as sns # Jupyterを利用していたら、以下のおまじないを書くとnotebook上に図が表示 %matplotlib inline # グラフのスタイルを指定 sns.set(style = 'whitegrid', context = 'notebook') # 変数のペアの関係をプロット sns.pairplot(df, size = 2.5) plt.show()
今回は、説明変数と目的変数の関係、つまり図の一番下(もしくは一番右)に着目します。
とはいえ、これでは変数が多すぎて情報過多になっています(汗)。そのため、今回フォーカスしたい部分に制限した図を以下に表示します。
# カラムを制限してプロット sns.pairplot(df[['LSTAT', 'MEDV']], size = 2.5) plt.show()
このLSTATとMEDVの散布図を見ると、直線的というよりは若干曲線的な関係性がありそうです。そのため、LSTATに関して非線形効果を導入することがモデルの性能向上に役に立つかもしれません。今回は、単純化したい理由で、説明変数はLSTATのみに絞り、モデルに投入する式を変えた場合の精度の変化を見たいと思います。
非線形効果 : 多項式
線形(1次)からの素朴な拡張として、多項式効果を入れる、という方法があります。ここでは、2次(quadratic)と3次(cubic)の項をモデルに投入した場合の精度を評価してみます。
多項式項を追加する場合、モデル構築の大枠は変更しなくても大丈夫です。
# 説明変数 X = df.loc[:, ['LSTAT']].values # 目的変数 y = df.loc[:, 'MEDV'].values # モデルのインスタンス生成 from sklearn.linear_model import LinearRegression mod = LinearRegression()
変更を受けるのは、変数加工の部分です。つまり、自分で2次、3次の変数を作成し、それを代入する、ということが実際にやることです。幸い(?)、scikit-learnには非線形の変数を加工する機能がありますので、ここではそれを利用します。
多項式の変数加工は、以下のメソッドを呼び出すことで行えます。
from sklearn.preprocessing import PolynomialFeatures # 2次(までの)変数を作成するインスタンス quadratic = PolynomialFeatures(degree = 2) # 3次(までの)変数を作成するインスタンス cubic = PolynomialFeatures(degree = 3) # 変数作成 X_quad = quadratic.fit_transform(X) X_cubic = cubic.fit_transform(X)
これで多項式変数の作成ができました。
次に、線形回帰モデル、2次の項を追加、3次の項まで追加したモデルを作成します。
# モデル式用に変数を作成 X_fit = np.arange(X.min(), X.max(), 1)[:, np.newaxis] # 線形回帰モデル、予測値、R^2を評価 mod_lin = mod.fit(X, y) y_lin_fit = mod_lin.predict(X_fit) r2_lin = mod.score(X, y) # 2次の項を追加、予測値、R^2を評価 mod_quad = mod.fit(X_quad, y) y_quad_fit = mod_quad.predict(quadratic.fit_transform(X_fit)) r2_quad = mod.score(X_quad, y) # 3次の項を追加、予測値、R^2を評価 mod_cubic = mod.fit(X_cubic, y) y_cubic_fit = mod_cubic.predict(cubic.fit_transform(X_fit)) r2_cubic = mod.score(X_cubic, y)
各モデルの結果をプロットします。参考のため、モデル式も重ねてプロットします。
# データ点をプロット plt.scatter(X, y, label = 'Traning points', color = 'lightgray') # 線形モデルのモデル式 plt.plot(X_fit, y_lin_fit, label = 'linear (d = 1), $R^2=%.2f$' % r2_lin, color = 'blue', lw = 2, linestyle = ':') # 2次 plt.plot(X_fit, y_quad_fit, label = 'quadratic (d = 2), $R^2=%.2f$' % r2_quad, color = 'red', lw = 2, linestyle = '-') # 3次 plt.plot(X_fit, y_cubic_fit, label = 'cubic (d = 3), $R^2=%.2f$' % r2_cubic, color = 'green', lw = 2, linestyle = '--') plt.xlabel('LSTAT') plt.ylabel('MEDV') plt.legend(loc = 'upper right') plt.show()
結果を見ると、3次の項まで追加した場合が最も当てはまりのよいモデルになっていることがわかります。
学習用データと検証用データで分割
前回の精度検証と同様、今回も学習用データと検証用データに分割して性能を判断してみます。
# 必要なメソッドのインポート from sklearn.cross_validation import train_test_split # 学習用70%、検証用30%に分割 (X_train, X_test, y_train, y_test) = train_test_split(X, y ,test_size = 0.3, random_state = 666) (X_quad_train, X_quad_test, y_train, y_test) = train_test_split(X_quad, y , test_size = 0.3, random_state = 666) (X_cubic_train, X_cubic_test, y_train, y_test) = train_test_split(X_cubic, y , test_size = 0.3, random_state = 666)
それぞれのモデルでのMSEとは以下になります。
線形モデル
# モデル作成 mod.fit(X_train, y_train) y_train_pred = mod.predict(X_train) y_test_pred = mod.predict(X_test) # 必要なメソッドのインポート 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を出力 print('R^2 Train : %.3f, Test : %.3f' % (mod.score(X_train, y_train), mod.score(X_test, y_test))) # 出力 >>>print('MSE Train : %.3f, Test : %.3f' % (mean_squared_error(y_train, y_train_pred), mean_squared_error(y_test, y_test_pred))) MSE Train : 42.126, Test : 30.120 >>>print('R^2 Train : %.3f, Test : %.3f' % (mod.score(X_train, y_train), mod.score(X_test, y_test))) R^2 Train : 0.532, Test : 0.578
2次の項を追加
# モデル作成 mod.fit(X_train, y_train) y_train_pred = mod.predict(X_quad_train) y_test_pred = mod.predict(X_quad_test) # 必要なメソッドのインポート 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を出力 print('R^2 Train : %.3f, Test : %.3f' % (mod.score(X_quad_train, y_train), mod.score(X_quad_test, y_test))) # 出力 >>>print('MSE Train : %.3f, Test : %.3f' % (mean_squared_error(y_train, y_train_pred), mean_squared_error(y_test, y_test_pred))) MSE Train : 32.211, Test : 26.241 >>>print('R^2 Train : %.3f, Test : %.3f' % (mod.score(X_quad_train, y_train), mod.score(X_quad_test, y_test))) R^2 Train : 0.642, Test : 0.632
3次の項を追加
# モデル作成 mod.fit(X_train, y_train) y_train_pred = mod.predict(X_cubic_train) y_test_pred = mod.predict(X_cubic_test) # 必要なメソッドのインポート 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を出力 print('R^2 Train : %.3f, Test : %.3f' % (mod.score(X_cubic_train, y_train), mod.score(X_cubic_test, y_test))) # 出力 >>>print('MSE Train : %.3f, Test : %.3f' % (mean_squared_error(y_train, y_train_pred), mean_squared_error(y_test, y_test_pred))) MSE Train : 29.914, Test : 26.968 >>>print('R^2 Train : %.3f, Test : %.3f' % (mod.score(X_cubic_train, y_train), mod.score(X_cubic_test, y_test))) R^2 Train : 0.668, Test : 0.622
検証用データも含めて判断すると
という傾向があるように見えます。
他の非線形効果
多項式以外の投入方法として、指数、対数、平方根をとるといった方法があります。今回のデータだと、LSTATの対数をとり、MEDVの平方根をとる、という操作を行うとよりモデルの当てはまりがより向上します。結果をプロットすると以下のようになります。
# 非線形変換 X_log = np.log(X) y_sqrt = np.sqrt(y) # 線形回帰モデル作成、R^2計算 X_fit = np.arange(X_log.min(), X_log.max() + 1, 1)[:, np.newaxis] mod_log = mod.fit(X_log, y_sqrt) y_sqrt_fit = mod_log.predict(X_fit) r2_sqrt = mod.score(X_log, y_sqrt) # プロット plt.scatter(X_log, y_sqrt, label = 'Traning points', color = 'lightgray') plt.plot(X_fit, y_sqrt_fit, label = 'linear (d = 1), $R^2=%.2f$' % r2_sqrt, color = 'blue', lw = 2, linestyle = ':') plt.xlabel('log[LSTAT]') plt.ylabel('sqrt[MEDV]') plt.legend(loc = 'upper right') plt.show()
となっており、今回考慮したモデルの中では最も高い数値になっています。
まとめ
非線形効果を入れると、一般には当てはまりがよくなるが、過学習しやすくなる傾向があります。
過学習が強いのは、予測を目的とすることが多い実社会の分析では、大きなデメリットです。早い段階で線形回帰モデル程度の精度で費用対効果的に十分である、といった見立てができているのならば、無理に非線形効果を入れずに線形モデルで施策を進めるというのも手だと思います。
*1:例えば、ビジネス要件を鑑みて、プロトタイプモデルでどれくらいの精度ならモデルの磨き込みはこれくらいの期間で、最終的にはこれくらいの精度を目指す、運用に載せてどうだったらどのようにモデルを改善する、とかを議論しておくなどです。ざっくりいうと、「正解」を探すのではなく、「納得解」を見つけるためにどうするかを決めておく感じだと思います。
*2:ここはあくまで私個人の意見です。
*3:サンプリングの影響は、外れ値があることが影響しているかもしれません。外れ値の除去や、外れ値の影響を抑えるモデル構築といった方向もモデルの性能向上方針の一つです。というより、実際に自分だったらまず回帰分析前の基礎集計の時点で外れ値の対応を検討するだろうなと思います。