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 を使う魅力の一つになると思います。
