Spark MllibでのClick through rate(CRT) 予測

1年以上も前ですが、edXのScalable Machine learningを受講していた時のコードが出てきたので、見直してみました。この講座のネタが、KaggleのCriteoのコンペのデータ
を基にしたCTR予測でして、One-hot-encodingやhushingした高次元のデータを
Mllibのロジスティック回帰モデル
)で予測をするといったものです。今回は予測精度そのものよりも、ハッシュ関数で次元圧縮したものが、OHE化したデータと比べてどの程度予測精度に差が出るのかを見ていきます。


Sparkもかなり前のバージョンだったので、今ならまた違う実装の仕方があるような気もしますが、出来合いの機能に頼らずに実装したことで理解が深まった記憶があるので、そのまんま。追加機能の勉強兼ねた比較はまた別のタイミングで行いたいなと思います。 SparkはVirtualbox上にたてたUbuntuにインストールしたものをシングルノードで使用しています。

おおまかな手順

①データの準備
元データを訓練、評価、テストデータにそれぞれ分割
②元データにOHEを適用。ここで特徴量は20万以上になります。
③OHE後のデータでロジスティック回帰で分析(model1)
④ハッシュトリックで次元量の削減(およそ3,000程度)
⑤ハッシュ後のデータでロジスティック回帰分析
-Mllibデフォルト設定のモデル(model2)
-グリッドサーチでパラメータ設定したモデル(model3)
⑥Loglosによるモデル評価
-OHEデータモデル(model1) vs Hushingデータモデル(model2)
-Hushingデータ・デフォルトモデル(model2) vs Hushingデータ・GSモデル(model3)

コード

まずはデータを読み込んで分割、眺めてみます。

import numpy as np
from pyspark.mllib.linalg import SparseVector
from pyspark.mllib.regression import LabeledPoint

rawData = (sc.textFile("dac_sample.txt").map(lambda x: x.replace('\t', ',')))

OneSample = rawData.take(1)
print OneSample


weights = [.8, .1, .1]
seed = 42
rawTrainData, rawValidationData, rawTestData = rawData.randomSplit(weights, seed)

#Cache each datasets as it is used repeatedly
rawTrainData.cache()
rawValidationData.cache()
rawTestData.cache()

nTrain = rawTrainData.count()
nVal = rawValidationData.count()
nTest = rawTestData.count()
print nTrain, nVal, nTest, nTrain + nVal + nTest


[u'0,1,1,5,0,1382,4,15,2,181,1,2,,2,68fd1e64,80e26c9b,fb936136,7b4723c4,25c83c98,7e0ccccf,de7995b8,1f89b562,a73ee510,a8cd5504,b2cb9c98,37c9c164,2824a5f6,1adce6ef,8ba8b39a,891b62e7,e5ba7672,f54016b9,21ddcdc9,b1252a9d,07b5194c,,3a171ecb,c5c50484,e8b83407,9727dd16']
79911 10075 10014 100000


OHEを分割したデータに当てます。

def createOneHotDict(inputData):
    OHEDict = (inputData
               .flatMap(lambda x: x)
               .distinct()
               .zipWithIndex()
               .collectAsMap())
    return OHEDict


def parsePoint(point):
    items = point.split(',')
    return [(i, item) for i, item in enumerate(items[1:])]
    
def oneHotEncoding(rawFeats, OHEDict, numOHEFeats):
    sizeList = [OHEDict[f] for f in rawFeats if f in OHEDict]
    sortedSizeList = sorted(sizeList)
    valueList = [1 for f in rawFeats if f in OHEDict ]
    return SparseVector(numOHEFeats, sortedSizeList, valueList)


def parseOHEPoint(point, OHEDict, numOHEFeats):
    parsedPoints = parsePoint(point)
    items = point.split(',')
    label = items[0]
    features = oneHotEncoding(parsedPoints, OHEDict, numOHEFeats)
    return LabeledPoint(label, features)

parsedFeat = rawTrainData.map(parsePoint)
ctrOHEDict = createOneHotDict(parsedFeat)
numCtrOHEFeats = len(ctrOHEDict.keys())


OHETrainData = rawTrainData.map(lambda point: parseOHEPoint(point, ctrOHEDict, numCtrOHEFeats)).cache()
OHEValidationData = rawValidationData.map(lambda point: parseOHEPoint(point, ctrOHEDict, numCtrOHEFeats)).cache()
OHETestData = rawTestData.map(lambda point: parseOHEPoint(point, ctrOHEDict, numCtrOHEFeats)).cache()

print ('Feature size after OHE:\n\tNumber of features = {0}'
       .format(numCtrOHEFeats))

print OHETrainData.take(1)

Feature size after OHE:
	Number of features = 233286
[LabeledPoint(0.0, (233286,[386,3077,6799,8264,8862,11800,12802,16125,17551,18566,29331,33132,39525,55794,61786,81396,82659,93573,96929,100677,109699,110646,112132,120260,128596,132397,132803,140620,160675,185498,190370,191146,195925,202664,204273,206055,222737,225958,229942],[1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0]))]


OHE後のデータにてロジスティック回帰を実施。とりあえずパラメーター設定はデフォルトのままやってみます。

from pyspark.mllib.classification import LogisticRegressionWithSGD

numIters = 100
stepSize = 1.0
regParam = 0.01
regType = 'l2'
includeIntercept = True

model0 = LogisticRegressionWithSGD.train(data = OHETrainData, 
                                         iterations = numIters, 
                                         step = stepSize, 
                                         regParam = regParam,
                                         regType = regType,
                                         intercept = includeIntercept)

#Compute loglos
from math import log, exp 

def computeLogLoss(p, y):
    epsilon = 10e-12
    if y == 1:
        return -log(epsilon + p) if p == 0 else -log(p)
    elif y == 0:
        return -log(1 - p + epsilon) if p == 1 else -log(1 - p)

def getP(x, w, intercept):
    rawPrediction = x.dot(w) + intercept

    # Bound the raw prediction value
    rawPrediction = min(rawPrediction, 20)
    rawPrediction = max(rawPrediction, -20)
    return 1.0 / (1.0 + exp(-rawPrediction))

def evaluateResults(model, data):
    return data.map(lambda x: computeLogLoss(getP(x.features, model.weights, model.intercept), x.label)).sum() / data.count()

