データサイエンティスト(仮)

元素粒子論博士。今はデータサイエンティスト(仮)。

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)とR^2がどうなるかを見てみます。

# 予測値を計算
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

今まで考えてきた手法と比較すると、かなり性能が上がっています。残念ながら、過学習する傾向は残っています。それでも検証用データに対するR^2は約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()

f:id:tekenuko:20160920220004p:plain

図を見ると、検証データよりも学習データに関するほうが散らばりは少なく、実際に学習データにより適合したモデルであることがわかります。さらに、学習データに関して、右下にランダムに見えない振る舞いが見られ、ランダムフォレストでも完全に捉えきれていない情報がある可能性があります。

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

さらにモデルの性能を上げる方法として、ハイパーパラメータの調整があります。ハイパーパラメータとは、モデルに事前に手で投入するパラメータのことで、これを変化させると精度が上がる可能性があります。ハイパーパラメータの調整方法は、大雑把に以下の方法があります。

  • グリッドサーチ:決められた範囲を総当り探索
  • ランダムサーチ:投入するパラメータをランダムに変えて試す

今回は、グリッドサーチを試してみます。グリッドサーチは、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とR^2は以下のようになります。

# 予測値を計算
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、そこをさらに細かく探す、といった方法を取るのがよいでしょう。

*1:条件式の組み合わせでモデルがされるので解釈しやすい、SQLでモデルが表現できるため、施策に落とすときに使いやすいといった利点もあります。

*2:このあたりは、バイアス-バリアンス分解とアンサンブル学習の一つであるバギングを知っているとわかる性質です。

*3:潤沢に時間があったらわかりやすく紹介するとかしたいですね。

*4:ここは手でやったり、ベイズ的最適化でデータからあたりをつけるといったアプローチが考えられます。