XGBoost

R と Python で XGBoost (eXtreme Gradient Boosting) を試してみたのでメモ。

Boosting

バギング (Bootstrap aggregating; bagging) が弱学習器を独立的に学習していくのに対して, ブースティング (Boosting) は弱学習器の学習を逐次的に行います。
各ステップごとに前の弱学習器の結果を参考にしながら新たな弱学習器の重みを変化させ間違ったものをうまく識別できるように損失関数を最小化します。学習速度はバギングに比べ直列性のため遅くなる傾向があります。

バギングの代表例はランダムフォレスト (Random Forest), ブースティングの代表例は AdaBoost [1] で, 重みを適応的 (adaptive) に更新していく手法 [2]です。OpenCVのオブジェクト検出で使われています。

オブジェクト検出では画像の局所的な部分における明暗関係などの特徴量について, 弱識別器を作り評価値が高い (エラー率の低い) 識別器とその時の重みを採用し, サンプルを更新していきます。最後に弱識別器の出力に α で重み付けして, 結合することで最終的な識別器とします。この学習器を使って, 検出ウィンドウのラスタスキャンで入力をずらしていき, 識別器を通した出力が閾値を超えた時に検出とします。

Gradient Boosting

Gradient Boosting は Boosting の一種で, 損失関数の勾配を用います。詳細は Gradient boosting machines, a tutorial [Natekin and Knoll, 2013] が参考になります。

XGBoost (eXtreme Gradient Boosting) [3] は勾配ブースティングのC++実装で, 大規模データに対してもスケールできるように高速化の工夫がされています。

XGBoostインストール

環境はOSX 10.10.2, R 3.1です。
現時点では CRAN に {xgboost} パッケージが見つからなかったため source からインストールします。

$ git clone https://github.com/dmlc/xgboost.git
$ R
> install.packages("devtools")
> require(devtools)
> install_local('xgboost/', subdir = 'R-package')

RでXGBoostを試してみる

定番のirisデータセットで XGBoost を試してみます。

require(xgboost)

# k-fold
data(iris)
num <- sample(nrow(iris), 150)
k <- 2
data <- split(iris[num,], 1:k)

model <- xgboost(
    params=list(objective="multi:softmax"),
    eval_metric="mlogloss",
    num_class=3,
    data=as.matrix(data$"1"[,1:4]),
    label=as.integer(data$"1"[,5])-1,
    nrounds=500
  )

pred <- predict(model, as.matrix(data$"2"[,1:4]))
table(as.integer(data$"2"[,5])-1, pred)
# pred
#   0  1  2
# 0 28  0  0
# 1  0 24  3
# 2  0  2 18

sum(ifelse(as.integer(data$"2"[,5])-1 == pred, 0, 1)) / (nrow(iris)/k)
# [1] 0.06666667

imp <- xgb.importance(names(iris), model=model)
print(imp)
# Feature        Gain      Cover  Frequence
# 1: Petal.Length 0.882104244 0.52692315 0.44000000
# 2:  Petal.Width 0.103582149 0.25266661 0.24153846
# 3:  Sepal.Width 0.008033523 0.06575587 0.09384615
# 4: Sepal.Length 0.006280085 0.15465438 0.22461538

Code はGitHubにあります。

PythonでXGBoostを試してみる

PythonでXGBoostを使うには XGBoost Python Package を参考にインストールします。手っ取り早く使いたかったので pip のヴァージョン指定でインストールしました。

$ pip install xgboost==0.4a30
# Updateの場合は以下
# $ pip install -U pip
# $ sudo pip install -U xgboost
# $ pip list | grep xgb
# xgboost (0.6a2)

定番の iris から。

# -*- coding: utf-8 -*-

import numpy as np
import scipy as sp
import xgboost as xgb
from sklearn import datasets
from sklearn.metrics import confusion_matrix
from sklearn.grid_search import GridSearchCV
from sklearn.grid_search import RandomizedSearchCV

iris = datasets.load_iris()
train_x = iris.data[0::2, :]
train_y = iris.target[0::2]
test_x = iris.data[1::2, :]
test_y = iris.target[1::2]

np.random.seed(131)

# Grid Search
params = {'learning_rate': [0.05, 0.1],
          'max_depth': [3, 5],
          'subsample': [0.9, 0.95],
          'colsample_bytree': [0.5, 1.0]
          }

clf = xgb.XGBClassifier()

cv = GridSearchCV(clf,
                  params,
                  cv=10,
                  scoring="log_loss",
                  n_jobs=1,
                  verbose=2)
cv.fit(train_x, train_y)
predict = cv.predict(test_x)

print confusion_matrix(test_y, predict)
print cv.best_estimator_