logLossValLR0 = evaluateResults(model0, OHEValidationData)
print ('Validation Logloss for model1:\n\t = {0}'
       .format(logLossValLR0))


素うどんのロジスティック回帰でどこまで予測できてるか、ROC曲線を描いてみます。ぼちぼちの精度でできてるようでした。

labelsAndScores = OHEValidationData.map(lambda lp: (lp.label, getP(lp.features, model0.weights, model0.intercept)))
labelsAndWeights = labelsAndScores.collect()
labelsAndWeights.sort(key=lambda (k, v): v, reverse=True)
labelsByWeight = np.array([k for (k, v) in labelsAndWeights])

length = labelsByWeight.size
truePositives = labelsByWeight.cumsum()
numPositive = truePositives[-1]
falsePositives = np.arange(1.0, length + 1, 1.) - truePositives

truePositiveRate = truePositives / numPositive
falsePositiveRate = falsePositives / (length - numPositive)

import matplotlib.pyplot as plt
%matplotlib inline

fig = plt.figure()
ax = plt.subplot(111)

ax.set_xlim(-.05, 1.05), ax.set_ylim(-.05, 1.05)
ax.set_ylabel('True Positive Rate (Sensitivity)')
ax.set_xlabel('False Positive Rate (1 - Specificity)')
plt.plot(falsePositiveRate, truePositiveRate, color='#8cbfd0', linestyle='-', linewidth=3.)
plt.plot((0., 1.), (0., 1.), linestyle='--', color='#d6ebf2', linewidth=2.)

f:id:tkzs:20160523053011p:plain


次に元データに対してハッシュトリックを当てます。
OHE後の次元数は20万程度、今回ハッシュ関数のバケット数を2の15乗で32768個としてます。

from collections import defaultdict
import hashlib

def hashFunction(numBuckets, rawFeats, printMapping=False):
    mapping = {}
    for ind, category in rawFeats:
        featureString = category + str(ind)
        mapping[featureString] = int(int(hashlib.md5(featureString).hexdigest(), 16) % numBuckets)
    sparseFeatures = defaultdict(float)
    for bucket in mapping.values():
        sparseFeatures[bucket] += 1.0
    return dict(sparseFeatures)


def parseHashPoint(point, numBuckets):
    parsedPoints = parsePoint(point)
    items = point.split(',')
    label = items[0]
    features = hashFunction(numBuckets, parsedPoints, printMapping=False)
    return LabeledPoint(label, SparseVector(numBuckets, features))

numBucketsCTR = 2 ** 15
hashTrainData = rawTrainData.map(lambda x: parseHashPoint(x, numBucketsCTR))
hashTrainData.cache()
hashValidationData = rawValidationData.map(lambda x: parseHashPoint(x, numBucketsCTR))
hashValidationData.cache()
hashTestData = rawTestData.map(lambda x: parseHashPoint(x, numBucketsCTR))
hashTestData.cache()

print ('Feature size after hushing:\n\tNumber of features = {0}'
       .format(numBucketsCTR))

print hashTrainData.take(1)

Feature size after hushing:
	Number of features = 32768
[LabeledPoint(0.0, (32768,[1305,2883,3807,4814,4866,4913,6952,7117,9985,10316,11512,11722,12365,13893,14735,15816,16198,17761,19274,21604,22256,22563,22785,24855,25202,25533,25721,26487,26656,27668,28211,29152,29402,29873,30039,31484,32493,32708],[1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0]))]

ハッシュ関数をあてた後のデータにてロジスティック分析。
デフォルト設定のものと、グリッドサーチでパラメータを探したものと二つモデルを作ります。

まずは素うどんの方から。

numIters = 100
stepSize = 1.0
regParam = 0.01
regType = 'l2'
includeIntercept = True

model1 = LogisticRegressionWithSGD.train(data = hashTrainData, 
                                         iterations = numIters, 
                                         step = stepSize, 
                                         regParam = regParam,
                                         regType = regType,
                                         intercept = includeIntercept)

logLossValLR1 = evaluateResults(model1, hashValidationData)
print ('Validation Logloss for model1:\n\t = {0}'
       .format(logLossValLR1))

Validation Logloss for model1:
	 = 0.482691122185


次にグリッドサーチを実施してモデル構築します。

numIters = 500
stepSizes = [0.1, 0.5, 1, 3, 5, 7]
regParams = [1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2]
regType = 'l2'
includeIntercept = True


# Initialize variables using values from initial model training
bestModel = None
bestLogLoss = 1e10
logLoss = np.zeros((len(stepSizes),len(regParams)))

for i in xrange(len(stepSizes)):
    for j in xrange(len(regParams)):
        model = LogisticRegressionWithSGD.train(data = hashTrainData, 
                                                iterations = numIters, 
                                                step = stepSizes[i], 
                                                regParam = regParams[j],
                                                regType = regType,
                                                intercept = includeIntercept)
        
        logLoss[i, j] = evaluateResults(model, hashValidationData)
        
        if (logLoss[i, j] < bestLogLoss):
            bestModel = model
            bestStepsize = stepSizes[i]
            bestParams = regParams[j]
            bestLogLoss = logLoss[i,j]
            

print ('best parameters:\n\tBest Stepsize = {0:.3f}\n\tBest Rarams = {1:.3f}'
       .format(bestStepsize, bestParams))

print bestParams,logLoss

%matplotlib inline
from matplotlib.colors import LinearSegmentedColormap


numRows, numCols = len(stepSizes), len(regParams)
logLoss = np.array(logLoss)
logLoss.shape = (numRows, numCols)

fig = plt.figure()
ax = plt.subplot(111)


ax.set_xticklabels(regParams), ax.set_yticklabels(stepSizes)
ax.set_xlabel('Regularization Parameter'), ax.set_ylabel('Step Size')

colors = LinearSegmentedColormap.from_list('blue', ['#0022ff', '#000055'], gamma=.2)
image = plt.imshow(logLoss,interpolation='nearest', aspect='auto',
                    cmap = colors)
pass

best parameters:
	Best Stepsize = 7.000
	Best Rarams = 0.000
