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.)
次に元データに対してハッシュトリックを当てます。
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]]
最期に下記の軸による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)から線形トレンドがありそうなので、こいつを取り除きます。
## 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")
時系列プロットから定常性が確認でき偏自己相関から、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()関数は残差をプロットしてくれます。
プロットを見る限り、残差に相関はなさそうなのでモデルが妥当だったということでしょう。
データ②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")
図から、先ほどの例とは異なり季節性が出ていることがわかります。
この季節性を取り除くために、下記のモデルを使います。
> 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")
プロットから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")
残差プロットや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
で、さくっとやってみます。
> 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だと選ばれました。
次に予測値を返させて、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
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)について、二次計画法で解く。
非線形計画法の自前での実装はムリ!って言われたので、おとなしくソルバーを使いました。無償かつポピュラーっぽいcvxoptというライブラリーが良さ気なのでインストール。
上記のようにドキュメントに従ってcvxoptでモデル化するには(7.10)のままだと見づらいので式変形。
表記的には気持ち悪い形ですが、ひとまず(7.10)'とcvxoptを見比べながらモデル化できたらそのままcvxoptに解いてもらうとラグランジュ乗数が求まります。
②SVMでは境界面の決定にサポートベクトルのみ関与し、これはKKT条件で絞り込まれます。
なんとなくまだ腹落ちしてないところもありますが、境界面上にあるベクトルはになるので、(7.16)より、でないといけないと解釈。これを利用してサポートベクトルを抽出します。
③で、求まったサポートベクトルを使って(7.18)のを求めます。
④(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")
結果
【論文メモ】人間を理解するためのロボティクス
石井加代子 (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について調べてたら、開発者本人から学ぶ的な動画があったので観てみた。
時間にして約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〜)
本の決定木があったとして、がその木からの予測値とすると、モデルは。(P31)
番目のデータ、の予測値は、であり、同様に番目のステップでの予測値は、となる(p32)。損失関数は、最小二乗法、2値分類と多クラス分類のためなどそれぞれ目的に応じて使う。また過学習を防ぐために正則化項を導入し、その式は
ここで、は葉の数、はj番目の葉における重みの値である(P34)。
損失関数と正則化項をまとめると下記のようになる。
XGboostでは、このを最適化するために、勾配降下法を用いる。
目的関数を下記のように置き直し、
1次導関数、2次導関数それぞれをそれぞれ下記のように求めると(P37)、
テイラー展開することで、
と近似できる。ここでは、
である(P38)。で、この目的関数を最適化させてを求めることがゴール。じゃあこのを出力するために、①決定木の形をどうやって求める?&②予測の値をどうやって出力させる?を考える(P43)。
①の問題が解決できてるとして、決定木は、で表される。ここでは、番目の葉にデータを当てはめる関数である。この式は予測値を返す過程を、
・データを葉に割り当てる
・対応する重みを割り当てるj番目の葉に割り当てられたデータにインデックスをに従って割り振る。
ここで、目的関数を書き換えると、
同じ葉にいる全てのデータは同じ予測値を共有するので、この式は葉ごとの予測値の総和を返す。
ここでは2次方程式となるので、を最適化するは簡単に求まり、
となるから、それに応じたの値は下記となる(P48)。
これで重みの値が求まるので、次にどのように木の形を決めるかを考える。
通常は木の決め方をジニ係数などでGainを計算するが、勾配ブースティングでは、
として、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"にすることで、予測値のベクトルを返してくれるようになります。
長くなったので、「じゃあパラメーターチューニングは具体的にどうやる?」については次回にします。