時系列解析① 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"にすることで、予測値のベクトルを返してくれるようになります。

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

Pythonで実装 PRML 第5章 ニューラルネットワーク

PRML5章から図5.3を再現するために、ニューラルネットワークを実装してみます。
先に申し上げておきますと、偉そうに実装とか言ってるものの、コードから再現した図は歯切れの悪いものとなっております。。。

まず、図5.3(b),(c),(d)に関してはPRMLの中の図に比べて予測精度がいまいち良くない印象。さらに、図5.3(a)に関しては全く見当はずれな予測が返ってくるという状況。試行錯誤しましたが、力不足でして、どなたか間違い気づかれましたらご指摘下さい。

ニューラルネットワークや、誤差伝播法(Backpropagation)そのものの解説はPRMLやはじパタなどに任せるとして、実装に必要な部分だけざっと確認したいと思います。

実装の大まかな流れ

①ニューラルネットを経たアウトプットは(5.9)で表される。PRML文中の式は活性化関数h()にシグモイド関数を想定しているが、図5.3ではtanh()が指定されている点に注意。


 y_k({\bf x}, {\bf w}) = \sigma(\sum_{j=0}^M, w^{(2)}_{kj} h(\sum_{i=0}^D w^{(1)}_{ji} x_i)) (5.9)

②ノード間の重み{\bf w}を学習するにあたり、各ノードでのアウトプットと実測値との差を求める。まずは隠れユニットの出力は(5.63)、出力ユニットの出力は(5.64)。


 z_j = {\rm tanh} (\sigma(\sum_{i=0}^D, w^{(1)}_{ji} x_i)) (5.63)


 y_k = \sum_{j=0}^M, w^{(2)}_{kj} z_i (5.64)


③次に出力層での誤差\delta_kを求める。

\delta_k = y_k - t_k (5.65)

④次に隠れ層での誤差\delta_jを求める。

\delta_j = (1-{z_j}^2) \sum_{k=1}^K w_{kj} \delta_k (5.65)

⑤(5.43)、(5.67)を用いてノード間の重みを更新する。


{\bf w}^{\rm \tau+1} = {\bf w}^{\rm \tau} - \mu \nabla  E({\bf{w}})(5.43)

コード

import matplotlib.pyplot as plt
from pylab import *
import numpy as np
import random

def heaviside(x):
    return 0.5 * (np.sign(x) + 1)

def NN(x_train, t, n_imput, n_hidden, n_output, eta, W1, W2, n_loop):
    for n in xrange(n_loop):
        for n in range(len(x_train)):
            x = np.array([x_train[n]])
            
            #feedforward
            X = np.insert(x, 0, 1) #Insert fixed term

            A = np.dot(W1, X) #(5.62)
            Z = np.tanh(A)  #(5.63)
            Z[0] = 1.0
            Y = np.dot(W2, Z) #(5.64)

   
            #Backprobagation
            D2 = Y - t[n]#(5.65)
            D1 = (1-Z**2)*W2*D2 #(5.66)
    
            W1 = W1- eta*D1.T*X #(5.67), (5.43)
            W2 = W2- eta*D2.T*Z #(5.67), (5.43)
    return  W1, W2

def output(x, W1, W2):
    X = np.insert(x, 0, 1) #Insert fixed term
            
    A = np.dot(W1, X) #(5.62)
    Z = np.tanh(A)  #(5.63)
    Z[0] = 1.0 #Insert fixed term
    Y = np.dot(W2, Z) #(5.64)
    return Y, Z

if __name__ == "__main__":
    #Set form of nueral network 
    n_imput = 2
    n_hidden = 4
    n_output = 1
    eta = 0.1
    W1 = np.random.random((n_hidden, n_imput))
    W2 = np.random.random((n_output, n_hidden))
    n_loop = 1000
    
    
    #Set train data
    x_train = np.linspace(-4, 4, 300).reshape(300, 1)
    y_train_1 = x_train * x_train
    y_train_2 = np.sin(x_train)
    y_train_3 = np.abs(x_train)
    y_train_4 = heaviside(x_train)

    #x_real = np.linspace(-4, 4, 2000).reshape(2000, 1)
    
    W1_1, W2_1= NN(x_train, y_train_1, n_imput, n_hidden, n_output, eta, W1, W2, n_loop) 
    W1_2, W2_2= NN(x_train, y_train_2, n_imput, n_hidden, n_output, eta, W1, W2, n_loop)
    W1_3, W2_3= NN(x_train, y_train_3, n_imput, n_hidden, n_output, eta, W1, W2, n_loop)
    W1_4, W2_4= NN(x_train, y_train_4, n_imput, n_hidden, n_output, eta, W1, W2, n_loop)

    Y_1 = np.zeros((len(x_train), n_output))
    Z_1 = np.zeros((len(x_train), n_hidden))

    Y_2 = np.zeros((len(x_train), n_output))
    Z_2 = np.zeros((len(x_train), n_hidden))

    Y_3 = np.zeros((len(x_train), n_output))
    Z_3 = np.zeros((len(x_train), n_hidden))

    Y_4 = np.zeros((len(x_train), n_output))
    Z_4 = np.zeros((len(x_train), n_hidden))

    for n in range(len(x_train)):
        Y_1[n], Z_1[n] =output(x_train[n], W1_1, W2_1)
        Y_2[n], Z_2[n] =output(x_train[n], W1_2, W2_2)
        Y_3[n], Z_3[n] =output(x_train[n], W1_3, W2_3)
        Y_4[n], Z_4[n] =output(x_train[n], W1_4, W2_4)
    
    
    plt.plot(x_train, Y_1, "r-")
    plt.plot(x_train, y_train_1, "bo", markersize=3)
    for i in range(n_hidden):
        plt.plot(x_train, Z_1[:,i], 'm--')
    xlim([-1,1])
    ylim([0, 1])
    title("Figure 5.3(a)")
    show()
    
    plt.plot(x_train, Y_2, "r-")
    plt.plot(x_train, y_train_2, "bo", markersize=2)
    for i in range(n_hidden):
        plt.plot(x_train, Z_2[:,i], 'm--')
    xlim([-3.14,3.14])
    ylim([-1, 1])
    title("Figure 5.3(b)")
    show()
    
    
    plt.plot(x_train, Y_3, "r-")
    plt.plot(x_train, y_train_3, "bo", markersize=4)
    for i in range(n_hidden):
        plt.plot(x_train, Z_3[:,i], 'm--')
    xlim([-1,1])
    ylim([0, 1])
    title("Figure 5.3(c)")
    show()
    
    
    plt.plot(x_train, Y_4, "r-")
    plt.plot(x_train, y_train_4, "bo" ,markersize=2)
    for i in range(n_hidden):
        plt.plot(x_train, Z_4[:,i], 'm--')
    xlim([-2,2])
    ylim([-0.05, 1.05])
    title("Figure 5.3(d)")
    show()

結果

f:id:tkzs:20151003070323p:plainf:id:tkzs:20151003070332p:plainf:id:tkzs:20151003070334p:plainf:id:tkzs:20151003070341p:plain