1e-07 [[ 0.51910013  0.51910015  0.51910029  0.51910175  0.51911634  0.51926454]
 [ 0.48589729  0.4858974   0.48589847  0.48590921  0.48605802  0.48715278]
 [ 0.47476165  0.47476184  0.47476375  0.47478289  0.4750353   0.47731296]
 [ 0.46130496  0.46130542  0.46131004  0.46137622  0.46208485  0.46799094]
 [ 0.45587263  0.45587339  0.45588094  0.45600577  0.45715016  0.465792  ]
 [ 0.45268179  0.45268281  0.45270706  0.45286577  0.45438559  0.46488834]]

f:id:tkzs:20160523053047p:plain


最期に下記の軸によるLoglosによるモデル評価です。
-OHEデータモデル(model1) vs Hushingデータモデル(model2)
-Hushingデータ・デフォルトモデル(model2) vs Hushingデータ・GSモデル(model3)

Loglossでの比較では、次元量を圧縮したHushing後のデータを使っても予測精度が落ちていませんね、面白いです。

logLossTestLR0 = evaluateResults(model0, OHETestData)
logLossTestLR1 = evaluateResults(model1, hashTestData)
logLossTestLRbest = evaluateResults(bestModel, hashTestData)

print ('OHECoding & Hashed Features Test Logloss:\n\tOHEModel = {0:.3f}\n\thashModel = {1:.3f}'
       .format(logLossTestLR0, logLossTestLR1))

print ('Hashed Features Validation Test Logloss:\n\tBaseModel = {0:.3f}\n\tBestModel = {1:.3f}'
       .format(logLossTestLR1, logLossTestLRbest))

OHECoding & Hashed Features Test Logloss:
	OHEModel = 0.490
	hashModel = 0.490
Hashed Features Validation Test Logloss:
	BaseModel = 0.490
	BestModel = 0.460

時系列解析① ARモデル

時系列データの解析を二つのデータに対して実施します。
基本的には、RawDataの確認、トレンドと季節性の排除、モデルのフィッティングの順にやっていきます。

データ①simulated data

Rで生成したデータについて分析してみます。
データを作成するコードは下記。
データのプロットと、自己相関関数を併せて出しておきます。

#### Generate time seriese data with a linear trend
>time <- 1:200
>data.ar <- arima.sim(model=list(ar=c(0.6, 0.3)), n=200, sd=1)
>data <- data.ar + 30 + 0.05 * time
>par(mfrow=c(2,1))
>plot(time, data, type="l", main="Raw data plot")
>acf(data, main="Sample autocorrelation function")

出てくるデータはこんな感じです。
グラフとコレログラム(Correlogram)から線形トレンドがありそうなので、こいつを取り除きます。

f:id:tkzs:20160217233012p:plain

## Remove the trend
>linear.model <- lm(data~time)
>summary(linear.model)

Call:
lm(formula = data ~ time)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.7143 -1.0238 -0.0246  1.0379  3.7570 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 27.098000   0.207381   130.7   <2e-16 ***
time         0.065674   0.001789    36.7   <2e-16 ***

>residual.series <- data - linear.model$fitted.values
>par(mfrow=c(3,1))
>plot(residual.series, main="Residual series")
>acf(residual.series, main="ACF of Residual series")
>pacf(residual.series, main="PACF of Residual series")

f:id:tkzs:20160217233040p:plain

時系列プロットから定常性が確認でき偏自己相関から、AR(2)が妥当なオーダーだと判断出来ますので、早速モデル化。

## Fit an AR(2) model to the data
model.ar <- arima(residual.series, order=c(2,0,0))

par(mfrow=c(3,1))
plot(model.ar$residuals, main="Residual series")
acf(model.ar$residuals, main="ACF of Residual series")
pacf(model.ar$residuals, main="PACF of Residual series")

ARIMA()関数は残差をプロットしてくれます。
プロットを見る限り、残差に相関はなさそうなのでモデルが妥当だったということでしょう。

f:id:tkzs:20160217233721p:plain

データ②Temperature Data

湖の温度を測ったデータを使っての時系列分析です。
データは私のGitHubに置いてあります。
GitHub:
https://github.com/tkazusa/Sample_code/tree/master/Data


> #Data import
> tempeature <- read.csv("temperature.csv", header=TRUE)
> t <- 1:nrow(tempeature)
> x <- tempeature$Temperature
> par(mfrow=c(2,1))
> plot(t,x, type="l", main="Raw data plot")
> acf(x,main="Temperature data ACF")

図から、先ほどの例とは異なり季節性が出ていることがわかります。

f:id:tkzs:20160218213525p:plain

この季節性を取り除くために、下記のモデルを使います。

s_t = \beta_0 + \beta_1 sin(\frac{2 \pi t}{12}) + \beta_2 cos(\frac{2 \pi t}{12})

> model.season <- lm(x~sin(2*pi*t/365) + cos(2*pi*t/365))
> summary(model.season)

Call:
lm(formula = x ~ sin(2 * pi * t/365) + cos(2 * pi * t/365))

Residuals:
     Min       1Q   Median       3Q      Max 
-10.3755  -2.1802  -0.0533   2.4349   9.0749 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)           6.4263     0.1029   62.47   <2e-16 ***
sin(2 * pi * t/365)  -2.5477     0.1455  -17.51   <2e-16 ***
cos(2 * pi * t/365)  -3.7933     0.1455  -26.07   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.404 on 1092 degrees of freedom
Multiple R-squared:  0.4746,	Adjusted R-squared:  0.4736 
F-statistic: 493.2 on 2 and 1092 DF,  p-value: < 2.2e-16

> residual.series <- x - model.season$fitted.values
> par(mfrow=c(3,1))
> plot(residual.series, main="Residual series", type="l")
> acf(residual.series, main="ACF of Residual series")
> pacf(residual.series, main="PACF of Residual series")

f:id:tkzs:20160218214622p:plain

プロットからResidualデータは定常性があり、PACFからAR(1)モデルが妥当だと判断できますので、モデル化します。

>## Fit an AR(1) model
> model.ar <- arima(residual.series, order=c(1,0,0))
> model.ar

Call:
arima(x = residual.series, order = c(1, 0, 0))

Coefficients:
         ar1  intercept
      0.4861     0.0026
s.e.  0.0264     0.1746

