読者です 読者をやめる 読者になる 読者になる

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

元素粒子理論屋。いっぱしのデータ分析屋になるために修行中です。

FastBDTの計算時間が速いかを確認してみる

はじめに

最近、ブースティング系のアルゴリズムでXGboostより速いものが実装されているようです。
github.com


論文は以下になります。
[1609.06119] FastBDT: A speed-optimized and cache-friendly implementation of stochastic gradient-boosted decision trees for multivariate classification

論文を見ると、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つのタウ粒子という、電子の親戚みたいな素粒子に崩壊する現象(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とは意味が異なります。この辺は、説明するととても長くなるのは、このへんでお茶をにごしておきます(汗)。

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

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()

f:id:tekenuko:20160919205618p:plain

線形回帰や非線形回帰と比較して、カクカクした直線が引かれます。その「カクカクする」性質のため、ある程度データの傾向を補足できるのが決定木のメリットです。一方で、予測値が微分不可能になってしまう制約があるため、解析的には扱いづらい手法になっています。

木の深さと予測値の関係

先程は、木の深さを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()

f:id:tekenuko:20160919205536p:plain

図から、木の深さが3より小さいとデータの傾向を表現しきれておらず、逆に3より大きくなると予測値が非常に複雑な振る舞いをしてしまう(過学習する)ことがわかります。今の場合は、木の深さ3がバランスがよいと判断されるでしょう。

Next Step

今回は前回までと毛色の異なった決定木の手法を紹介しました。ざっくりデータの傾向を説明できる利点がある一方、過学習の危険性を内包した手法でした。さらに、ある工夫をするより性能が高いモデルを構築できるようになるのですが、それは次回以降に紹介していこうと思います*1

*1:というよりそれらの手法をやりたかったからあえてまず決定木を紹介したのでした。

Pythonでデータ分析:非線形効果を導入

導入

前回、ボストン近郊の住宅情報のデータを用いて線形回帰モデルを作りました。
tekenuko.hatenablog.com

今回は、モデルの性能を上げる可能性の一つとして、多項式や指数・対数などの非線形効果をモデルに投入した場合の振る舞いを見ようと思います。

参考

前回と同様、以下の書籍を参考(写経?)にします。

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)

非線形効果を入れる前に

非線形効果を入れる際に、(私個人の考えとして)注意する点が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()

f:id:tekenuko:20160919162755p:plain

今回は、説明変数と目的変数の関係、つまり図の一番下(もしくは一番右)に着目します。
とはいえ、これでは変数が多すぎて情報過多になっています(汗)。そのため、今回フォーカスしたい部分に制限した図を以下に表示します。

# カラムを制限してプロット
sns.pairplot(df[['LSTAT', 'MEDV']], size = 2.5)
plt.show()

f:id:tekenuko:20160919163312p:plain

この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()

f:id:tekenuko:20160919181952p:plain

結果を見ると、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と R^2は以下になります。

線形モデル
# モデル作成
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

検証用データも含めて判断すると

  • 線形モデルは最も汎化性能は高い(?)が、当てはまりが悪い
  • 3次の項を追加すると過学習が高まるが、当てはまりは良い
  • 汎化性能と過学習の観点では、この中では2次のモデルが一番バランスが良い(?)

という傾向があるように見えます。

結局どうなのか

さらに調べるとわかるのですが、先程の結論(的なもの)には落とし穴があります。実は、先程の結果は、サンプリングにかなり依存しており(つまりseed値を変えたときの学習用データと検証用データのR^2の差の変化が激しい)、2次の項を加えたモデルは学習用データでのR^2のほうが検証用データでのR^2より大きいといった傾向が見えます。逆に、線形モデルはサンプリングを変えても多くの場合、過学習が弱い傾向があります。というわけで、サンプリングも含めると以下の結論になると思われます*3

  • 線形モデルは最も汎化性能は高いが、当てはまりが悪い
  • 非線形項(2次、3次)を追加すると過学習が強くなるが、当てはまりは良くなる

他の非線形効果

多項式以外の投入方法として、指数、対数平方根をとるといった方法があります。今回のデータだと、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()

f:id:tekenuko:20160919193235p:plain

R^2 = 0.69となっており、今回考慮したモデルの中では最も高い数値になっています。

まとめ

非線形効果を入れると、一般には当てはまりがよくなるが、過学習しやすくなる傾向があります。
過学習が強いのは、予測を目的とすることが多い実社会の分析では、大きなデメリットです。早い段階で線形回帰モデル程度の精度で費用対効果的に十分である、といった見立てができているのならば、無理に非線形効果を入れずに線形モデルで施策を進めるというのも手だと思います。

