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

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

複数の棒グラフを表示させるのはpandasが便利

経緯

ある対象に対して、複数のアプローチの結果を可視化したいとき、棒グラフで並べて比較する方法があります。これをmatplotlib.pyplot.bar()で描いていましたが、棒の太さやら目盛の調整が大変でした。matplotlibは柔軟な可視化ができる反面、匠の技が要求されることもあるため、使いこなすのは結構難しいです。より簡単にできる方法はないかと考えていたところ、PandasのDataFrameにmatplotlibのラッパーがあったことを思い出しました。というわけで、復習がてら載せておきます。

可視化する対象

幾つかのモデルの回帰係数の値を可視化します。使用データは、ボストン近郊の住宅情報、モデルは線形回帰、Lasso、Ridge、Elastic Netの4種類です。Lasso、Ridge、Elastic Netは罰則付き回帰と呼ばれる手法で、目的関数に罰則項というモデルの複雑さを調整する項を入れることで、過学習を緩和させることが目的の一つです。手法によって、特定の重要でない変数に対する回帰係数が0になったり、係数が全体的に縮小する、などバリエーションがあります。今回は、代表的な3手法の推定結果にどのような違いが出るかを可視化してみます。

可視化

ライブラリのインポート

必要なライブラリたちをインポートしておきます。matplotlibを呼び出していますが、これはグラフの整形に使っています。

# データ加工
import numpy as np
from pandas import DataFrame
# サンプルデータ
from sklearn.datasets import load_boston
# 学習用、検証用データに分割
from sklearn.cross_validation import train_test_split
# 比較対象:線形回帰
from sklearn.linear_model import LinearRegression
# 罰則付き回帰、正則化項をクロスバリデーションで求めるライブラリを使用
from sklearn.linear_model import LassoCV
from sklearn.linear_model import RidgeCV
from sklearn.linear_model import ElasticNetCV
# MSE、R^2
from sklearn.metrics import mean_squared_error
from sklearn.metrics import r2_score
# 可視化
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
# きれいに書くおまじない
plt.style.use('ggplot')
font = {'family': 'meiryo'}
matplotlib.rc('font', **font)

サンプルデータ

サンプルデータは、ボストン近郊の住宅情報です。そろそろ使いすぎてデータマートを作る操作が滑らかになってきました。

# データのロード、データマート作成
boston = load_boston()
df = DataFrame(boston.data, columns = boston.feature_names)
df['MEDV'] = boston.target
# 説明変数、目的変数
X = df.iloc[:, :-1].values
y = df.loc[:, 'MEDV'].values
# 学習用、検証用データに分割
(X_train, X_test, y_train, y_test) = train_test_split(X, y, test_size = 0.3,  random_state = 666)

モデル作成、プロット

matploblibの場合、プロットしたい対象を一つずつ追加していくという方法が取られます。PandasのDataFrameのplotメソッドは、DataFrameのデータを可視化することができます。棒グラフは、plot(kind = 'bar')と指定するかplot.bar()とすると指定できます。DataFrameに対してこのメソッドを使用すると、indexが横軸で、columnsごとにデータを表示します。今回、4つのモデルの係数に関する棒グラフを表示したいので、DataFrameにはindexには各変数を、columnsには各モデルをラベルするようにデータを格納します。

コード例と結果は以下となります。

# モデルインスタンス生成
Linear = LinearRegression()
Lasso = LassoCV()
Ridge = RidgeCV()
ElasticNet = ElasticNetCV()
# 学習
Linear.fit(X_train, y_train)
Lasso.fit(X_train, y_train)
Ridge.fit(X_train, y_train)
ElasticNet.fit(X_train, y_train)
# 予測値
y_train_Linear_pred = Linear.predict(X_train)
y_test_Linear_pred = Linear.predict(X_test)
y_train_Lasso_pred = Lasso.predict(X_train)
y_test_Lasso_pred = Lasso.predict(X_test)
y_train_Ridge_pred = Ridge.predict(X_train)
y_test_Ridge_pred = Ridge.predict(X_test)
y_train_ElasticNet_pred = ElasticNet.predict(X_train)
y_test_ElasticNet_pred = ElasticNet.predict(X_test)
# データフレームに格納
df_coef = DataFrame(data = np.transpose(np.array([Linear.coef_, Lasso.coef_, Ridge.coef_, ElasticNet.coef_])),  
                    index = df.iloc[:, :-1].columns, 
                    columns = ['Linear MSE train: %.3f, test: %.3f, R^2 train: %.3f, test: %.3f' 
                               %(mean_squared_error(y_train, y_train_Linear_pred),  mean_squared_error(y_test, y_test_Linear_pred),
                                 r2_score(y_train, y_train_Linear_pred), r2_score(y_test, y_test_Linear_pred)),
                               'Lasso MSE train: %.3f, test: %.3f, R^2 train: %.3f, test: %.3f' 
                               %(mean_squared_error(y_train, y_train_Lasso_pred),  mean_squared_error(y_test, y_test_Lasso_pred),
                                 r2_score(y_train, y_train_Lasso_pred), r2_score(y_test, y_test_Lasso_pred)), 
                               'Ridge MSE train: %.3f, test: %.3f, R^2 train: %.3f, test: %.3f' 
                               %(mean_squared_error(y_train, y_train_Ridge_pred),  mean_squared_error(y_test, y_test_Ridge_pred),
                                 r2_score(y_train, y_train_Ridge_pred), r2_score(y_test, y_test_Ridge_pred)), 
                               'ElasticNet MSE train: %.3f, test: %.3f, R^2 train: %.3f, test: %.3f' 
                               %(mean_squared_error(y_train, y_train_ElasticNet_pred),  mean_squared_error(y_test, y_test_ElasticNet_pred),
                                 r2_score(y_train, y_train_ElasticNet_pred), r2_score(y_test, y_test_ElasticNet_pred)) ])