best_param = {'learning_rate': cv.best_estimator_.learning_rate,
              'max_depth': cv.best_estimator_.max_depth,
              'gamma': cv.best_estimator_.gamma,
              'silent': cv.best_estimator_.silent,
              'objective': cv.best_estimator_.objective,
              'num_class': 3
              }

dtrain = xgb.DMatrix(np.array(train_x), label=train_y)
num_round = 10
bst = xgb.train(best_param, dtrain, num_round)
bst.save_model('./xgb/iris01.model')

irisは分類問題だったので classifier を使いましたが, 回帰問題では regressor を使うことができます。

Kaggle の Let’s try to run gs on xgb regressor (Apache 2.0 open source license.)を参考に試してみます。

上記では, スコアリングにジニ係数を用いています。

import numpy as np
import scipy as sp
import pandas as pd
import xgboost as xgb
from sklearn import pipeline, metrics, grid_search

def gini(y_true, y_pred):
    """ Simple implementation of the (normalized) gini score in numpy.
        Fully vectorized, no python loops, zips, etc. Significantly
        (>30x) faster than previous implementions

        Credit: https://www.kaggle.com/jpopham91/
    """

    # check and get number of samples
    assert y_true.shape == y_pred.shape
    n_samples = y_true.shape[0]

    # sort rows on prediction column
    # (from largest to smallest)
    arr = np.array([y_true, y_pred]).transpose()
    true_order = arr[arr[:,0].argsort()][::-1,0]
    pred_order = arr[arr[:,1].argsort()][::-1,0]

    # get Lorenz curves
    L_true = np.cumsum(true_order) / np.sum(true_order)
    L_pred = np.cumsum(pred_order) / np.sum(pred_order)
    L_ones = np.linspace(0, 1, n_samples)

    # get Gini coefficients (area between curves)
    G_true = np.sum(L_ones - L_true)
    G_pred = np.sum(L_ones - L_pred)

    # normalize to true Gini coefficient
    return G_pred/G_true

def normalized_gini(y_true, y_pred):
    ng = gini(y_true, y_pred)/gini(y_true, y_true)
    return ng

def fit(train, target):

    # set up pipeline
    est = pipeline.Pipeline([
            ('xgb', xgb.XGBRegressor(silent=True)),
        ])

    # create param grid for grid search
    params = {
        'xgb__learning_rate': [0.05, 0.1, 0.3, ],
        'xgb__min_child_weight': [5, 6, 7, ],
        'xgb__subsample': [0.5, 0.7, 0.9, 0.95,],
        'xgb__colsample_bytree': [0.5, 0.7, 0.9, ],
        'xgb__max_depth': [1, 3, 5, 7, 9, 11, ],
        'xgb__n_estimators': [10, 50, 100, ],
    }

    # set up scoring mechanism
    gini_scorer = metrics.make_scorer(normalized_gini, greater_is_better=True)

    # initialize gridsearch
    gridsearch = grid_search.RandomizedSearchCV(
        estimator=est,
        param_distributions=params,
        scoring=gini_scorer,
        verbose=10,
        n_jobs=-1,
        cv=10,
        n_iter=5,
    )

    # fit gridsearch
    gridsearch.fit(train, target)
    print('Best score: %.3f' % gridsearch.best_score_)
    print('Best params:')
    for k, v in sorted(gridsearch.best_params_.items()):
        print("\t%s: %r" % (k, v))

    # get best estimator
    return gridsearch.best_estimator_

def predict(est, testX):
    return est.predict(testX)

def save(predY, testY):
    diff = predY - testY
    pred = pd.DataFrame({'predict': predY, 'actual': testY, 'diff': diff})
    pred.to_csv('../data/xgb.csv')

    print('MSE: %.3f' % np.mean(np.absolute(diff)))
    print('Predictions saved to xgb.csv')

if __name__ == '__main__':
    np.random.seed(131)

    # load data
    df = pd.read_csv("data.csv", header=0)

    # create dataset
    trainX = data[0:800,[4,5,6,7,9,10,11,12]]
    trainY = data[0:800,8]
    testX = data[801:1198,[4,5,6,7,9,10,11,12]]
    testY = data[801:1198,8]

    # fit model
    est = fit(trainX, trainY)

    predY = predict(est, testX)
    save(predY, testY)

[1] A Decision-Theoretic Generalization of on-Line Learning and an Application to Boosting [Freund and Schapire, 1995]
[2] 出力が離散値の場合は, Discrete AdaBoost, 実数値の場合は Real AdaBoost。実数の場合は重み α は実質的に要らなくなる。
[3] XGBoost: A Scalable Tree Boosting System [Chen and Guestrin, 2016]
[4] 勾配ブースティングについてざっくりと説明する

タグ: XGBoost

【Boosting】R と Python で XGBoost