sigma^2 estimated as 8.829:  log likelihood = -2746.36,  aic = 5498.72
> par(mfrow=c(3,1))
> plot(model.ar$residuals, main="Residual series")
> acf(model.ar$residuals, main="ACF of Residual series")
> pacf(model.ar$residuals, main="PACF of Residual series")

f:id:tkzs:20160218220208p:plain

残差プロットやPACFより綺麗にモデル化できていると判断できます。

OLS、リッジ回帰、ラッソでのパフォーマンス比較

OLS、特徴選択、リッジ回帰、ラッソの4つの方法でTrainデータからモデルをFittingし、Testデータを用いて、平均2乗誤差(MSE)を推定して比べるといったことをします。

まずはデータの準備

データはStanford大学の統計学部からprostateのデータをダウンロード、正規化、訓練データ、テストデータに分けて使います。

> # Download prostate data
> con = "http://statweb.stanford.edu/~tibs/ElemStatLearn/datasets/prostate.data"
> prostate=read.csv(con,row.names=1,sep="\t")
> 
> # Scale data and prepare train/test split
> prost.std <- data.frame(cbind(scale(prostate[,1:8]),prostate$lpsa))
> names(prost.std)[9] <- 'lpsa'
> data.train <- prost.std[prostate$train,]
> data.test <- prost.std[!prostate$train,]
> y.test <- data.test$lpsa
> n.train <- nrow(data.train)

OLS

単純な最小二乗法でMSEを求めてみます。

> m.ols <- lm(lpsa ~ . , data=data.train)
> summary(m.ols)

Call:
lm(formula = lpsa ~ ., data = data.train)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.64870 -0.34147 -0.05424  0.44941  1.48675 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.46493    0.08931  27.598  < 2e-16 ***
lcavol       0.67953    0.12663   5.366 1.47e-06 ***
lweight      0.26305    0.09563   2.751  0.00792 ** 
age         -0.14146    0.10134  -1.396  0.16806    
lbph         0.21015    0.10222   2.056  0.04431 *  
svi          0.30520    0.12360   2.469  0.01651 *  
lcp         -0.28849    0.15453  -1.867  0.06697 .  
gleason     -0.02131    0.14525  -0.147  0.88389    
pgg45        0.26696    0.15361   1.738  0.08755 .  
---
Signif. codes:  
0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7123 on 58 degrees of freedom
Multiple R-squared:  0.6944,	Adjusted R-squared:  0.6522 
F-statistic: 16.47 on 8 and 58 DF,  p-value: 2.042e-12

> y.pred.ols <- predict(m.ols,data.test)
> summary((y.pred.ols - y.test)^2)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.000209 0.044030 0.116400 0.521300 0.445000 4.097000 
> 
> mse.ols=sum((y.pred.ols - y.test)^2)
> mse.ols
[1] 15.63822

特徴選択

次に特徴量を減らしてみます。
leap関数でCp(Process Capability Index)をもとに特徴選択を行った上で、MSEを求めます。
whichで確認すると、gleasonを削ったモデルが良さそうだと判断されてます。

> library(leaps)
> l <- leaps(data.train[,1:8],data.train[,9],method='Cp')
> plot(l$size,l$Cp)
> 
> # Select best model according to Cp
> bestfeat <- l$which[which.min(l$Cp),]
> bestfeat
    1     2     3     4     5     6     7     8 
 TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE 
> 
> # Train and test the model on the best subset
> m.bestsubset <- lm(lpsa ~ .,data=data.train[,bestfeat])
> summary(m.bestsubset)

Call:
lm(formula = lpsa ~ ., data = data.train[, bestfeat])

Residuals:
     Min       1Q   Median       3Q      Max 
-1.65425 -0.34471 -0.05615  0.44380  1.48952 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.46687    0.08760  28.161  < 2e-16 ***
lcavol       0.67645    0.12384   5.462 9.88e-07 ***
lweight      0.26528    0.09363   2.833   0.0063 ** 
age         -0.14503    0.09757  -1.486   0.1425    
lbph         0.20953    0.10128   2.069   0.0430 *  
svi          0.30709    0.12190   2.519   0.0145 *  
lcp         -0.28722    0.15300  -1.877   0.0654 .  
pgg45        0.25228    0.11562   2.182   0.0331 *  
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7064 on 59 degrees of freedom
Multiple R-squared:  0.6943,	Adjusted R-squared:  0.658 
F-statistic: 19.14 on 7 and 59 DF,  p-value: 4.496e-13

> y.pred.bestsubset <- predict(m.bestsubset,data.test[,bestfeat])
> summary((y.pred.bestsubset - y.test)^2)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.000575 0.038000 0.122400 0.516500 0.449300 4.002000 
> 
> mse.bestsubset=sum((y.pred.bestsubset - y.test)^2)
> mse.bestsubset
[1] 15.4954

リッジ回帰

次にリッジ回帰でやってみます。
リッジ回帰とラッソについて、なんだっけという場合は下記がご参考。
tjo.hatenablog.com

qiita.com

で、さくっとやってみます。

> library(MASS)
> m.ridge <- lm.ridge(lpsa ~ .,data=data.train, lambda = seq(0,20,0.1))
> plot(m.ridge)
> 
> #by using the select function
> select(m.ridge)
modified HKB estimator is 3.355691 
modified L-W estimator is 3.050708 
smallest value of GCV  at 4.9 

下記のような図と共に、ベストなλは4.9だと選ばれました。
f:id:tkzs:20160126151347j:plain

次に予測値を返させて、MSEを求めます。

> y.pred.ridge = scale(data.test[,1:8],center = m.ridge$xm, scale = m.ridge$scales)%*% m.ridge$coef[,which.min(m.ridge$GCV)] + m.ridge$ym
> summary((y.pred.ridge - y.test)^2)
       V1          
 Min.   :0.000003  
 1st Qu.:0.043863  
 Median :0.136450  
 Mean   :0.494410  
 3rd Qu.:0.443495  
 Max.   :3.809051  
> mse.ridge=sum((y.pred.ridge - y.test)^2)
> mse.ridge
[1] 14.8323

ラッソ