# プロット
df_coef.plot.bar(figsize = (15, 8), width = 1)
plt.title('罰則付き回帰で推定された回帰係数の比較', size = 20)
plt.xlabel('特徴量')
plt.ylabel('回帰係数の大きさ')
plt.legend(loc = 'best')

f:id:tekenuko:20160924214300p:plain

いい感じに横に並べることができました。ちなみに、棒の太さはwidthでいじることができ、デフォルトは0.8です。これを大きくしすぎると、棒グラフどうしがかぶってしまって何のことやら…という状態になるので、振る舞いを見てちょうどよい値を選択すると良いと思います。

回帰係数の結果ですが、Lassoは重要でない変数が0になる、Ridgeは線形回帰と少し似ているが、全体的に係数のサイズが小さめ、Elastic Netはそれらの折衷案、のような振る舞いをしています。また、MSEとR^2を見ると、LassoとElastic Netは学習データの精度を少し犠牲にして汎化性能が向上しているようには見えます。ただし、全体的に線形回帰と比較して精度が落ちているため、これらの間で何を採用するかは悩ましいところです。

他のプロット

pandasのplotだけでも工夫するといろいろといじれます。例えば、以下が参考になると思います。

sinhrks.hatenablog.com


matplotlibを使いこなして職人になるのもよいですが、時間の制約があった場合はこういった便利なラッパーを使うとよい場面があるかもしれません*1

*1:自分は割りと原理主義者なのでmatplotlibを調整できるようになりたいですが(笑)。

アソシエーション分析のちょい復習

動機

最近こんな本を購入しました。
www.shoeisha.co.jp
第I部の「データアナリティクスの基礎」は、データ分析のプロジェクトの全体像に触れていて、データ分析のプロジェクトをどう進めていくかについて整理するのには良いのかな、と思いました*1。後半部分の分析内容は結構あっさりしていて、すでに専門的なスキルを持っている方向けではなさそうです。

この書籍の第8章で、Rを使ったアソシエーション分析を取り扱っています。これは、マーケティングの分野でよく使用される手法です。顧客の購買履歴から、「商品Aを買った人は商品Bを買う確率が高い」といった関係性を見つけ出します。自分は、業務では予測分析がメインなので、そういえば久しくこういった手法に触れてなかったな…と感じました。というわけで、復習の一環で簡単にデモしてみようと思います*2

アソシエーション分析の説明は、例えば以下のサイトが参考になるかと思います。
sinhrks.hatenablog.com
sinhrks.hatenablog.com
実務では、今回使用するサンプルデータのような形式に簡単にまとめられない場合があるので、前処理も含めた解説は非常にためになります。

本当はPythonでやろうと思ったけど、自分の環境ではうまく構築できなかった

Pythonでは、Orangeというものでアソシエーション分析が可能ということで、pipでOrange(何も指定しないとVer. 2.7.8)をインストールしようとしました。しかし、setup.pyのビルドでコケてしまい、結局解消せず*3。大本のページ(Orange Data Mining)を見ると、最新版はVer. 3.3.7とのこと。pip経由だと、Orange3と指定するとインストールできそうなので実行。こちらはインストールに成功しましたが、アソシエーション分析に必要なメソッドが無い。OrangeでPythonからAprioriを動かす - Qiitaを見ると、どうやら本当に無い模様。

