xgboost package のR とpython の違い
python と xgboost で検索をかけられている方も多く見受けられるので、R とほぼ重複した内容になりますが、記事にまとめておきます。python のxgboost のインストール方法はgithub を参考にされると良いと思います。github.com
R とpython のxgboost を使う際に感じる違い
R の利点
- 視覚化(visualization) が強い
- 自動化が簡単
- early stopping が簡単に使える
python の利点
- ハイパーパラメータのチューニングに hyperopt package が使用できる
現状として、R のpackage を使う方がメリットが大きいと思います。
まず、R の方から見ていきます。python でも主要な機能は実装されていますが、変数重要度を求めたときの視覚化が未実装で(計画はあるみたいです)、変数との対応も分かりにくいです。
実際にiris data で変数重要度を求めると、R の場合には以前のブログ
Xgboost のR における具体例 (クラス分類) - puyokwの日記
に記載したように、どの変数がどのくらい貢献しているのかが分かりやすく、print でもplot でも表示されます。python では、
import numpy as np import xgboost as xgb import pandas as pd from sklearn import datasets iris=datasets.load_iris() # python は0 base, R は1 base # 偶数番目をトレーニングデータ trainX=iris.data[1::2,:] trainY=iris.target[1::2] # 奇数番目を検証用データ testX=iris.data[2::2,:] testY=iris.target[2::2] # DMatrix 型に変換 trainX=pd.DataFrame(trainX) trainX.columns=iris.feature_names dtrain = xgb.DMatrix(trainX.as_matrix(),label=trainY.tolist()) testX=pd.DataFrame(testX) testX.columns=iris.feature_names dtest=xgb.DMatrix(testX.as_matrix()) np.random.seed(131) # シードを固定 # パラメータの設定 params={'objective': 'multi:softprob', 'eval_metric': 'mlogloss', 'eta': 0.3, # default=0.3 'max_depth': 6, # 木の深さの最大値 [default=6] 'min_child_weight': 1, # [default=1] 'subsample': 1, # [default=1] 'colsample_bytree': 1, # [default=1] 'num_class': 3 } # トレーニング bst=xgb.train(params,dtrain,num_boost_round=18) # 変数重要度を求める imp=bst.get_fscore() print(imp) # 以下のものが出力される {'f0': 15, 'f1': 7, 'f2': 57, 'f3': 28}
これは、f(列番号)となっています。。
一応、算出はできるもののR の方が使い勝手がいいように感じます。
R のもう一つの利点の自動化は、python でも交差検定できますが、
cv=xgb.cv(params,dtrain,num_boost_round=40,nfold=10) # 出力 # (中略) #tree prunning end, 1 roots, 0 extra nodes, 0 pruned nodes ,max_depth=0 #tree prunning end, 1 roots, 4 extra nodes, 0 pruned nodes ,max_depth=2 #tree prunning end, 1 roots, 2 extra nodes, 0 pruned nodes ,max_depth=1 #tree prunning end, 1 roots, 0 extra nodes, 0 pruned nodes ,max_depth=0 #tree prunning end, 1 roots, 2 extra nodes, 0 pruned nodes ,max_depth=1 #tree prunning end, 1 roots, 2 extra nodes, 0 pruned nodes ,max_depth=1 #[39] cv-test-mlogloss:0.150059+0.225514 cv-train-mlogloss:0.038407+0.003179 print(cv) # 出力 #'[0] cv-test-mlogloss:0.761828+0.039446 cv-train-mlogloss:0.752364+0.004754', #(中略) #'[39] cv-test-mlogloss:0.150059+0.225514 cv-train-mlogloss:0.038407+0.003179' bst=xgb.train(params,dtrain,18) pred=bst.predict(dtest) pred=pd.DataFrame(pred) print confusion_matrix(testY, pred.idxmax(axis=1)) # [[25 0 0] # [ 0 23 2] # [ 0 0 25]]
cv に含まれているのはただの文字列なので、これを処理するのは面倒に感じます。とりあえず、予測までしておきました。
R では以前の記事では紹介していなかったのですが簡単に最小値を探せます。(最後の2,3 行以外は以前の記事とほぼ同じコードです)
(2015/09/06: コードをより簡単な書き方に変更)
dim(iris) # 行数:150, 列数:5 odd.n<-2*(1:75)-1 iris.train<-iris[odd.n,] # 奇数を訓練データ iris.test<-iris[-odd.n,] # 偶数を検証データ library(xgboost) y <- iris.train[,5] # 目的変数 y <- as.integer(y)-1 #xgboost で既定されいるクラスは 0 base x <- iris.train[,1:4] # 説明変数 set.seed(131) # 固定シードで試す param <- list("objective" = "multi:softprob", # 多クラスの分類で各クラスに所属する確率を求める "eval_metric" = "mlogloss", # 損失関数の設定 "eta" = 0.3, # default= 0.3 "max_depth" = 6, # 木の深さの最大値 [default=6] "min_child_weight" = 1, # [default=1] "subsample" = 1, # [default=1] "colsample_bytree" = 1, # [default=1] "num_class" = 3 # class がいくつ存在するのか ) # 最適な木の数を探す cv.nround <- 100 #search bst.cv <- xgb.cv(param=param, data = data.matrix(x), label = y, nfold = 10, nrounds=cv.nround) min(bst.cv$test.mlogloss.mean) # [1] 0.108445 which.min(bst.cv$test.mlogloss.mean) # [1] 27
のようにして、mlogloss であれば、最小値とそのときのnround の値は取得できます。この値を用いて、検証データに適用すればOK です。逆に、最小値が先ほどnrounds の設定した値とほぼ同じであれば、nrounds の数を大きくしてもう一度交差検定するプログラムを書けば自動化できます。パラメータのチューニングも同様の方針でできます。
early stopping は以下のように簡単に使えます。(パラメータにearly.stop.round を追加)
set.seed(131) # 再現性の確保のため再度 bst.cv <- xgb.cv(param=param, data = data.matrix(x), label = y, nfold = 10, nrounds=cv.nround, early.stop.round=10) # 出力は以下のようになります # 一部省略 # [23] train-mlogloss:0.032116+0.002245 test-mlogloss:0.109858+0.128204 # [24] train-mlogloss:0.031495+0.002149 test-mlogloss:0.109487+0.128595 # [25] train-mlogloss:0.030918+0.002035 test-mlogloss:0.108704+0.127878 # [26] train-mlogloss:0.030378+0.001937 test-mlogloss:0.108445+0.127594 # [27] train-mlogloss:0.029885+0.001852 test-mlogloss:0.108808+0.127736 # [28] train-mlogloss:0.029433+0.001771 test-mlogloss:0.110306+0.129390 # [29] train-mlogloss:0.029008+0.001705 test-mlogloss:0.110110+0.129282 # [30] train-mlogloss:0.028601+0.001637 test-mlogloss:0.110021+0.129223 # [31] train-mlogloss:0.028222+0.001576 test-mlogloss:0.109943+0.129480 # [32] train-mlogloss:0.027860+0.001519 test-mlogloss:0.110860+0.130982 # [33] train-mlogloss:0.027510+0.001469 test-mlogloss:0.110752+0.131017 # [34] train-mlogloss:0.027191+0.001417 test-mlogloss:0.111216+0.131821 # [35] train-mlogloss:0.026877+0.001363 test-mlogloss:0.112265+0.133605 # [36] train-mlogloss:0.026584+0.001322 test-mlogloss:0.112999+0.135183 # Stopping. Best iteration: 27 # 予測する nround <- which.min(bst.cv$test.mlogloss.mean) bst <- xgboost(param=param,data=data.matrix(x),label=y,nrounds=nround) pred <- predict(bst,data.matrix(iris.test[,1:4])) pred <- matrix(pred,3,length(pred)/3) # 今回は3クラスあるので pred <- t(pred) table(iris.test[,5],max.col(pred)+1) # 2 3 4 # setosa 25 0 0 # versicolor 0 23 2 # virginica 0 0 25
直近 early.stop.round 個の間test-mlogloss がすべて上昇していればそこで打ち切ります。
[2015/08/09 追記あり] python でearly stopping を用いるのは、やや面倒ですが、hyperopt と組み合わせて用いた方法を紹介します。 kaggle の過去のコンペでhyperopt の使い方が紹介されていました。
Kaggle-stuff/hyperopt_xgboost.py at master · bamine/Kaggle-stuff · GitHub
こちらはそれを参考にしたコードです。以前にup したコードでは汎化性能などに問題があったため、改善しました。
import numpy as np import xgboost as xgb import pandas as pd from sklearn import datasets, cross_validation from sklearn.metrics import confusion_matrix, log_loss from sklearn.cross_validation import KFold from sklearn import preprocessing from hyperopt import hp from hyperopt import fmin, tpe, hp, STATUS_OK, Trials iris=datasets.load_iris() trainX=iris.data[0::2,:] trainY=iris.target[0::2] testX=iris.data[1::2,:] testY=iris.target[1::2] np.random.seed(131) def score(params): print "Training with params : " print params Sum=0.0 ite=0.0 for train, test in skf: X_train, X_test, y_train, y_test = trainX[train], trainX[test], trainY[train], trainY[test] dtrain = xgb.DMatrix(X_train, label=y_train) dvalid = xgb.DMatrix(X_test, label=y_test) watchlist = [(dtrain, 'train'),(dvalid, 'eval')] model = xgb.train(params, dtrain, num_boost_round=150,evals=watchlist,early_stopping_rounds=10) predictions = model.predict(dvalid).reshape((X_test.shape[0], 3)) ite+=model.best_iteration Sum+=model.best_score Sum/=len(skf) ite/=len(skf) print "\tAverage of best iteration {0}\n".format(ite) print "\tScore {0}\n\n".format(Sum) return {'loss': Sum, 'status': STATUS_OK} # hp.quniform('変数名',最小値,最大値,間隔) の形式 def optimize(trials): space = { 'eta' : hp.quniform('eta', 0.1, 0.5, 0.05), 'max_depth' : hp.quniform('max_depth', 1, 10, 1), 'min_child_weight' : hp.quniform('min_child_weight', 1, 10, 1), 'subsample' : hp.quniform('subsample', 0.5, 1, 0.05), 'gamma' : hp.quniform('gamma', 0, 1, 0.05), 'colsample_bytree' : hp.quniform('colsample_bytree', 0.5, 1, 0.05), 'num_class' : 3, 'eval_metric': 'mlogloss', 'objective': 'multi:softprob', 'lambda': 1e-5, 'alpha': 1e-5, 'silent' : 1 } best = fmin(score, space, algo=tpe.suggest, trials=trials, max_evals=250) print best skf = cross_validation.StratifiedKFold(trainY, n_folds=3, shuffle=True,random_state=123) trials = Trials() optimize(trials) # 最後に最適値が出力される #{'colsample_bytree': 1.0, 'min_child_weight': 1.0, 'subsample': 0.85, #'eta': 0.35, 'max_depth': 3.0, 'gamma': 0.15} params={'objective': 'multi:softprob', 'eval_metric': 'mlogloss', 'eta': 0.35, 'max_depth': 3, 'gamma': 0.15, 'min_child_weight': 1, 'subsample': 0.85, 'colsample_bytree': 1, 'lambda': 1e-5, 'alpha': 1e-5, 'num_class': 3 } score(params) # Average of best iteration 23.3333333333 # Score 0.0825176666667 trainX = pd.DataFrame(trainX) testX = pd.DataFrame(testX) dtrain = xgb.DMatrix(trainX.as_matrix(),label=trainY.tolist()) dtest=xgb.DMatrix(testX.as_matrix()) bst=xgb.train(params,dtrain,num_boost_round=23) pred=bst.predict(dtest) pred=pd.DataFrame(pred) print confusion_matrix(testY, pred.idxmax(axis=1)) #[[25 0 0] # [ 0 23 2] # [ 0 0 25]]
early stopping を用いる際にwatchlist が必要となりますが、early stopping の対象とするのはwatchlist の一番最後の項です。また、early stopping を使用すると、(上記のコードでは)model.best_iteration,model.best_score で最も良いスコアとそのときのiteration が分かります。
結果は、特に変わりませんでしたが、hyperopt を使ってパラメータのチューニングを行うことが、python でxgboost を使う魅力の一つになると思います。
Windows でLasagne (NeuralNet)
(情報が古いため近日中に新しく更新する予定です。最新版ではおそらく動かないと思います。追記:2016/06/02)
kaggle 勢でdeeplearning はよく用いられていますが、windows ユーザーは何かとインストールなどで苦労する点が多いように思えます。R ではh2O のdeeplearning ほぼ1択です。python では様々ありますが、個人的にはLasagne に落ち着いたのでこれのインストールと基本的な使い方のコードを紹介したいと思います。
まずLasagne の説明とコードは
Welcome to Lasagne — Lasagne 0.2.dev1 documentation
github.com
が参考になります。
windows ではWinPython を使うのが最も簡単に環境構築できました。インストールの注意点はPython 2.7の32bit 版(2015/07/31現在の最新版はWinPython-2.7.10.1)を使用することです(64bit 版も頑張ればできるかも)。
まず、環境変数の設定をします。
C:\<該当するパスを入れて下さい>\WinPython-32bit-2.7.10.1\python-2.7.10; C:\<該当するパスを入れて下さい>\WinPython-32bit-2.7.10.1\python-2.7.10\Scripts;
github からのダウンロードとLasagne のインストールです。
git clone https://github.com/Lasagne/Lasagne.git cd Lasagne pip install -r requirements.txt python setup.py install
WinPython-32-bit-2.7.10.1\python-2.7.10 の配下にあるpython.exe を直接実行してimport theano とするとエラーがでるのですが、なぜかWinPython-32-bit-2.7.10.1 にあるWinPython Command Prompt.exe を立ち上げて、こちらからpython.exe を実行し、import theano とすると無事にインポートできます。
WinPython Command Prompt.exe -> python.exe
の順に試して頂ければよいかと思います。
import theano import lasagne
エラーが出なければOKです。
サンプルコードとして、examples の中にmnist.py が含まれていますので、とりあえずこれを実行してみると雰囲気が分かるとは思います。個人的には、このコードは分かりにくいのであまり参考にするものではないかなという印象です。kaggle でよく使われているのはnolearn の方なので、こちらをインストールします。
git clone https://github.com/dnouri/nolearn.git cd nolearn pip install -r requirements.txt python setup.py develop
nolearn の詳細はgithub に記載されています。
github.com
基本的なコードの書き方で参考になりそうなのは、
Neural Net では
kaggle_otto/nn_adagrad.py at master · ahara/kaggle_otto · GitHub
convolutional Neural Net では
kfkd-tutorial/kfkd.py at master · dnouri/kfkd-tutorial · GitHub
があります。
最近の更新で、やや不安定になった感じ(random.seedが効かない、あるいはconvolutional Neural Net が動作しないなど)はありますが、
kaggle のdigit recognizer で使用しているコードの一部を抜粋して、それぞれのコードを紹介します。ここで使用するコードはgithub
GitHub - puyokw/kaggle_digitRecognizer: Scripts for digit Recognizer in kaggle
上にあります。
# multi layer perceptron layers0 = [('input', InputLayer), ('dense0', DenseLayer), ('dropout', DropoutLayer), ('dense1', DenseLayer), ('dropout', DropoutLayer), ('dense2', DenseLayer), ('dropout', DropoutLayer), ('output', DenseLayer)] net0 = NeuralNet(layers=layers0, input_shape=(None, num_features), dense0_num_units=1024, dropout_p=0.5, dense1_num_units=1024, dense2_num_units=2048, output_num_units=num_classes, output_nonlinearity=softmax, update=adagrad, eval_size=0.01, verbose=1, max_epochs=20, update_learning_rate=theano.shared(float32(0.01)), on_epoch_finished=[ AdjustVariable('update_learning_rate', start=0.01, stop=0.001),]) net0.fit(X, y) y_prob = net0.predict_proba(X_test)
# convolutional Neural Net net = NeuralNet( layers=[ ('input', InputLayer), ('conv1', Conv2DLayer), ('pool1', MaxPool2DLayer), ('dropout1', DropoutLayer), ('conv2', Conv2DLayer), ('pool2', MaxPool2DLayer), ('dropout2', DropoutLayer), ('hidden3', DenseLayer), ('dropout3', DropoutLayer), ('hidden4', DenseLayer), ('dropout4', DropoutLayer), ('output', DenseLayer), ], input_shape=(None, 1, 28, 28), conv1_num_filters=32, conv1_filter_size=(3, 3), pool1_pool_size=(2, 2), dropout1_p=0.2, conv2_num_filters=64, conv2_filter_size=(3, 3), pool2_pool_size=(2, 2), dropout2_p=0.2, hidden3_num_units=1024, dropout3_p=0.5, hidden4_num_units=1024, dropout4_p=0.5, output_num_units=num_classes, output_nonlinearity=softmax, update=adagrad, update_learning_rate=theano.shared(float32(0.03)), eval_size=0.01, on_epoch_finished=[ AdjustVariable('update_learning_rate', start=0.01, stop=0.0001), ], max_epochs=30, verbose=1, ) net.fit(X, y) y_prob = net.predict_proba(X_test)
これらが各々のニューラルネットの構成を指定している部分です。自分のノートPC ではconvolutional Neural Netの 1epoch に20分くらいかかっていました。gpu を使うと10倍以上速くなるらしいですが、金銭的な余裕がないのでgpu を使用した場合は試せていません。
次元削除 ( t-SNE )
今回は、kaggle のOtto Group Production Classification Challenge の上位の方々が次元削除の手法としてt-SNE(t-distributed stochastic neighbor embedding) を使用されていたので調べてみようと思いました。個人的には、pca(主成分分析) ぐらいしか思い付かなかったのですが、それぞれ比較しながら見ていきます。
t-sne の詳細についてこちらを参考にするといいかと思います。
http://jmlr.org/papers/volume9/vandermaaten08a/vandermaaten08a.pdf
こちらに書かれているようにt-SNE は高次元のものを2 または3 次元に写像するように作られています。とりあえず、R のtsne package を試してみます。(あとでより高速なRtsne package を用います)
library(tsne) colors = rainbow(length(unique(iris$Species))) names(colors) = unique(iris$Species) set.seed(1) iris.tsne = tsne(iris[,1:4], max_iter=200) # 時間短縮のため max_iter=200 plot(iris.tsne, t='n', main="tsne") text(iris.tsne, labels=as.integer(as.factor(iris$Species)), col=colors[iris$Species]) # PCA と比較 dev.new() iris.pca = princomp(iris[,1:4])$scores[,1:2] plot(iris.pca, t='n', main="pca") text(iris.pca, labels=as.integer(as.factor(iris$Species)),col=colors[iris$Species])
これを実行すると次のような図を得ます。
まず、t-sne
次に、pca
これを見た感じだとpca の方が良さそうです。
次に、先ほどの例ではt-SNE の良さがあまり分からなかったので、別の例として手書き文字(数字)の認識を例としてみていきたいと思います。t-SNE とpca を比較したスクリプトがkaggle の"script of the week" で取り上げていただきました！ので、そちらにリンクを張っておきます。
https://www.kaggle.com/puyokw/digit-recognizer/clustering-in-2-dimension-using-tsnewww.kaggle.com
結果はこのようになっています。
t-sne
pca
こちらでは、t-SNE がかなりきれいにクラス分類されているのが分かります。ただこの場合でも、pca はダメという分けではなく、単にどのような処理をしたいかに合わせてpca とt-SNE を選択する必要があるということです。
R にはt-SNE の高速版のRtsne package があり、これは主成分分析をしてからt-SNE を実行します(prcomp を使用しています)。
Rtsne についての論文はこちらになります。
http://arxiv.org/pdf/1301.3342v2.pdf
上記のページが意外と重かったので、アブストが書かれているページにもリンクしておきます。
[1301.3342] Barnes-Hut-SNE
ここに書かれているように、この方法は Barnes-Hut-SNE と呼ばれていてこのアルゴリズムではO(nlogn) で実行できます。元のアルゴリズムはO(n^2) です。これはvantage-point trees やBarnes-Hut のアルゴリズムを使用しています。 使い方はtsne とほぼ同じ感じです。
library(Rtsne) set.seed(1) # 再現性の確保 # verbose=TRUE で途中経過が出力されます iris.tsne <- Rtsne(as.matrix(iris[,1:4]), check_duplicates = FALSE, verbose=TRUE) colors = rainbow(length(unique(iris$Species))) plot(iris.tsne$Y, t='n', main="Rtsne") text(iris.tsne$Y, labels=as.integer(as.factor(iris$Species)), col=colors[iris$Species])
Rtsne のパラメータ中にtheta(=0.5 :default) が含まれています(上記のコードには入れていません)が、これを0 に近づけるほど精度は良くなりますがその代わりに実行時間が長くなります。
追記(2016/06/04)
一応、python での例も載せておきます。
import numpy as np import pandas as pd from sklearn.manifold import TSNE from sklearn import datasets import matplotlib.pyplot as plt # import iris data iris = datasets.load_iris() X = iris.data Y = iris.target model = TSNE(n_components=2, perplexity=50, n_iter=500, verbose=3, random_state=1) X = model.fit_transform(X) # Plot the points plt.figure(2, figsize=(8, 6)) plt.clf() plt.scatter(X[:, 0], X[:, 1], c=Y, cmap=plt.cm.Paired) plt.xlabel('tsne1') plt.ylabel('tsne2') plt.show()
追記(2016/05/23)
t-sne の利点:
他の視覚化よりも綺麗に図示されることが多い。(特にクラス分類)
t-sne の欠点:
計算に非常に時間がかかる
メモリの消費量が多い(特に、python)
kaggle における実例：
otto: Otto Group Product Classification Challenge | Kaggle
digit recognizer(上記のもの): Digit Recognizer | Kaggle
spring leaf marketing response(コードのみ): Solution Sharing - Springleaf Marketing Response | Kaggle
こちらのサイトに有名なデータセットに対するいくつかの例がExmaples の項目に載っています。t-SNE – Laurens van der Maaten
追記(2015/08/02)：実際に、t-SNE で3次元にしたものと元のデータを結合してrandomForest などを用いたコード(kaggle のdigit recognizerで用いたもの)をgithub 上にupload しました。精度は若干向上します。
GitHub - puyokw/kaggle_digitRecognizer: Scripts for digit Recognizer in kaggle
