stacked generalization

[概要]
最近のkaggle のコンペのwinning solution で、stacked generalization がよく使われています。これの元になった論文は、1992 年のWolpert さんによるものです。
triskelion さんのブログKaggle Ensembling Guide | MLWave
の中でもこの手法についての説明があります。

様々な学習器を上手く組み合わせて、より精度の良いモデルを作ろうというのが基本的な考え方です。具体的には次の図のような感じです。

level 0 は、元となるデータです。またこの場合における各学習器はgeneralizer と呼ばれています。level 0 のデータにgeneralizer を適用して生成されたデータがlevel 1 のデータとなります。
その後も、同様に名づけられています。

[過去のコンペ]
まずは、多層パーセプトロンのように単純に組み合わせた場合の例を挙げます。
Otto のコンペでは、1 位、2 位、4 位の方はこの手法を使われていました。これは多クラス分類で、9 つのクラスに分ける問題で、評価関数はmlogloss です。
4 位の方はgithub にコードをupload されていますので、気になる方は参照してください。
1 位の方のモデル(kaggle のforum にて):
1st PLACE - WINNER SOLUTION - Gilberto Titericz & Stanislav Semenov - Otto Group Product Classification Challenge | Kaggle
2 位の方のモデル(kaggle blog にて):
Otto Product Classification Winner’s Interview: 2nd place, Alexander Guschin ¯\_(ツ)_/¯ | No Free Hunch
4 位の方のコード(github):
GitHub - diefimov/otto_2015

応用版は、さらに複雑なネットワークを形成したモデルです。
Dato のコンペは2 値分類で、評価関数はauc です。
Dato Winners’ Interview: 1st place, Mad Professors | No Free Hunch

こちらでは、階層をまたいだグラフになっています。

[訓練データの作り方]
まず検証用のデータは、いつも通りに訓練データを用いてモデルを作成して検証用のデータに適用します。
次に、1 段階上で使用する訓練データを作成します。ここでは、nfold =2 の場合で概略を説明します。

cross validation の要領で、validation data の予測値を集めたものが次の段階で使用する訓練データになります。
nfold を増やすと、次の段階で使用するトレーニングデータの平均値が良くなります。しかしその反面、分散(標準偏差)が増えたり、計算に時間が掛かったりします。
分散が増えると、元のデータと異なる傾向のものになってしまう危険性があります。すなわち、cv score は改善しているのに、提出してみると酷いスコアが出てしまいます><
分散が大きいときの対策としては、nfold を小さくする、あるいはbagging をすることが考えられます。

自明なことですが、下位のモデルが改善すれば、それが伝播してより上位のgeneralizer が生成するモデルも良くなります。
またstacker model は、stacked generalization を使用するつもりがなくても、xgboost + NeuralNet のような単純なensemble をする場合に、どのくらいの割合で混ぜるのかを調べるのにも有用です。

チューニングのやり方は他のモデルと同様に、

• cv score

• 変数重要度(xgboost やrandomForest などでは)

を見ながら行っていきます。

[実際に試してみました]
ここでは、Otto のコンペのデータで実際にstacked generalization を使用してやり直した際に使用したモデルです。
上記の図を再掲しますが、これがそのときに使用したモデルです。

コードはほぼcross validation を自分で書くだけのようなもので、やり方だけ確認して自分で組んでみるのが良いかなと思います。ちなみに、2nd level のxgboost による変数重要度の上位30 項目は次のようになっています。

各々の学習器の何が効いているのかが、ちらっと垣間見えます。また上記の図では、モデルを簡略化して書いていまして、実際のところ、1st level の学習器の詳細は次のようになっています。

• randomForest (n_estimators=300) + calibration

• extraTrees (n_estimators=300) + calibration

• xgboost (max_depth = 12, eta = 0.01)

• kNearestNeighbor (k=2,4,8,16,32,64,128,256,512,1024,2048)

• NeuralNet ( hidden-dropout の順に、中間層が次のような3 パターンを試しています。

512-0.5-512-0.5,
512-0.5-512-0.5-512-0.5,
512-0.5-512-0.5-512-0.5-512-0.5 )

• svm(cost = 10)

今回の評価関数はmlogloss で、各クラスの境界付近を各学習器がどのように分けているのかが非常に重要で、mlogloss が0.40 を下回るにはそこを意識する必要があったように思えました。
ちなみに、今回試したモデルでのスコアは次のようになってました。

Public LB score: 0.39756, Private LB score: 0.39869 でした。

ここで、使用したコードはgithub
github.com
にup していますので、気になる方はご覧ください。

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 を使用した場合は試せていません。