*1:例えば、ビジネス要件を鑑みて、プロトタイプモデルでどれくらいの精度ならモデルの磨き込みはこれくらいの期間で、最終的にはこれくらいの精度を目指す、運用に載せてどうだったらどのようにモデルを改善する、とかを議論しておくなどです。ざっくりいうと、「正解」を探すのではなく、「納得解」を見つけるためにどうするかを決めておく感じだと思います。

*2:ここはあくまで私個人の意見です。

*3:サンプリングの影響は、外れ値があることが影響しているかもしれません。外れ値の除去や、外れ値の影響を抑えるモデル構築といった方向もモデルの性能向上方針の一つです。というより、実際に自分だったらまず回帰分析前の基礎集計の時点で外れ値の対応を検討するだろうなと思います。

Pythonでデータ分析:線形回帰モデル

導入

データ分析にて、最も基本的な回帰分析から始めていきます*1。回帰分析とは、説明したい変数(目的変数)とそれを説明するための変数(説明変数)の間の関係を求める手法です。機械学習の手法の区分としては、教師あり学習(解答に相当する教師データを用いてモデルを構築)に属するものです。

回帰分析の目的は大きく2種類あります。

  • 理解:文字通り、変数間の構造・関係の理解
    • 例:ある病気に対して、どのような遺伝子が寄与しているかをマイクロアレイデータなどのデータを用いて明らかにする
  • 予測:変数間の関係をもとに、将来予測をする
    • 例:過去の広告の投下量と売上実績から、これらの関係を説明するモデルを構築し、それをもとに将来の広告への投資計画を立てる

ビジネスの世界では、予測を目的に使用されることがほとんどです。それは、モデル構築によって要因を明らかにするだけでなく、その後の施策でそのモデルの知見を使い、価値を生み出すことが求められるからです。逆に、アカデミックでは理解が目的の研究も一定数あります。実務でデータ分析をやっている方は、例えば、マーケティングサイエンス学会に参加するとこのあたりの思想の違いを感じ取れて面白いと思います。

使用するデータセット

ボストン近郊の住宅情報のデータを使用します。こちらは、無償で提供されており、例えばUCI Machine Learning Repository(UCI Machine Learning Repository: Housing Data Set)からダウンロードすることが可能です。ここでは、scikit-learnからロードする形で導入します。

## 必要なライブラリのインポート
from pandas import DataFrame
from sklearn.datasets import load_boston

## データの用意
# ボストンデータをインポート
boston = load_boston()
# 説明変数たちをDataFrameへ変換
df = DataFrame(boston.data, columns = boston.feature_names)
# 目的変数をDataFrameへ追加
df['MEDV'] = np.array(boston.target)

DataFrameという形式で保持するために、pandasを導入しています。最初の数行を表示してみます。

# 最初の5行を表示
df.head()

# コンソールでの表示結果
>>> df.head()
      CRIM    ZN  INDUS  CHAS    NOX     RM   AGE     DIS  RAD    TAX  
0  0.00632  18.0   2.31   0.0  0.538  6.575  65.2  4.0900  1.0  296.0   
1  0.02731   0.0   7.07   0.0  0.469  6.421  78.9  4.9671  2.0  242.0   
2  0.02729   0.0   7.07   0.0  0.469  7.185  61.1  4.9671  2.0  242.0   
3  0.03237   0.0   2.18   0.0  0.458  6.998  45.8  6.0622  3.0  222.0   
4  0.06905   0.0   2.18   0.0  0.458  7.147  54.2  6.0622  3.0  222.0   

   PTRATIO       B  LSTAT  MEDV  
0     15.3  396.90   4.98  24.0  
1     17.8  396.90   9.14  21.6  
2     17.8  392.83   4.03  34.7  
3     18.7  394.63   2.94  33.4  
4     18.7  396.90   5.33  36.2 

全部で14個のカラムがありますが、今回は最初の13個を説明変数、14個目を目的変数として使用することを想定しています。それぞれのカラムにある情報は以下です。

  • CRIM : 犯罪発生率(人口単位)
  • ZN : 25000平方フィート以上の住宅区間の割合
  • INDUS : 非小売業の土地面積の割合(人口単位)
  • CHAS : チャールズ川沿いかどうかフラグ(川沿いなら1、そうでなければ0)
  • NOX : 窒素感化物の濃度(pphm単位)
  • RM : 一戸あたりの平均部屋数
  • AGE : 1940年よりも前に建てられた家屋の割合
  • DIS : ボストンの主な5つの雇用圏までの重み付き距離
  • RAD : 幹線道路へのアクセス指数
  • TAX : 10000ドルあたりの所得税
  • PTRATIO : 教師あたりの生徒の数(人口単位)
  • B : 1000(Bk - 0.63)^2として計算される量(Bkはアフリカ系アメリカ人居住者の割合(人口単位))
  • LSTAT : 低所得者の割合
  • MEDV : 住宅価格の中央値(単位は1000ドル)