というわけで、いったんPythonで進めるのは保留にして、Rでデモをやるだけにします。

Rの{arules} + {arulesVIz}を使って分析

Rでアソシエーション分析を行う場合、例えば{arules}というパッケージがあります。今回は、このパッケージを使用していきます。

{arules}のインストール

if(!require(arules)){
  install.packages("arules", dependencies = TRUE)
}
require(arules)

データセットの読み込み

店舗の30日間のPOSデータを使用します。9835トランザクション、196カテゴリからなるデータで、「milk(牛乳)」、「root vegitables(根菜類)」といった商品の併売情報が入っています。

data(Groceries)
Groceries

# 出力
transactions in sparse format with
 9835 transactions (rows) and
 169 items (columns)

先頭の数行を確認してみます。併売情報を出力するのに、inspect()を用います。

inspect(head(Groceries))

# 出力
  items                     
1 {citrus fruit,            
   semi-finished bread,     
   margarine,               
   ready soups}             
2 {tropical fruit,          
   yogurt,                  
   coffee}                  
3 {whole milk}              
4 {pip fruit,               
   yogurt,                  
   cream cheese ,           
   meat spreads}            
5 {other vegetables,        
   whole milk,              
   condensed milk,          
   long life bakery product}
6 {whole milk,              
   butter,                  
   yogurt,                  
   rice,                    
   abrasive cleaner}

このようなデータを作ったり、アクセスしたりするのは一筋縄ではいかないのですが、ここではうまくトランザクションデータが手に入っているという状況で先に進みます。

summaryで基本統計量の確認

summary(Groceries)

# 出力
transactions as itemMatrix in sparse format with
 9835 rows (elements/itemsets/transactions) and
 169 columns (items) and a density of 0.02609146 

most frequent items:
      whole milk other vegetables       rolls/buns             soda           yogurt 
            2513             1903             1809             1715             1372 
         (Other) 
           34055 

element (itemset/transaction) length distribution:
sizes
   1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17 
2159 1643 1299 1005  855  645  545  438  350  246  182  117   78   77   55   46   29 
  18   19   20   21   22   23   24   26   27   28   29   32 
  14   14    9   11    4    6    1    1    1    1    3    1 

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.000   2.000   3.000   4.409   6.000  32.000 

includes extended item information - examples:
       labels  level2           level1
1 frankfurter sausage meat and sausage
2     sausage sausage meat and sausage
3  liver loaf sausage meat and sausage

出力結果の意味は以下です。

  • most frequent items : 頻度(=商品が売れた数)
  • element (itemset/transaction) length distribution : 1トランザクションで購入した商品の数
  • Min, 1st Qu., Median, Mean, 3rd Qu., Max : 最小値、第一四分位数、中央値、平均値、第三四分位数、最大値

商品ごとのヒストグラム

先程の要約だと、商品ごとの頻度の情報が十分ではありませんでした。ここでは、購入された商品に対して相対出現頻度によるヒストグラムを出力します。これは、itemFrequencyPlot()で実行できます。全商品を可視化すると複雑になるので、引数のsupportで支持度4%以上の商品のみに限定しています。

# 目盛りを無駄にメイリオに設定
quartzFonts(Hir = c("HiraMaruPro-W4", "HiraMaruProN-W4", 
                    "Meiryo", "HiraMinProN-W3"))
itemFrequencyPlot(Groceries, support = 0.04, cex.names = 0.7, 
                  col = "red", family = "Hir", font = 3)

f:id:tekenuko:20160924005258p:plain

アソシエーション分析の実行

apriori()で実行できます。support(支持度)とconfidence(確信度)は、それぞれ手で0.5%以上、1%以上をしきい値としました。

rules <- apriori(Groceries, parameter = list(support = 0.005, confidence = 0.01))

# 表示結果
Apriori

Parameter specification:
 confidence minval smax arem  aval originalSupport support minlen maxlen target   ext
       0.01    0.1    1 none FALSE            TRUE   0.005      1     10  rules FALSE

Algorithmic control:
 filter tree heap memopt load sort verbose
    0.1 TRUE TRUE  FALSE TRUE    2    TRUE

Absolute minimum support count: 49 

set item appearances ...[0 item(s)] done [0.00s].
set transactions ...[169 item(s), 9835 transaction(s)] done [0.00s].
sorting and recoding items ... [120 item(s)] done [0.00s].
creating transaction tree ... done [0.00s].
checking subsets of size 1 2 3 4 done [0.00s].
writing ... [2138 rule(s)] done [0.00s].
creating S4 object  ... done [0.00s].

