Xgboost のR における具体例 (クラス分類)
前回、Xgboost のパラメータについて列挙しましたが、あれだけ見ても実際にどう使うのかよく分かりません。そこで今回はR で、とりあえず iris data を用いてその使い方を見ていきたいと思います。
まず、iris data の奇数番目を訓練データ、偶数番目を検証データとします。
dim(iris) # 行数:150, 列数:5 odd.n<-2*(1:75)-1 iris.train<-iris[odd.n,] # 奇数を訓練データ iris.test<-iris[-odd.n,] # 偶数を検証データ
次に、xgboost に適した形式に変換します。
library(xgboost) y <- iris.train[,5] # 目的変数 y <- as.integer(y)-1 #xgboost で既定されいるクラスは 0 base train.x<-iris.train[,1:4] x <- rbind(train.x,iris.test[,-5]) # xgboost を使うときのため x <- as.matrix(x) trind <- 1:length(y) # 先程定義したx の中の訓練データを指すのに使う teind <- (nrow(train.x)+1):nrow(x) # 先程定義したx の中の検証用データを指すのに使う
必要最低限のパラメータの設定は次のようになります。
set.seed(131) # 固定シードで試す param <- list("objective" = "multi:softprob", # 多クラスの分類で各クラスに所属する確率を求める "eval_metric" = "mlogloss", # 損失関数の設定 "num_class" = 3 # class がいくつ存在するのか )
実際に、予測をする前に最適な木の数を探します。
k<-round(1+log2(nrow(train.x))) cv.nround <- 100 #search bst.cv <- xgb.cv(param=param, data = x[trind,], label = y, nfold = k, nrounds=cv.nround)
これを実行すると、次のような結果が出てきます。
[0] train-mlogloss:0.752857+0.007342 test-mlogloss:0.772249+0.027355 [1] train-mlogloss:0.542712+0.008457 test-mlogloss:0.569818+0.047944 [2] train-mlogloss:0.404424+0.008669 test-mlogloss:0.441932+0.068164 [3] train-mlogloss:0.306743+0.007542 test-mlogloss:0.351062+0.084861 (略) [25] train-mlogloss:0.032532+0.003092 test-mlogloss:0.122064+0.128217 [26] train-mlogloss:0.031981+0.002944 test-mlogloss:0.121880+0.128820 [27] train-mlogloss:0.031499+0.002809 test-mlogloss:0.121960+0.128635 (略) [97] train-mlogloss:0.023760+0.001694 test-mlogloss:0.134890+0.156262 [98] train-mlogloss:0.023744+0.001690 test-mlogloss:0.135084+0.156458 [99] train-mlogloss:0.023728+0.001686 test-mlogloss:0.135296+0.156967
この場合、訓練データにおける最適な木の数は26+1(0 baseのため)となります。
これを基に予測をしてみます。
set.seed(131) nround <- 27 # モデルの構築 bst <- xgboost(param=param, data = x[trind,], label = y, nrounds=nround) pred <- predict(bst,x[teind,]) # モデルを使って予測値を算出 pred <- matrix(pred,3,length(pred)/3) # 今回は3クラスあるので pred <- t(pred) colnames(pred)<-c("setosa","versicolor","virginica")
結果を見てみると、
head(pred,3) setosa versicolor virginica [1,] 0.9838002 0.008915161 0.007284546 [2,] 0.9870363 0.008944485 0.004019236 [3,] 0.9870363 0.008944485 0.004019236
となっていて、各データが"setosa","versicolor","virginica"のどれに属するかの確率が求められています。この中の最大値をとるクラスに、そのデータは属すると思われます。
このような確率に興味がなく、どこのグループに属するかだけを知りたいなら次のように"multi:softmax"を使います。
param <- list("objective" = "multi:softmax", # multi:softmax に変更! "eval_metric" = "mlogloss", "num_class" = 3 ) set.seed(131) nround <- 27 bst <- xgboost(param=param, data = x[trind,], label = y, nrounds=nround) pred <- predict(bst,x[teind,])
データを整えて、結果を見てみます。
for(i in 1:length(pred)){ if(pred[i]==0) {pred[i]="setosa"} else if(pred[i]==1) {pred[i]="versicolor"} else {pred[i]="virginica"} } table(iris.test[,5],pred) # 結果 pred setosa versicolor virginica setosa 25 0 0 versicolor 0 23 2 virginica 0 0 25
これが良いのかどうかいまいち分からないので、とりあえずランダムフォレストと比較してみます。
odd.n<-2*(1:75)-1 iris.train<-iris[odd.n,] # 奇数を訓練データ iris.test<-iris[-odd.n,] # 偶数を検証データ # randomForest library(randomForest) set.seed(131) train.x<-iris.train[,1:4] train.y<-as.factor(iris.train[,5]) model.rf<-tuneRF(train.x,train.y,doBest=T) pred<-predict(model.rf,iris.test[,-5]) table(iris.test[,5],pred) # 結果 pred setosa versicolor virginica setosa 25 0 0 versicolor 0 24 1 virginica 0 2 23
これだけでは、なんとも言えなさそうですが、xgboost の性能の高さは認めらるかなと思います。
追記(2015/05/06): 変数重要度を求めることができることに気が付いたので追記しておきます。
# 変数重要度を求める imp<-xgb.importance(names(iris),model=bst) print(imp) # print の結果 Feature Gain Cover Frequence 1: Petal.Length 0.663081352 0.63384270 0.4453125 2: Petal.Width 0.318246635 0.24789925 0.2656250 3: Sepal.Length 0.011982305 0.09317686 0.2187500 4: Sepal.Width 0.006689708 0.02508119 0.0703125 # print の結果ここまで xgb.plot.importance(imp) # これをグラフで表示
変数重要度を図示したものがこちらになります。
ランダムフォレストで変数重要度がどうなるのかも気になるので試してみます。
# 変数重要度を求める(和が1になるように調整) print( model.rf$importance /sum(model.rf$importance) ) # print の結果 MeanDecreaseGini Sepal.Length 0.01306678 Sepal.Width 0.01388757 Petal.Length 0.42804382 Petal.Width 0.54500183 # print の結果ここまで varImpPlot(model.rf) # 図示(略)
xgboost とランダムフォレストでずいぶん変数重要度が異なっていました。異なる基準を採用しているのかまでは、まだ調べられていません。(分かり次第また更新しようと思います。)
追記(2015/6/22): 決定木ベースの手法なので、決定木も出力できた方が良いかと思ったので、追記しておきます。
# 決定木を表示 xgb.plot.tree(feature_names=names(iris[,-5]),model=bst, n_first_tree=2)
これを実行すると次のような図がhtml 形式のファイルで出力されます。
追記(2016/11/14): xgb.plot.tree を実行しても木が出力されないエラーに遭遇したのでコメントを付記しておきます。現在cran にあるパッケージではエラーがでるようでgithub からインストールすれば、正しく木が表示されました。
追記(2015/04/30):今回使用したコードはpublic gist にup しておきます。
introduction to xgboost · GitHub
今回用いたiris のデータでは、パラメータ調整してもほとんどその効果が見られないため、改めて別のデータで試してみるか検討してみます。(やる場合にはGW中にup する予定です) 2015/04/30 0:40 追記
xgboost のパラメータ
xgboost を使う上で、日本語のサイトが少ないと感じましたので、今回はパラメータについて、基本的にこちらのサイトの日本語訳です。
xgboost を実行する前に、共通変数、ブースター変数、タスク変数の3つをセットしなければなりません。
- 共通変数は、木あるいは線形モデルのブースティングモデルに共通の変数です
- ブースター変数は、木あるいは線形モデルの各々に固有の変数です
- タスク変数は、どのように学習させるかを決めるもので、例えば、回帰タスクではランキングタスクの変数が異なっています
- この3つの変数に加えて、コンソール変数があり、これはxgboost のコンソール版の動作に関連しています。(たとえば、作製したモデルを保存するときなど)
Rパッケージにおける変数
R のパッケージでは、. (ドット) を_ (アンダーバー)に置き換えて利用できます。例えば、max.depth はmax_depth とすると使えます。
共通変数
- booster [デフォルト=gbtree]
- gbtree もしくは gblinear のいずれのブースターを使うかを決めます。詳細はこちらにあります(英語: Boosters · dmlc/xgboost Wiki · GitHub ) 。
- silent [デフォルト=0]
- 0 にすると起動中のメッセージを出力し、1 にすると出力しません
- nthread [設定しなければ、デフォルト値として使用可能なスレッドの最大数になります]
- xgboost を使用時の並列処理を行うスレッドの数
- num_pbuffer [xgboost が自動的に設定するため、ユーザーが設定する必要はありません]
- 予測バッファのサイズで、たいていトレーニングデータ数で設定されます。ブースティングの最後の段階で予測の結果を保存するために使用されます。
- num_feature [xgboost が自動的に設定するため、ユーザーが設定する必要はありません]
- ブースティングに使用する特徴次元の数で、その特徴次元の最大数に設定されます
ブースター変数
xgboost-unity より、接頭辞 bst:
はブースター変数に付けなくても大丈夫です。変数に接頭辞 bst: を付けても付けなくても同じです。(たとえば、 bst:eta も eta も正しい変数で同じものを指しています)
木ブースター固有の変数
- eta [default=0.3]
ここで、強固と使っている意味は、外れ値に強くなるという意味と汎化性能が低くなるという意味の両方を込めています(良い点と悪い点)。以下も同様です。
- gamma
- 損失還元(loss reduction)の最小値で、木の葉ノードをさらに分割する際に必要な値。値を大きくすると、そのアルゴリズムはより強固になります。
- max_depth [default=6]
- 木の深さの最大値
- min_child_weight [default=1]
- 子ノードにおける必要な最小の重み。木の分割段階で、ある葉ノードにおける重みの合計値がmin_child_weight 未満であれば、それ以上分割しません。たとえば線形回帰の場合には、各ノードに必要な例の最小値に対応しています。値を大きくするとそのアルゴリズムはより強固になります。
- max_delta_step [default=0]
- 値を0に設定した場合、特になにもしません。もし正の値にした場合には、ステップの更新でより強固に(柔軟性がない) します。このパラメータは必要でないことが多いが、クラスのバランスが著しく悪いときのロジスティック回帰で役立つかもしれません。値を1~10に設定すると良いでしょう。
- subsample [default=1]
- サブサンプルを生成する際のトレーニングデータの抽出割合。たとえば、0.5に設定すると、XGBoost はデータの半分をランダムに選んで木を成長させることで、オーバーフィッティングを防ぎます。
- colsample_bytree [default=1]
- 各木を作成するときの列におけるサブサンプルの割合
線形ブースター固有の変数
タスク変数
- objective [ default=reg:linear ]
- 学習方法を明確にします。そのオプションは以下のようになります。
- "reg:linear" --線形回帰
- "reg:logistic" --ロジスティック回帰
- "binary:logistic" --2-クラス分類向けのロジスティック回帰で、確率を出力
- "binary:logitraw" --2-クラス分類向けのロジスティック回帰で、ロジスティック回帰式に代入する前の値
- "multi:softmax" --XGBoost をソフトマックス関数を使用した多クラス分類に設定します。変数num_class (クラスの数)も設定する必要があります。
- "multi:softprob" --softmax と同じであるが、各データがどのクラスに属するかの予測確率が出力されます。
- "rank:pairwise" --順位付けをペアワイズ損失が最小になるようにXGBoost を設定します
- base_score [ default=0.5 ]
- すべての例に対する初めの予測スコアでグローバルバイアスである
- eval_metric [ default according to objective ]
- 検証データに対する評価基準で、デフォルトの基準はobjective によって決まります。回帰であればrmse、クラス分類であればerror、順位付けではmean average precision がデフォルトとなります。
- 選択肢は以下の通りです
- "rmse": 2乗平均平方根誤差
- "logloss": 負の対数尤度
- "error": 2-クラス分類のエラー率。(分類の誤ったケース数) / (すべてのケース数) で計算される。 たとえば予測において、この評価値は0.5 よりも大きければ予測は良い例であり、0.5 よりも小さい場合には悪い例となる。
- "merror": 多クラス分類のエラー率。(誤認したケースの数) / (すべてのケース数).
- "mlogloss": 多クラスの対数損失
- "auc": Area under the curve の略でROC曲線下の面積で、性能の良さを表します
- "ndcg":Normalized Discounted Cumulative Gain の略で順位付け問題の性能評価に用いられる
- "map": Mean average precision の略
- seed [ default=0 ]
- 乱数シード値
コンソール変数
以下の変数はコンソール版のxgboost にのみ適用されます(一部省略)
- use_buffer [ default=1 ]
- テキストからの入力のためにバイナリバッファを作るかどうかで、たいていの場合データの読み込みのスピードが改善します
- num_round
- ブースティングを行う回数
- data
- トレーニングデータのパス
- test:data
- 予測する用のテストデータのパス
*1:強固になるの意味は、外れ値に強くなるという意味と汎化性能が低くなるの両方の意味を込めています。
トレーニングデータが小さい ><
kaggle の Restaurant Revenue Predicton に参加していますが、このコンテストのデータセットは、トレーニングデータが非常に小さく(n=137)、テストデータが大きい(n=10,000)のが特徴です。
以前にもこのような傾向のコンテストがあったらしいです。そのときに優勝された方の手法 で特徴的に思えたのが、「サポートベクトルマシンで cost の値を非常に大きくしている」ところです。実際にこのコンテストでも、使用されている方は見受けられます(Forum を見ている感じで)。
また、トレーニングデータが小さい場合に気を付けるポイントは、交差検定の部分にもあります。
交差検定の k の値はどれくらいにすればいいのか - ほくそ笑む
こちらにもありますように普通の場合には、k=サンプル数 のときが最も良い評価が出ると考えられます。ところが、トレーニングデータが小さい場合に、同様のことをすると過適合してしまいます。kaggle の Forum でもこのことが話題に上がっていました。
How to correct cross-validate? - Restaurant Revenue Prediction | Kaggle
そこでまず、k=10 (gbm のデフォルト値でもある) を試す、あるいは先ほどのサイトにもあったように k=lg( データ数 ) とするのが良さそうです。(lg は底が2の対数)
Forum で参考文献が挙げらていたのは、小規模のデータセットに対するSVMのパラメータのチューニングについては
http://robotics.stanford.edu/~ang/papers/cv-final.pdf
があるらしいです。また過適合の際のモデル選択については
http://jmlr.csail.mit.edu/papers/volume11/cawley10a/cawley10a.pdf
が参考に挙げられていました。
追記(2016/01/31): このデータセットのやり直しをしました!
今回のモデルは、randomForest の10 回平均(10 個のsingle model の平均) で行いました。score は
Public LB score: 1713786.67770, Private LB score: 1753006.78906
です。
randomForest の10 回平均についてですが、このデータセットではシード値による予測値のブレが非常に大きかったため、10 回平均を取ることにしました。こうすることで、実際にcv score が安定して手元でモデルの比較を行うことが以前よりはできました。
randomForest を行う前に、前処理として欠損値補完を行っております。
P14 ~ P18, P24 ~ P27, P30 ~ P37 が同時にゼロになるところがあり、そこを欠損値とみなしてmice package で補完しました。ここはwining solution を参考にしました。
ただし、wining solution では、Time / Age related information を抽出して、これを欠損値補完に役立てたらしいのですが、いまいち意味が分からなくてそこの実装ができませんでした。また、wininng solution では、gbm のsingle model を使用している点も異なっています。