あまり使っている見ませんが、larsライブラリー(https://cran.r-project.org/web/packages/lars/index.html)を使います。

> library(lars)
> m.lasso <- lars(as.matrix(data.train[,1:8]),data.train$lpsa,type="lasso")
> plot(m.lasso)
> # Cross-validation
> r <- cv.lars(as.matrix(data.train[,1:8]),data.train$lpsa)
> bestfraction <- r$index[which.min(r$cv)]
> bestfraction
[1] 0.5252525

f:id:tkzs:20160126210507j:plain

fractionは0.525あたりと図と出力から判断出来ます。

> # Observe coefficients
> coef.lasso <- predict(m.lasso,as.matrix(data.test[,1:8]),s=bestfraction,type="coefficient",mode="fraction")
> 
> # Prediction
> y.pred.lasso <- predict(m.lasso,as.matrix(data.test[,1:8]),s=bestfraction,type="fit",mode="fraction")$fit
> summary((y.pred.lasso - y.test)^2)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.000144 0.072470 0.128000 0.454000 0.394400 3.565000 
> 
> mse.lasso=sum((y.pred.lasso - y.test)^2)
> mse.lasso
[1] 13.61874
> 

平均二乗誤差の比較

> # Compare them 
> mse.ols
[1] 15.63822
> mse.bestsubset
[1] 15.4954
> mse.ridge
[1] 14.8323
> mse.lasso
[1] 13.61874

今回のデータに関しては、ラッソで行った分析が一番MSEを最小に出来ました。

PRML 7章 サポートベクトルマシン

パターン認識と機械学習、第7章「Sparse Kernel Machines」に入ってSVMの実装をしてみます。
図7.2の再現を行うにあたり、カーネルがガウジアンカーネルが選択されているようなので、今回の実装もこれに従います。

いつも通り境界面の求め方など、SVMそのものの解説はPRMLやはじパタにお任せするとして、実装にあたっての流れだけ簡単に示します。

実装の大まかな流れ

①ラグランジュ関数の双対表現 (7.10)について、二次計画法で解く。


\tilde{L}(a) = \sum_{n=1}^N a_n - \frac{1}{2} \sum_{n=1}^N \sum_{m=1}^N  a_na_m t_n t_m k({\bf x_n}, {\bf x_m})    (7.10)


非線形計画法の自前での実装はムリ!って言われたので、おとなしくソルバーを使いました。無償かつポピュラーっぽいcvxoptというライブラリーが良さ気なのでインストール。

f:id:tkzs:20150916005515p:plain

上記のようにドキュメントに従ってcvxoptでモデル化するには(7.10)のままだと見づらいので式変形。


 - \tilde{L}(a) = \frac{1}{2} {\bf A}^{\mathrm{T}} ({\bf T}^{\mathrm{T}} k({\bf x_n}, {\bf x_m}) {\bf T}){\bf A} - {\bf A} (7.10)'

表記的には気持ち悪い形ですが、ひとまず(7.10)'とcvxoptを見比べながらモデル化できたらそのままcvxoptに解いてもらうとラグランジュ乗数aが求まります。

②SVMでは境界面の決定にサポートベクトルのみ関与し、これはKKT条件で絞り込まれます。
なんとなくまだ腹落ちしてないところもありますが、境界面上にあるベクトルはt_ny({\bf x})-1=0になるので、(7.16)より、a>0でないといけないと解釈。これを利用してサポートベクトルを抽出します。

③で、求まったサポートベクトルを使って(7.18)のbを求めます。

④(7.13)に求めた値を代入して、境界面を求めます。

コード

今回、図を描くのにえらい苦労しまいた。どうしても境界面が上手くプロットできず、aidiaryさんのこちらの記事を全面的に参考にさせて頂きました。ありがとうございました。


import math
import matplotlib.pyplot as plt
from scipy import stats
from scipy.stats.kde import gaussian_kde
from pylab import *
import numpy as np
import random
import cvxopt
from cvxopt import matrix, solvers
%matplotlib inline

def kernel(x, y):
    sigma = 5.0
    return np.exp(-norm(x-y)**2 / (2 * (sigma ** 2)))

#(7.10)' (Quadratic Programming)
def L(t, X, N):
    K = np.zeros((N, N))
    for i in range(N):
        for j in range(N):
            K[i, j] = t[i] * t[j] * kernel(X[i], X[j])
    Q = matrix(K)
    p = matrix(-1 * np.ones(N))             
    G = matrix(np.diag([-1.0]*N))       
    h = matrix(np.zeros(N))             
    A = matrix(t, (1,N))                
    b = matrix(0.0)                     
    sol = solvers.qp(Q, p, G, h, A, b)
    a = array(sol['x']).reshape(N)
    return a

#(7.13)
def y_x(a, t, X, N, b, x):
    sum = 0
    for n in range(N):
        sum += a[n] * t[n] * kernel(x, X[n])
    return sum + b

#(7.18)
def b(a, t, X, S):
    sum_A = 0
    for n in S:
        sum_B = 0
        for m in S:
            sum_B += a[m] * t[m] * kernel(X[n], X[m])
        sum_A += (t[n] - sum_B)
    return sum_A/len(S)

if __name__ == "__main__":
    N = 36
    mu_blue = [1,-1]
    cov = [[0.1,0.05], [0.05,0.1]]
    
    x_blue,y_blue = np.random.multivariate_normal(mu_blue, cov, N/2).T
    
    x_red = [0.3, 0.8, 0.9, 0.95, 1.1, 1.3, 1.6, 1.9, 1.75, 1.8, 2.0, 2.1, 2.3, 2.25, 2.4, 2.7, 3.0, 3.2]
    y_red = [-0.2, 0.1, 0.25, 0.14, -0.1, 1.6, 1.2, 0.6, 0.8, -0.6, -0.8, -0.75, 1.2, -1.15, -0.12, -0.3, -0.4, 1.4]
    
    t_blue = np.ones((1, N/2))
    t_red = -1*np.ones((1, N/2))

    blue = vstack((x_blue, y_blue))
    red = vstack((x_red, y_red))

    X = np.concatenate((blue, red), axis=1).T
    t = np.concatenate((t_blue, t_red), axis=1).T
    
    #(7.10)' (Quadratic Programming)
    a = L(t, X, N)

    #Extract Index of support vectors from (7.14) 
    S = []
    for n in range(len(a)):
        if a[n] < 0.0001: continue
        S.append(n)

    #(7.18)
    b = b(a, t, X, S)

    
    #Plot train data sets
    plt.scatter(x_blue,y_blue,color='b',marker='x')
    plt.scatter(x_red,y_red,color='r',marker='x')
    
    # Enphasize suport vectors
    for n in S:
        plt.scatter(X[n,0], X[n,1], color='g', marker='o')
    
    # Plot the decision surface
    X1, X2 = meshgrid(linspace(-10,10,100), linspace(-10,10,100))
    w, h = X1.shape
    X1.resize(X1.size)
    X2.resize(X2.size)
    Z = array([y_x(a, t, X, N, b, array([x1,x2])) for (x1, x2) in zip(X1, X2)])
    X1.resize((w, h))
    X2.resize((w, h))
    Z.resize((w, h))
    CS = contour(X1, X2, Z, [0.0], colors='k', linewidths=1, origin='lower')
    xlim(0, 4)
    ylim(-2, 2)
    title("Figure 7.2")

結果

f:id:tkzs:20150916011816p:plain

【論文メモ】人間を理解するためのロボティクス

石井加代子 (2007). 人間を理解するためのロボティクス. Cognitive studies 14(1), 11-30, 日本認知科学会 (https://www.jstage.jst.go.jp/article/jcss/14/1/14_1_11/_pdf)

「全ての作業をプログラムす ることは不可能であるから,ロボットが作業を習得 する必要があり,ヒトと同じ形態のロボットのほう が,遠隔操作・教育し易いから」

「人類の進化に大きく貢献したのは,学習する能力では なく,教える能力である」と考えられることから, ⃝1 ロボットはまず,“ 学習する ”のでなく,“ 養育・ 教育される ”必要がある,⃝2 ロボットが他のロボッ トに,情報を単に移転するのではなく,“ 教えるこ とができる ”ようになる必要がある,⃝3 他のロボッ トに“ 教えることができる ”ということを,当該ロ ボットの“ 学習が成立した ”状態の証しとする,⃝4 実用目的の明確な研究でなければならない.(p.14)

→逆強化学習に加えて、ロボットによる教育への発展性について。

【発見的(ヒューリスティク)知識】ヒトは, 現実世界の問題を解決するために,時々刻々必 要な有限個の情報を迅速に選択している.機 械にはこれができない〈フレーム問題〉.ヒト も新規で複雑な状況では,ヒューリスティク な解決が困難になるが,フレーム問題が“ な いかのように振舞い ”,立ち往生しないです ますことはできる.(p.16)

→フレームそのものの選択の仕方と、そのフレームを評価/決定する時間の間に環境が変化してしまうことについてどう解決しようか

①まず養育者が笑いかける.他者を真似するという 生得的デバイスのある赤ん坊は,これを真似る.す ると,行動に誘発された快感を体験する.快感という報酬ゆえに,笑う行動を繰り返すようになる. ② 養育者は,赤ん坊の笑いを「赤ん坊自身の情動の 表現である(例:何かが面白くて笑っている)」か のように,応答する.養育者と赤ん坊の間で,笑い という行為が交わされる.養育者と赤ん坊以外の事 物をめぐって(指差し・掌握・揺り動かしなど)この相互作用が続行する.養育者は,子どもが言葉を 理解する前から,快感や笑いに関する様々な言語表 現を赤ん坊に対して語りかけ,赤ん坊の動きや発声 を,「赤ん坊が言語表現とその意味を理解したかの ように」応答することによって,赤ん坊を言語ゲー ムに引き込んでゆく.他者の表情の知覚(知覚モデ ル)と,自分の情動表現(行動モデル)が緊密な関 係で形成される.
③ 他者の笑う顔をみると,自分の行動モデルの逆モ デルを使って,他者が快いと感じているのだと推測 する(他者の感情理解).言語ゲームを介して,こ れに「笑う」・[面白い]・「楽しい」・「嬉しい」など の言語表現が結合される.
④自分が笑ったとき(行動モデル),他者理解(知 覚モデル)とそれに対する言語表現を利用して,「自 分は楽しくて笑っている」というような自己の感情 理解が成立する. (p.23)

→実験心理学における感情の獲得過程。行動の結果に報酬を与えるだけでなく、行動そのものからも報酬を得る必要がある?

集団の利益を最適化するアルゴリズムは,ある特 定の個人が道徳的「悪をなす」ことを禁ずるもので はなく,多くの人が「道徳的悪をされる」ことを抑制するものである.他者の「悪をされたくない」という心的内容を推測できない,あるいは推測や理
解に基づいて行動制御できないならば,特定個人の 「悪をなす」行為は,自己制御できない.協調行動をシミュレートするロボット集団で,心の理論と倫
理的行動の関係を調べることができるだろう. (p.24)

→「自分がされて嫌なことを他人にしない」が集団の利益最大化のアルゴリズム。禁止項目をリスト化するんじゃダメ。

PythonでXGBoostをちゃんと理解する(3) hyperoptでパラメーターチューニング

xgboostのハイパーパラメーターを調整するのに、何が良さ気かって調べると、結局「hyperopt」に落ち着きそう。
対抗馬はSpearmintになりそうだけど、遅いだとか、他のXGBoost以外のモデルで上手く調整できなかった例があるとかって情報もあって、時間の無い今はイマイチ踏み込む勇気はない。

Hyperparameter Optimization using Hyperopt - Otto Group Product Classification Challenge | Kaggle
Optimizing hyperparams with hyperopt - FastML

前回辺りにアルゴリズム振り返って、チューニングには特別気をつけなきゃいけないことも無さそうなので、ガリガリとコード書いて動かしてみます。

hyperoptはつまるところ最適化問題のソルバーで、目的関数を一定の条件下で最大or最小化するというもの。
もともとはベイズ的最適化のためにデザインされてたけど、今はランダムサーチとTree of Parzen Estimators (TPE)だけがインプリされてる状況。

hyperopt本体の挙動を視覚的に理解するには、下記の記事がオススメ。districtdatalabs.silvrback.com

で、id:puyokwさんも記事の中で挙げてらっしゃったけど、KaggleのOttoのコンペにXGBoostにhyperoptを使ったコードが公開されてました。

xgboost package のR とpython の違い - puyokwの日記
optimizing hyperparameters of an xgboost model on otto dataset · bamine/Kaggle-stuff@17f7828 · GitHub

で、もうやりたかったことこれだよね、ってpuyokwさんと丸かぶりなアイデアだったので参考と言うなのほぼ写経。ありがとうございます。
元のコードに交差検定加えたいなと思ったけど、xgb.CVは使い勝手悪いし、sklearnでやってます。

import numpy as np
import xgboost as xgb
import pandas as pd

from sklearn import datasets, cross_validation
from sklearn.cross_validation import train_test_split
from sklearn.metrics import confusion_matrix, log_loss
from sklearn.cross_validation import KFold
from sklearn import preprocessing

from hyperopt import fmin, tpe, hp, STATUS_OK, Trials

iris=datasets.load_iris()
X = iris.data
y = iris.target

#Split data set to train and test data
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.5, random_state=25)

np.random.seed(25)


def score(params):
    print "Training with params : "
    print params
    N_boost_round=[]
    Score=[]
    skf = cross_validation.StratifiedKFold(y_train, n_folds=10, shuffle=True,random_state=25)
    for train, test in skf:
        X_Train, X_Test, y_Train, y_Test = X_train[train], X_train[test], y_train[train], y_train[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)
        N = model.best_iteration
        N_boost_round.append(N)
        score = model.best_score
        Score.append(score)
    Average_best_num_boost_round = np.average(N_boost_round)
    Average_best_score = np.average(Score)
    print "\tAverage of best iteration {0}\n".format(Average_best_num_boost_round)
    print "\tScore {0}\n\n".format(Average_best_score)
    return {'loss': Average_best_score, 'status': STATUS_OK}


def optimize(trials):
    space = {
        "objective": "multi:softprob",
        "eval_metric": "mlogloss",
        
        #Control complexity of model
        "eta" : hp.quniform("eta", 0.2, 0.6, 0.05),
        "max_depth" : hp.quniform("max_depth", 1, 10, 1),
        "min_child_weight" : hp.quniform('min_child_weight', 1, 10, 1),
        'gamma' : hp.quniform('gamma', 0, 1, 0.05),
        
        #Improve noise robustness 
        "subsample" : hp.quniform('subsample', 0.5, 1, 0.05),
        "colsample_bytree" : hp.quniform('colsample_bytree', 0.5, 1, 0.05),
        
        'num_class' : 3,
        'silent' : 1}
    best = fmin(score, space, algo=tpe.suggest, trials=trials, max_evals=250)
    print "best parameters",best

trials = Trials()
optimize(trials)

#出力
#best parameters {'colsample_bytree': 0.9500000000000001, 'min_child_weight': 1.0, #'subsample':0.9500000000000001, 'eta': 0.6000000000000001, 'max_depth': 3.0, 'gamma': 0.0}

#Adapt best params
params={'objective': 'multi:softprob',
        'eval_metric': 'mlogloss',
        'colsample_bytree': 0.9500000000000001, 
        'min_child_weight': 1.0, 
        'subsample': 0.9500000000000001, 
        'eta': 0.6000000000000001, 
        'max_depth': 3.0, 
        'gamma': 0.0,
        'num_class': 3
        }

score(params)

#出力
#Training with params : 
#{'subsample': 0.9500000000000001, 'eta': 0.6000000000000001, 'colsample_bytree': 0.9500000000000001, 'gamma': #0.0, 'eval_metric': 'mlogloss', 'objective': 'multi:softprob', 'num_class': 3, 'max_depth': 3.0, #'min_child_weight': 1.0}
#	Average of best iteration 37.7
#	Score 0.1092628

X_train = pd.DataFrame(X_train)
X_test = pd.DataFrame(X_test)
dtrain = xgb.DMatrix(X_train.as_matrix(),label=y_train.tolist())
dtest=xgb.DMatrix(X_test.as_matrix())

bst=xgb.train(params,dtrain,num_boost_round=38)
pred=bst.predict(dtest)
pred=pd.DataFrame(pred)
print confusion_matrix(y_test, pred.idxmax(axis=1))

#出力:Confusion_matrix
#[[26  0  0]
# [ 0 24  2]
# [ 0  2 21]]

ちなみに、パラメーターを全部デフォルトにした素うどんXGBoostでやってみても、結果は一緒でした。
あやめのデータだとこんなもんですかね。

params_Default={'objective': 'multi:softprob',
        'eval_metric': 'mlogloss',
        'eta': 0.3, 
        'max_depth': 6, 
        'min_child_weight': 1, 
        'subsample': 1, 
        'colsample_bytree': 1, 
        'num_class': 3,
        }
cv=xgb.cv(params_Default,dtrain,num_boost_round=200,nfold=10)
print(cv)

bst=xgb.train(params_Default,dtrain,num_boost_round=13)
pred=bst.predict(dtest)
pred=pd.DataFrame(pred)
print confusion_matrix(y_test, pred.idxmax(axis=1))

#出力:Confusion_matrix
#[[26  0  0]
# [ 0 24  2]
# [ 0  2 21]]

最近、やりたいことはすでにほかの人がやっちゃってるパターン多くて恥ずかしい限りですね。
本当はこの後、xgboostで得た情報をどう特徴選択などなどでどう使うかまでやろうかと思ったけど、そのネタはアンサンブルとして扱った方が良いのかなと思いxgboostはとりあえず一旦終了。

次は、Spermintとhyperoptの比較あたりやってみようかな。
ChainerとかLasagneとかのニューラルネット周りもやりたいけど。

PythonでXGBoostをちゃんと理解する(2) ライブラリ作者から学ぶ

XGBoostについて調べてたら、開発者本人から学ぶ的な動画があったので観てみた。

www.youtube.com

時間にして約1時間半、英語が苦手でなくて時間がある方は直接見て頂くと面白いかも。

目次はこんな感じ。
・Introduction
・Basic Walkthrough
・Real World Application
・Model Specification
・Parameter Introduction
・Advanced Features
・Gaggle Winning Solution


前半は「まぁみんな知ってるよね」ってことが多かったが、Model SpecificationとParameter Introductionの中のParameter Tuning、あとAdvanced Featuresが面白かったのでメモ。

Model Specification

よく見る勾配ブースティングの説明とは違うように見えたので、面白くて動画を見ながらスライドの内容をメモってみました。
所詮スライドなので、動画内での口頭での補足があるとよりわかりやすいですが、ざっと眺めてもそれなりにわかりやすいのでご参考にどうぞ。
日本語訳が雑でごめんなさい。

動画内スライド引用 (24:00〜)

K本の決定木があったとして、f_kがその木からの予測値とすると、モデルは\sum_{k=1}^K f_k。(P31)
 i番目のデータ、x_iの予測値は、\hat{y_i} = \sum_{k=1}^K f_k (x_i)であり、同様にt番目のステップでの予測値は、\hat{y_i}^{(t)} = \sum_{k=1}^t f_k (x_i)となる(p32)。損失関数Lは、最小二乗法、2値分類と多クラス分類のためなどそれぞれ目的に応じて使う。

また過学習を防ぐために正則化項を導入し、その式は

 \Omega = \gamma {\bf T} +\frac{1}{2}  \lambda\sum_{j=1}^T w_j^2
ここで、{\bf T}は葉の数、w_j^2はj番目の葉における重みの値である(P34)。
損失関数と正則化項をまとめると下記のようになる。

 Obj = L + \Omega

XGboostでは、このObjを最適化するために、勾配降下法を用いる。

目的関数を下記のように置き直し、

 {Obj}^{(t)} = \sum_{i=1}^N L(y_i, \hat{y_i}^{(t)}) + \sum_{i=1}^t \Omega(f_i) = \sum_{i=1}^N L(y_i, \hat{y_i}^{(t-1)} + f_t(x_t)) + \sum_{i=1}^t \Omega(f_i)


1次導関数、2次導関数それぞれをそれぞれ下記のように求めると(P37)、

{\partial_{\hat{y_i}^{(t)}} Obj^{(t)}}

{\partial_{\hat{y_i}^{(t)}}^2 Obj^{(t)}}

テイラー展開することで、

 {Obj}^{(t)} \simeq \sum_{i=1}^N [L(y_i, \hat{y_i}^{(t-1)}) + \frac{1}{2} h_i f_t^2(x_i)] + \sum_{i=1}^t \Omega(f_i)
と近似できる。ここでは、
g_i =\partial_{\hat{y_i}^{(t)}} l(y_i, \hat{y_i}^{(t-1)})
h_i =\partial_{\hat{y_i}^{(t)}}^2 l(y_i, \hat{y_i}^{(t-1)})
である(P38)。で、この目的関数{Obj}^{(t)}を最適化させてf_tを求めることがゴール。

じゃあこのf_tを出力するために、①決定木の形をどうやって求める?&②予測の値をどうやって出力させる?を考える(P43)。

①の問題が解決できてるとして、決定木は、f_t(x) = w_{q(x)}で表される。ここでq(x)は、q(x)番目の葉にデータを当てはめる関数である。この式は予測値を返す過程を、
・データxを葉qに割り当てる
・対応する重みw_{q(x)}を割り当てる

j番目の葉に割り当てられたデータにインデックスをI_j = \{i|q(x_i) = j\}に従って割り振る。

ここで、目的関数Objを書き換えると、
Obj^{(t)} = \sum_{i=1}^n [g_i f_i(x_i) + \frac{1}{2} h_i f_t^2(x_i)] + \gamma {\bf T} +\frac{1}{2} \lambda\sum_{j=1}^T w_j^2  =\sum_{j=1}^T [(\sum_{i \in I_j} g_i) w_i +  \frac{1}{2} (\sum_{i \in I_j} h_i+\lambda )]{w_i}^2 + \gamma {\bf T}

同じ葉にいる全てのデータは同じ予測値を共有するので、この式は葉ごとの予測値の総和を返す。

ここでw_jは2次方程式となるので、Objを最適化するw_jは簡単に求まり、

w_j^{\ast} = - \frac{\sum_{i \in I_j} g_i}{\sum_{i \in I_j} h_i+\lambda}

となるから、それに応じたObjの値は下記となる(P48)。
 Obj^{(t)} = - \frac{1}{2} \sum_{j = 1}^T  \frac{(\sum_{i \in I_j} g_i)^2}{\sum_{i \in I_j} h_i+\lambda}+ \gamma {\bf T}

これで重みの値が求まるので、次にどのように木の形を決めるかを考える。
通常は木の決め方をジニ係数などでGainを計算するが、勾配ブースティングでは、


gain = \frac{1}{2} [\frac{(\sum_{i \in I_L} g_i)^2}{\sum_{i \in I_L} h_i+\lambda} + \frac{(\sum_{i \in I_R} g_i)^2}{\sum_{i \in I_R} h_i+\lambda}- \frac{(\sum_{i \in I} g_i)^2}{\sum_{i \in I} h_i+\lambda} - \gamma ]
として、gainを計算する(P53)。
この式に従って木をMax_depthに達するまで茂らせていき、gainが負になる枝を剪定する。
この作業を指定の回数繰り返す。

Parameter Tunin(55:20〜)&Advanced Features(1:00:00〜)

パラメーターをチューニングするにあたっては、以下の3つに気をつけなさいと。

①Control Overfitting
過学習を抑えるために、パラメーターについては二つに分けて考えているようです。
 ・モデルの複雑さをコントロールするために
   −max_depth, min_child_weight,gamma
 ・ノイズに対する堅牢さを高めるために
  −sub_sample, colsample_bytree
まぁこの辺りは何も特別な事考えず、etaを加えて6個のパラメーターを結果を見ながら調整するんでいいんじゃないかなと思います。

②Deal with Imbalanced data
データのバランスをコントロールするために以下3つを試してみなさいと。
・scale_pos_weight使ってデータのバランスを調整
・AUCを評価指標に使う
・max_delta_stepを"1"辺りに設定する


③Trust the cross validation
動画の中でも、CVを信じなさいと言ってました。
で、CVを行ってnroundを決める場合には、十分な回数を指定しつつearly.stop.roundを使って終わらせる方法が取り上げられてました。
さらに、Overfittingが見られたら、etaを少なくして、nroundを大きくすることを同時にする。
また、個別のオススメの機能としては、xgb.cv()で"prediction = TRUE"にすることで、予測値のベクトルを返してくれるようになります。

長くなったので、「じゃあパラメーターチューニングは具体的にどうやる?」については次回にします。