この条件だと併売に関する2138個のルールが出来上がります。

特定の商品に注目した場合

アソシエーションのルールは、「X ならば Y」という形になっています。作成したルールを用いて、Yを固定したときにXがどのようなものがあるかを確認してみます。

# Y = "onions"(たまねぎ)にしてみる
onionRules <- subset(rules, subset = rhs %in% "onions")
# リフト値が高い順に表示
inspect(sort(onionRules, by = "lift"))

# 出力
     lhs                                   rhs      support     confidence lift     
1304 {root vegetables,other vegetables} => {onions} 0.005693950 0.12017167 3.8750440
1307 {other vegetables,whole milk}      => {onions} 0.006609049 0.08831522 2.8478038
316  {root vegetables}                  => {onions} 0.009456024 0.08675373 2.7974523
324  {other vegetables}                 => {onions} 0.014234875 0.07356805 2.3722681
308  {whipped/sour cream}               => {onions} 0.005083884 0.07092199 2.2869434
310  {citrus fruit}                     => {onions} 0.005592272 0.06756757 2.1787771
314  {tropical fruit}                   => {onions} 0.005693950 0.05426357 1.7497776
312  {bottled water}                    => {onions} 0.005897306 0.05335787 1.7205725
320  {yogurt}                           => {onions} 0.007219115 0.05174927 1.6687019
326  {whole milk}                       => {onions} 0.012099644 0.04735376 1.5269647
322  {rolls/buns}                       => {onions} 0.006812405 0.03703704 1.1942927
47   {}                                 => {onions} 0.031011693 0.03101169 1.0000000
318  {soda}                             => {onions} 0.005287239 0.03032070 0.9777183

これより「root vegetables(根菜類), other vegetables(他の野菜)」やroot vegetables(根菜類), whole vegetables(牛乳)」といった条件がonions(たまねぎ)に結びつくものがリフト値が高い(単に玉ねぎを買うよりよりもXを買ってさらにYを買う確率が高い)ことがわかります。

結果の可視化

{arulesViz}のインストール

if(!require(arulesViz)){
  install.packages("arulesViz", dependencies = TRUE)
}
require(arulesViz)

(図が単純になるように)ルールを再設定

rules <- apriori(Groceries, parameter = list(support = 0.004, 
                                             confidence = 0.5))

以下では、有用そうな図をいくつか紹介していきます。

散布図

散布図は、method = "scatterplot"で指定します。

graphics::plot(rules, method = "scatterplot")

f:id:tekenuko:20160924011554p:plain

横軸が支持度、縦軸が確信度で、点の色が赤に近いほどリフト値が高くなっています。これより、例えばリフト値が高い領域がわかり、パラメータ設定の参考にすることができます。

バブルチャート

method = "grouped"と指定することで、XとYの関係を直観的(丸の大きさと色)で表現できます。

graphics::plot(rules, method = "grouped")

f:id:tekenuko:20160924011607p:plain

横軸がアソシエーションの左辺(X)、縦軸が右辺(Y)、支持度の大きさは丸のサイズ、リフト値の大きさは色の濃さで表されています。

アソシエーショングラフ

有向グラフ形式(矢印で関係が表現されているグラフ)で表現することも可能です。表示を簡単にするため、リフト値の上位10項目のみに着目します。method = "graph"で有向グラフを表現できます。

# リフト値の高い順に10項目
rules_high_lift <- head(sort(rules, by = "lift"), 10)
graphics::plot(rules_high_lift, method = "graph", control = list(type = "items"))

f:id:tekenuko:20160924012017p:plain

herbsを購入した人は、root vegetablesを購入する可能性が高い(?)などが表示されます。
実はWindowとMacでこの表示が異なるそうで、Windowsはマウス操作で移動ができるそうです。

おわりに

アソシエーション分析を用いると、複雑なトランザクションデータだけを見ている場合にはわからなかったインサイトが得られる可能性があります。しかしながら、ルールがわかってもそれがビジネス戦略とあっているか*4はまた別の問題です*5。アソシエーション分析で課題解決を図る際は、きちんと目的を明確にし、分析が目的に沿っているかを精査しておく必要があります。自分が業務で携わるときは、このあたりを肝に命じて進めていこうと思います。

*1:一応自分はコンサルタントなので、たびたび全体感を整理しとくのは業務でも有用かなと思っています。

*2:というわけで非常にクオリティが低い内容です。。

*3:自分はPython3.5を使っているのが原因の一つっぽいですが、setup.pyをPython3に合うように書き換えても成功には至らず。。

*4:売上を上げたいなど。

*5:ルールが利益に結びつかない可能性があります。

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