線形回帰モデルの構築

通常はモデル構築に至るまでにデータ理解や基礎集計、データ加工といった非常に重要なステップがありますが、省略します。ここでは、目的変数をMEDV、説明変数をCRIM~LSTATとした回帰モデルを構築します。今回、目的変数は連続値とみなすことができるため、回帰モデルの中でも最も単純な、目的変数と説明変数間の関係が線形である場合を仮定します。
 MEDV \sim \beta _0 + \beta _1 \times ZN + \cdots + \beta _{13} \times LSTAT + \epsilon
\beta _0, \cdots ,\beta _{13}は回帰係数と呼ばれる量で、それぞれの変数がどれだけ目的変数に寄与するかを表しています。特に、\beta _0はすべての説明変数の寄与が0の場合も残る「ベース」的な量で、切片(Intercept)とも呼ばれます。 \epsilonは残差と呼ばれる量で、各データとモデル式の「ズレ」を表します。線形回帰モデルでは、各データ点の残差の二乗を足し合わせたものを最小にするという条件のもとで、係数を求めます。

scikit-learnによる係数の推定

ここでは、いったん全データを用いて線形回帰モデルを構築(=回帰係数の推定)をします。scikit-learnにあるLinearRegressionを用いると、簡単に係数の推定が行えます。

## 説明変数と目的変数
# 説明変数
X = df.loc[:, boston.feature_names].values
# 目的変数
y = df.loc[:, 'MEDV'].values
## 必要なライブラリのインポート
from sklearn.linear_model import LinearRegression
# オブジェクト生成
mod = LinearRegression(fit_intercept = True, normalize = False, copy_X = True, n_jobs = 1)
# fit関数でパラメータ推定
mod.fit(X, y)

LinearRegressionの引数を明示的に書いてありますが、それぞれは以下の意味を持っています。

  • fit_intercept : 切片項をモデル式に入れるかどうか。デフォルト値はTrue。
  • normalize : 説明変数を正規化するかどうか。デフォルト値はFalse。
  • copy_X : メモリ内でデータを複製してから実行するかどうか。 デフォルト値はTrue。
  • n_jobs : 計算に使うジョブ数。デフォルト値は1で-1にするとすべてのCPUを使って計算。

回帰係数を出力してみます。

# 回帰係数を出力
print(mod.coef_)
print(mod.intercept_)

# 出力結果
>>>print(mod.coef_)
[ -1.07170557e-01   4.63952195e-02   2.08602395e-02   2.68856140e+00
  -1.77957587e+01   3.80475246e+00   7.51061703e-04  -1.47575880e+00
   3.05655038e-01  -1.23293463e-02  -9.53463555e-01   9.39251272e-03
  -5.25466633e-01]
>>>print(mod.intercept_)
36.4911032804

学習データ、検証データに分割して性能評価

先程は、全データを用いて回帰係数の推定を行いました。しかし、それだけだと今持っているデータ「だけ」によく当てはまるモデルを作ってしまい、予測で良い性能を発揮しない場合があります。それを確認するための方法として、パラメータ推定の際に使用しなかったデータを用いて精度検証を行う、という方法があります。ここでは、持っているデータを学習用、検証用の2種類に分割して精度検証します。

上記のような精度検証は、scikit-learnを利用すると簡単に行うことができます。

# train_test_splitをインポート
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)
# 学習用データでパラメータ推定
mod.fit(X_train, y_train)
# 作成したモデルから予測(学習用、検証用モデル使用)
y_train_pred = mod.predict(X_train)
y_test_pred = mod.predict(X_test)

残差プロット

モデルの性能を見る大雑把な方法として、残差プロットがあります。これは、予測値と実際の観測値の差分を可視化し、モデルで拾いきれていない非線形性や外れ値がないか、といった観点から、残差がランダムに分布しているかをチェックする方法です。数値的に性能を評価する前に、可視化して傾向を見ておくのは重要です。残差をプロットするコードは、例えば以下になります。

# matplotlibをインポート
import matplotlib.pyplot as plt
# Jupyterを利用していたら、以下のおまじないを書くとnotebook上に図が表示
%matplotlib inline
# 学習用、検証用それぞれで残差をプロット
plt.scatter(y_train_pred, y_train_pred - y_train, c = 'blue', marker = 'o', label = 'Train Data')
plt.scatter(y_test_pred, y_test_pred - y_test, c = 'lightgreen', marker = 's', label = 'Test Data')
plt.xlabel('Predicted Values')
plt.ylabel('Residuals')
# 凡例を左上に表示
plt.legend(loc = 'upper left')
# y = 0に直線を引く
plt.hlines(y = 0, xmin = -10, xmax = 50, lw = 2, color = 'red')
plt.xlim([10, 50])
plt.show()

コードの実行結果は以下のようになります。
f:id:tekenuko:20160919144249p:plain

予測が完璧である場合は、残差はちょうど0になりますが、現実的にはそれはまずありません。良い回帰モデルの場合は、残差がちょうと0を中心にランダムに分布するような傾向が見えるはずです。今回のモデルに関しては、概ねランダムに分布しているようにも見えますが、学習用データを使った場合の残差に関してグラフの右下になにか直線的な傾向があり、モデルに取り切れていない情報がある可能性があります。

平均二乗誤差、決定係数

数値的にモデルの性能を評価する指標の例として、平均二乗誤差(Mean Squared Error : MSE)や決定係数( R^2)があります。これらは、データ点を iでラベル付けしたとすると、以下のように定義されるものです。
 MSE = \frac{1}{n} \sum _{i = 1}^n (y^{(i)} - y_{pred}^{(i)})^2
 R^2 = 1 - \frac{MSE}{Var(y)}
MSEは、誤差の二乗和の平均値で、学習用データ、検証用データともに小さければモデルの性能が良い、と判断されます。R^2に関しては、MSEが0の場合に1をとる、つまりモデルの性能が良いほど1に近い量になります。また、R^2は学習用データに関しては0から1の値を取りますが、検証用データの場合には負の値を取ることもあります。

MSEとR^2は、scikit-learnを用いると簡単に評価できます。これらを出力するコードは、以下のようになります。

# 平均二乗誤差を評価するためのメソッドを呼び出し
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 : 21.371, Test : 24.336
>>>print('R^2 Train : %.3f, Test : %.3f' % (mod.score(X_train, y_train), mod.score(X_test, y_test)))
R^2 Train : 0.763, Test : 0.659

ちなみに、R^2はsklearn.metrics.r2_scoreを呼び出して評価することも可能です。結果を見ると、学習用データのほうが検証用データよりも

  • MSEが小さい
  • R^2が大きい

という傾向があり、過学習気味のモデルであることがわかります。

モデルの磨き込み

上記の結果から、モデルの性能を上げるために、大きく2つの可能性があると考えられます。

  • モデルに非線形性を導入
  • 変数を制限する(=変数選択)

一つ目は、モデルの表現力を底上げする方向、二つ目は、過学習を緩和される方向です。次回以降は、これらの観点を取り入れて磨き込みをしていった場合にどうなるかを検討していきます。

*1:これは、分析プロジェクトでいうモデル構築のフェーズでの基本という意味です。

Pythonでデータ分析:導入

背景

最近、業務でPythonを使っているのですが、不慣れな部分もありRレベルで自在に使いこなせていないと感じています*1。そのため、基本的な部分からおさらいをしていこうと考えました。

目標

試してみたいデータがあった場合に、簡単な分析に関してはとっさに手が動くレベルをまず目指します。そのために、データ加工ではnumpyとpandas(あとseries)、モデル構築ではscikit-learn、可視化ではmatplotlibとseabornを交えて紹介していくことで、この目標を実現していこうと考えています*2

進め方

データ分析をやる上でもっとも基本的と思われる回帰分析から始めていこうと思います。基本的なことがらではありますが

  • 様々な拡張的手法のベースになっている
  • ビジネスの世界では解釈のしやすさなどの理由で回帰分析がまだまだ大事

といった理由から、このアウトプットは無駄にはならないと期待しています。

線形回帰のあとは、回帰問題だけではなく、機械学習の手法を幾つか取り上げてまとめていきます。また、単なるライブラリの使用だけでなく、可能な範囲で実装もできたらよいなと考えています。

参考文献

以下の本では、ライブラリの使い方だけでなく、実装も交えて機械学習の手法が紹介されているので、こちらを参考にしていこうと考えています*3
book.impress.co.jp

*1:定期的にPythonは勉強しているのですが、業務ではRを使用する率が圧倒的に高かったり、携わるのが上流工程で、実分析は別担当者が…(涙)といった案件もたびたびあり、いざPythonを使おうとなったときに結構忘れてたりします。

*2:あえてブログでアウトプットすることで、記憶の定着度の向上も狙いです

*3:あくまでトレーニングが目的なので、単なる写経になってしまう可能性もありますが、ご了承ください。