kamonohashi – ページ 9 – かものはしの分析ブログ

Japan.R 2016のスライドまとめ

まだ手に入れていないスライドもあるので随時更新しますが、Japan.R 2016(connpass)のスライドをまとめています。後日、登場したパッケージなどのサンプルコードも載せていく予定です。

目次

・石田基広さんのキーノート
・ホクソエムとは何だったのか（ホクソエムさん）
・Rと探索的データ分析で、国連での日本の立ち位置を可視化する（安田洋介さん）
・マウス操作でかんたん予測分析（鈴木了太さん）
・高速・省メモリにlibsvm形式でダンプする方法を研究してみた（@hskksk）
・Rでてんしょくかつどう（@Med_KU）
・RStudio vs Emacs（@y__mattu）
・randomforestで高次元の変数重要度見る（@siero5335）
・Rで本を作りたい（前田和寛さん）
・28歳でプログラミングを始めた話（市川太祐さん）
・LDA-Visパッケージのご紹介（@doradora09）
・【e2d3R】E2D3からDot-Bar-Chartのご紹介（楠本一哲さん）
・このIRのグラフがすごい！上場企業2016（@ito_yan）
・Rでカルマンフィルタをしたい（@tetsuroito）
・PPAP（仮）（@yutannihilation）
・スライド未公開、ユーザーの状態遷移に関する分析のお話（@sanoche16）
・私とR（高栁慎一さん）
・めくってもめくってもサンプル画像（服部恵美さん）
・木と電話と選挙(causalTree)（安井翔太さん）
・スライド未公開、dplyrの話（@tomomoto）
・てかLINEやってる？(仮)（@wonder_zone）
・心理学における「再現性」の問題とBayes Factor（@NSushi）

・石田基広さんのキーノート

スライド未公開です。

・Linux使い
・ヘブライ語の意味構造を代数学でやっていた
・S/R言語の生みの親はJohn Chambers
　以下の二つは最近書かれた本だそうです。
　Software for Data Analysis: Programming with R (Statistics and Computing)
　Extending R (Chapman & Hall/CRC The R Series)
・S→S-plus→Rの順番で発展
・purrrを最近使い始めたそうです。
・XLConnectパッケージを使って、大学教員の採点活動を効率化しているそうです。

・ホクソエムとは何だったのか（ホクソエムさん）

匿名技術者集団ホクソエムの2016年の成果
・densratio（ densratio: Density Ratio Estimation ）
・githubinstall
（ githubinstall: A Helpful Way to Install R Packages Hosted on GitHub ）
・healthplanet（ Wrapper package for healthplanet api ）
・RODBCDBI
（ RODBCDBI: Provides Access to Databases Through the ODBC Interface ）
・jpmesh（ jpmesh: Utilities for Japanese Mesh Code ）

起業されたとのことです。懸命に頑張って下さい！
株式会社ホクソエム

awesomeな人材が必要とのことで、awesomeな方はアプライしてみてはいかがでしょうか。

・Rと探索的データ分析で、国連での日本の立ち位置を可視化する（安田洋介さん）

スライド未公開です。
国連のデータを使って、Exploratoryを用いた探索的データ分析の実演をされていました。

・マウス操作でかんたん予測分析（鈴木了太さん）

R AnalyticFlow
Rで実践！データサイエンス～初めの一歩から高度な応用まで～

・高速・省メモリにlibsvm形式でダンプする方法を研究してみた（@hskksk）

高速・省メモリにlibsvm形式でダンプする方法を研究してみた from Keisuke Hosaka

・Rでてんしょくかつどう（@Med_KU）

Rmd でreveal.js のhtml スライドプレゼンテーション

・RStudio vs Emacs（@y__mattu）

RStudio vs Emacs Japan.R 2016

・randomforestで高次元の変数重要度見る（@siero5335）

Randomforestで高次元の変数重要度を見る #japanr LT from Akifumi Eguchi

・Rで本を作りたい（前田和寛さん）

Rで本を作りたい

・28歳でプログラミングを始めた話（市川太祐さん）

・医療関連のアプリ開発でデータサイエンスを駆使しようとしているそうです。

スライド未公開です。
スライドがシェアされ次第載せます。

・LDA-Visパッケージのご紹介（@doradora09）

20161127 doradora09 japanr2016_lt from Nobuaki Oshiro

・【e2d3R】E2D3からDot-Bar-Chartのご紹介（楠本一哲さん）

スライドは未公開です。
E2D3をRで表示する試みのようです。
Experiments with e2d3 in R

・このIRのグラフがすごい！上場企業2016（@ito_yan）

スライド未公開です。後日シェアしていただけるようです。

今年のプレゼンはネタを集めるにしても、これが限界でした…楽しんで頂けたでしょうか？スライドは後日アップします。 #JapanR

— いとーちゃん (@ito_yan) 2016年11月27日

・Rでカルマンフィルタをしたい（@tetsuroito）

・PPAP（仮）（@yutannihilation）

Tidyverseとは from yutannihilation

・スライド未公開、ユーザーの状態遷移に関する分析のお話（@sanoche16）

スライドがシェアされ次第載せます。

・私とR（高栁慎一さん）

RjpWiki
統計・データ解析
 統計解析フリーソフト R の備忘録頁 ver.3.1
seekR(R限定の検索エンジン)
からだにいいもの
アブラタニブログってなんでしょう。油谷さんのブログ？

・めくってもめくってもサンプル画像（服部恵美さん）

Rのサンプルコードはあるけれども、どんな図ができるのかはわからない。そこで、サンプルコードとグラフを大量にまとめているサイトを作ったそうです。検索性は未知数ですが、暇なときに眺めておきたいですね。
R Graphical Manual

・木と電話と選挙(causalTree)（安井翔太さん）

木と電話と選挙(causalTree) from Shota Yasui

・スライド未公開、dplyrの話（@tomomoto）

スライドがシェアされ次第載せます。

・てかLINEやってる？(仮)（@wonder_zone）

てかLINEやってる？ (Japan.R 2016 LT) #JapanR from cancolle

・心理学における「再現性」の問題とBayes Factor（@NSushi）

スライドは後日公開とのことです。

今日はjapan.Rお疲れ様でしたー！楽しかったですありがとうございましたー！
ネット環境がないので心理学の再現性およびベイズファクターおよびにこにースライドは明日アップしますー。

— numba shushi (@NSushi) 2016年11月27日

『マーケティング・サイエンス入門』に出てくる手法をRで実行してみる

友人に『マーケティング・サイエンス入門』がおすすめと言われて読んだんですが、やっぱり実行できないとモヤモヤしてしまいますよね。そこで、登場する手法に関連したRのコードやらを集めてみました。

・BASSモデル
・多次元尺度法
・因子分析
・ロジット&プロビット
・分散分析
・クラスター分析
・判別分析
・決定木
・コンジョイント分析
・RFM分析
・共分散構造分析

BASSモデル

市場全体の規模が動的にどのように変化するかを予測するために使われるモデル。
R を使ってバスモデルを当てはめてみた – 廿TT
こちらにRのコードや適用例がいくつか載っています。

早速、私も携帯電話の加入契約数の時系列データを用いて、コードを実行してみました。データは平成25年版の総務省の情報通信白書の表から得ました。（第2部　情報通信の現況・政策の動向）

> fit_logis$value
[1] 69833.02
> fit_bass$value
[1] 68588.92

> fit_logis$value

[1] 69833.02

> fit_bass$value

[1] 68588.92

当てはまりはわずかながら、BASSモデルの方が良いようです。

多次元尺度法

多次元尺度法で遊んでみる（オレ流 R入門）
こちらのブログで山手線の駅間の距離データの可視化がなされています。
各駅ごとの距離からなる行列さえ用意すれば、cmdscale()関数を実行することで可能なようです。

今回はContaminatedMixtパッケージに含まれているワインのデータセットを使って多次元尺度法を適用してみようと思います。

データはこんな感じです。

以下のコードで実行しました。

library(ContaminatedMixt)

data(wine)
wine.dist <- dist(wine[,-1])
wine.cmd  <- cmdscale(wine.dist)

plot(wine.cmd,type = "n")
wine.lab <- factor(wine$Type)
text(wine.cmd,labels = wine.lab,col = unclass(wine.lab))

library(ContaminatedMixt)

data(wine)

wine.dist <- dist(wine[,-1])

wine.cmd <- cmdscale(wine.dist)

plot(wine.cmd,type = "n")

wine.lab <- factor(wine$Type)

text(wine.cmd,labels = wine.lab,col = unclass(wine.lab))

Barbera（バルベーラ）・・・基本的にはタンニンをあまり含まず、酸味の強い色の濃い赤ワインで庶民的。
Barolo（バローロ）・・・アルコール度数が高く、非常に重厚な味わいのワインでワインの王様と呼ばれる。
Grignolino（グリニョリーノ）・・・僅かにタンニンを感じるサッパリとした辛口の赤ワインで庶民的。

庶民と王様のワインは成分においても違いがありそうですね。

因子分析

psychパッケージというものがあるようです。こちらのサイトを参考にして進めます。（スナック菓子の食感についてRで因子分析してみた）
今回は大好きなwiskyのデータセットを使ってみます。（ Classification of whiskies ）

library(psych)
library(GPArotation)

whiskies <- read.csv("whiskies.txt", row.names = 1, stringsAsFactors = FALSE)
rownames(whiskies) <- whiskies[,1]
whiskies <- whiskies[,-1]
whiskies <- whiskies[,-15]
whiskies <- whiskies[,-14]
whiskies <- whiskies[,-13]

MAPminres <- vss(whiskies, fm="minres")
print(MAPminres)

MAPml <- vss(whiskies, fm="ml")
print(MAPml)

parallel <- fa.parallel(whiskies)
print(parallel)

res <- fa(whiskies, nfactors=2, fm="minres", rotate="oblimin")
print(res, digits=3)

biplot(res, labels=rownames(whiskies))

library(psych)

library(GPArotation)

whiskies <- read.csv("whiskies.txt", row.names = 1, stringsAsFactors = FALSE)

rownames(whiskies) <- whiskies[,1]

whiskies <- whiskies[,-1]

whiskies <- whiskies[,-15]

whiskies <- whiskies[,-14]

whiskies <- whiskies[,-13]

MAPminres <- vss(whiskies, fm="minres")

print(MAPminres)

MAPml <- vss(whiskies, fm="ml")

print(MAPml)

parallel <- fa.parallel(whiskies)

print(parallel)

res <- fa(whiskies, nfactors=2, fm="minres", rotate="oblimin")

print(res, digits=3)

biplot(res, labels=rownames(whiskies))

グレンフィディックやカリラやタリスカーがイメージ通りにプロットされています。ラガブーリンやアードベッグがはみ出しているのが残念ですが。

ロジット・プロビット

これらの手法はビルトインの関数でできてしまいますが、せっかくウイスキーのデータがあるので、薬っぽさに繋がりそうな変数を見つけてみます。

library(dplyr)

whiskies2 <- whiskies %>% mutate(target=ifelse(whiskies$Medicinal>0,1,0))
whiskies2 <- whiskies2[,-which(colnames(whiskies2)=="Medicinal")]

#ロジット
model_logit <- glm(target ~ ., data=whiskies2, family=binomial(link="logit"))
summary(model_logit)

#プロビット
model_probit <-glm(target ~ ., data=whiskies2, family=binomial(link="probit"))
summary(model_probit)

library(dplyr)

whiskies2 <- whiskies %>% mutate(target=ifelse(whiskies$Medicinal>0,1,0))

whiskies2 <- whiskies2[,-which(colnames(whiskies2)=="Medicinal")]

#ロジット

model_logit <- glm(target ~ ., data=whiskies2, family=binomial(link="logit"))

summary(model_logit)

#プロビット

model_probit <-glm(target ~ ., data=whiskies2, family=binomial(link="probit"))

summary(model_probit)

推定結果はこちらです。スモーキーさが関係しているのは納得です。

> model_logit <- glm(target ~ ., data=whiskies2, family=binomial(link="logit"))
> summary(model_logit)

Call:
glm(formula = target ~ ., family = binomial(link = "logit"), 
    data = whiskies2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.0173  -0.6313  -0.3509   0.2206   2.4214  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  4.39412    2.47716   1.774   0.0761 .
Body        -0.06226    0.51985  -0.120   0.9047  
Sweetness   -0.59932    0.46458  -1.290   0.1970  
Smoky        1.35983    0.58327   2.331   0.0197 *
Tobacco      1.75440    1.10595   1.586   0.1127  
Honey       -0.39101    0.46604  -0.839   0.4015  
Spicy       -0.22916    0.44331  -0.517   0.6052  
Winey       -0.86015    0.49829  -1.726   0.0843 .
Nutty       -0.56436    0.42843  -1.317   0.1878  
Malty       -1.17076    0.66243  -1.767   0.0772 .
Fruity      -0.33978    0.43754  -0.777   0.4374  
Floral      -0.64322    0.47291  -1.360   0.1738  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 107.023  on 85  degrees of freedom
Residual deviance:  64.882  on 74  degrees of freedom
AIC: 88.882

Number of Fisher Scoring iterations: 6

> model_probit <-glm(target ~ ., data=whiskies2, family=binomial(link="probit"))
> summary(model_probit)

Call:
glm(formula = target ~ ., family = binomial(link = "probit"), 
    data = whiskies2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9365  -0.6399  -0.3491   0.2064   2.4188  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  2.53606    1.41370   1.794   0.0728 .
Body        -0.05535    0.28834  -0.192   0.8478  
Sweetness   -0.36200    0.26740  -1.354   0.1758  
Smoky        0.76890    0.32649   2.355   0.0185 *
Tobacco      1.01419    0.64201   1.580   0.1142  
Honey       -0.19843    0.26056  -0.762   0.4463  
Spicy       -0.13537    0.25214  -0.537   0.5913  
Winey       -0.50543    0.27073  -1.867   0.0619 .
Nutty       -0.30335    0.24152  -1.256   0.2091  
Malty       -0.65838    0.36494  -1.804   0.0712 .
Fruity      -0.20000    0.25160  -0.795   0.4267  
Floral      -0.35983    0.26557  -1.355   0.1754  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 107.02  on 85  degrees of freedom
Residual deviance:  65.06  on 74  degrees of freedom
AIC: 89.06

Number of Fisher Scoring iterations: 7

> model_logit <- glm(target ~ ., data=whiskies2, family=binomial(link="logit"))

> summary(model_logit)

Call:

glm(formula = target ~ ., family = binomial(link = "logit"),

data = whiskies2)

Deviance Residuals:

Min 1Q Median 3Q Max

-2.0173 -0.6313 -0.3509 0.2206 2.4214

Coefficients:

Estimate Std. Error z value Pr(>|z|)

(Intercept) 4.39412 2.47716 1.774 0.0761 .

Body -0.06226 0.51985 -0.120 0.9047

Sweetness -0.59932 0.46458 -1.290 0.1970

Smoky 1.35983 0.58327 2.331 0.0197 *

Tobacco 1.75440 1.10595 1.586 0.1127

Honey -0.39101 0.46604 -0.839 0.4015

Spicy -0.22916 0.44331 -0.517 0.6052

Winey -0.86015 0.49829 -1.726 0.0843 .

Nutty -0.56436 0.42843 -1.317 0.1878

Malty -1.17076 0.66243 -1.767 0.0772 .

Fruity -0.33978 0.43754 -0.777 0.4374

Floral -0.64322 0.47291 -1.360 0.1738

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 107.023 on 85 degrees of freedom

Residual deviance: 64.882 on 74 degrees of freedom

AIC: 88.882

Number of Fisher Scoring iterations: 6

> model_probit <-glm(target ~ ., data=whiskies2, family=binomial(link="probit"))

> summary(model_probit)

Call:

glm(formula = target ~ ., family = binomial(link = "probit"),

data = whiskies2)

Deviance Residuals:

Min 1Q Median 3Q Max

-1.9365 -0.6399 -0.3491 0.2064 2.4188

Coefficients:

Estimate Std. Error z value Pr(>|z|)

(Intercept) 2.53606 1.41370 1.794 0.0728 .

Body -0.05535 0.28834 -0.192 0.8478

Sweetness -0.36200 0.26740 -1.354 0.1758

Smoky 0.76890 0.32649 2.355 0.0185 *

Tobacco 1.01419 0.64201 1.580 0.1142

Honey -0.19843 0.26056 -0.762 0.4463

Spicy -0.13537 0.25214 -0.537 0.5913

Winey -0.50543 0.27073 -1.867 0.0619 .

Nutty -0.30335 0.24152 -1.256 0.2091

Malty -0.65838 0.36494 -1.804 0.0712 .

Fruity -0.20000 0.25160 -0.795 0.4267

Floral -0.35983 0.26557 -1.355 0.1754

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 107.02 on 85 degrees of freedom

Residual deviance: 65.06 on 74 degrees of freedom

AIC: 89.06

Number of Fisher Scoring iterations: 7

ちなみに、多項ロジットに関しては、mlogitパッケージを使えばできるようです。（多項ロジット（Multinomial Logit）, R – mlogit 使用メモ）大学院時代に多項ロジットはSTATAでよく使っていましたが、Rだとこのパッケージなんですかね。推定した係数の値の解釈が若干複雑だったりします。

分散分析

分散分析もビルトインの関数で実行することができます。今回はワインのデータを用いて、銘柄から30個ランダムサンプリングをした上で、アルコールに関して群間の母平均値が同じかどうかを確かめてみます。コードはこちらを参考にしました。（ R による分散分析（一元配置））

A <- wine %>% filter(Type=="Barbera")
sampleNum_A <- sample(nrow(A),30)
A <- A[sampleNum_A,]

B <- wine %>% filter(Type=="Barolo")
sampleNum_B <- sample(nrow(B),30)
B <- B[sampleNum_B,]

C <- wine %>% filter(Type=="Grignolino")
sampleNum_C <- sample(nrow(C),30)
C <- C[sampleNum_C,]

d <- data.frame(A = A$Alcohol, B = B$Alcohol, C = C$Alcohol )
library(reshape)
x <- melt(d, variable_name = "group")
res <- aov(value ~ factor(group), data = x)
summary(res)

A <- wine %>% filter(Type=="Barbera")

sampleNum_A <- sample(nrow(A),30)

A <- A[sampleNum_A,]

B <- wine %>% filter(Type=="Barolo")

sampleNum_B <- sample(nrow(B),30)

B <- B[sampleNum_B,]

C <- wine %>% filter(Type=="Grignolino")

sampleNum_C <- sample(nrow(C),30)

C <- C[sampleNum_C,]

d <- data.frame(A = A$Alcohol, B = B$Alcohol, C = C$Alcohol )

library(reshape)

x <- melt(d, variable_name = "group")

res <- aov(value ~ factor(group), data = x)

summary(res)

推定結果はこちらです。アルコールに関しては、3群間において差があるようです。

> summary(res)
              Df Sum Sq Mean Sq F value Pr(>F)    
factor(group)  2  38.12  19.059   91.22 <2e-16 ***
Residuals     87  18.18   0.209                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> summary(res)

Df Sum Sq Mean Sq F value Pr(>F)

factor(group) 2 38.12 19.059 91.22 <2e-16 ***

Residuals 87 18.18 0.209

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

クラスター分析

クラスター分析もビルトインの関数で実行可能です。ここでは参考文献（ K-means Clustering 86 Single Malt Scotch Whiskies ）のウイスキーのサンプルで取り上げられたK-mean法をそのまま紹介します。

whiskies_k <- scale(whiskies)
ssPlot <- function(data, maxCluster = 9) {
  # Initialize within sum of squares
  SSw <- (nrow(data) - 1) * sum(apply(data, 2, var))
  SSw <- vector()
  for (i in 2:maxCluster) {
    SSw[i] <- sum(kmeans(data, centers = i)$withinss)
  }
  plot(1:maxCluster, SSw, type = "b", xlab = "Number of Clusters", ylab = "Within groups sum of squares")
}
ssPlot(whiskies_k)

fit <- kmeans(whiskies_k, 4) 
whiskies <- data.frame(whiskies, fit$cluster)
whiskies$fit.cluster <- as.factor(whiskies$fit.cluster)

subset(whiskies, fit.cluster == 2)

whiskies_k <- scale(whiskies)

ssPlot <- function(data, maxCluster = 9) {

# Initialize within sum of squares

SSw <- (nrow(data) - 1) * sum(apply(data, 2, var))

SSw <- vector()

for (i in 2:maxCluster) {

SSw[i] <- sum(kmeans(data, centers = i)$withinss)

}

plot(1:maxCluster, SSw, type = "b", xlab = "Number of Clusters", ylab = "Within groups sum of squares")

}

ssPlot(whiskies_k)

fit <- kmeans(whiskies_k, 4)

whiskies <- data.frame(whiskies, fit$cluster)

whiskies$fit.cluster <- as.factor(whiskies$fit.cluster)

subset(whiskies, fit.cluster == 2)

気になるクラスターの結果ですが、どうやらアイラ島系のウイスキーのクラスターを作れたようです。

> subset(whiskies, fit.cluster == 2)
          Body Sweetness Smoky Medicinal Tobacco Honey Spicy Winey Nutty Malty Fruity
Ardbeg       4         1     4         4       0     0     2     0     1     2      1
Caol Ila     3         1     4         2       1     0     2     0     2     1      1
Clynelish    3         2     3         3       1     0     2     0     1     1      2
Lagavulin    4         1     4         4       1     0     1     2     1     1      1
Laphroig     4         2     4         4       1     0     0     1     1     1      0
Talisker     4         2     3         3       0     1     3     0     1     2      2
          Floral fit.cluster
Ardbeg         0           2
Caol Ila       1           2
Clynelish      0           2
Lagavulin      0           2
Laphroig       0           2
Talisker       0           2

> subset(whiskies, fit.cluster == 2)

Body Sweetness Smoky Medicinal Tobacco Honey Spicy Winey Nutty Malty Fruity

Ardbeg 4 1 4 4 0 0 2 0 1 2 1

Caol Ila 3 1 4 2 1 0 2 0 2 1 1

Clynelish 3 2 3 3 1 0 2 0 1 1 2

Lagavulin 4 1 4 4 1 0 1 2 1 1 1

Laphroig 4 2 4 4 1 0 0 1 1 1 0

Talisker 4 2 3 3 0 1 3 0 1 2 2

Floral fit.cluster

Ardbeg 0 2

Caol Ila 1 2

Clynelish 0 2

Lagavulin 0 2

Laphroig 0 2

Talisker 0 2

判別分析

MASSパッケージで実行可能です。線形識別関数の実行例がこちらの参考文献に載っていたので、ワインのデータで試してみます。（【Rによるデータサイエンス】線形判別分析）

library(MASS)
wine2 <- data.frame(wine[,-1],type=wine.lab)
even.n <- 2*(1:(nrow(wine2)/2))-1
wine2.train <- wine2[even.n,]
wine2.test  <- wine2[-even.n,]

#train
Z.lda<-lda(type~.,data=wine2.train)
Z.lda$means
Z.lda$scaling
Z.lda$means%*%Z.lda$scaling
apply(Z.lda$means%*%Z.lda$scaling,2,mean)
table(wine2.train[,ncol(wine2)],predict(Z.lda)$class)
plot(Z.lda,dimen=1)
plot(Z.lda,dimen=2)

#test
testResult <- predict(Z.lda,wine2.test[,-ncol(wine2)])

table(wine2.test[,ncol(wine2)],testResult$class)
testResult$class

plot(testResult$x,type="n")
text(testResult$x,labels=wine2.test$type)

library(MASS)

wine2 <- data.frame(wine[,-1],type=wine.lab)

even.n <- 2*(1:(nrow(wine2)/2))-1

wine2.train <- wine2[even.n,]

wine2.test <- wine2[-even.n,]

#train

Z.lda<-lda(type~.,data=wine2.train)

Z.lda$means

Z.lda$scaling

Z.lda$means%*%Z.lda$scaling

apply(Z.lda$means%*%Z.lda$scaling,2,mean)

table(wine2.train[,ncol(wine2)],predict(Z.lda)$class)

plot(Z.lda,dimen=1)

plot(Z.lda,dimen=2)

#test

testResult <- predict(Z.lda,wine2.test[,-ncol(wine2)])

table(wine2.test[,ncol(wine2)],testResult$class)

testResult$class

plot(testResult$x,type="n")

text(testResult$x,labels=wine2.test$type)

判別関数得点

さすがパッケージ用のデータセットだけあって、綺麗に分類できたようです。誤分類は2件だけです。

> table(wine2.test[,ncol(wine2)],testResult$class)
            
             Barbera Barolo Grignolino
  Barbera         24      0          0
  Barolo           0     29          0
  Grignolino       0      2         34

> table(wine2.test[,ncol(wine2)],testResult$class)

Barbera Barolo Grignolino

Barbera 24 0 0

Barolo 0 29 0

Grignolino 0 2 34

決定木

決定木はrpartパッケージで実行します。ウイスキーのデータを使って、薬っぽさを決める条件を探してみます。コードはこちらを参考にしました。（ R言語で決定木分析）

library(rpart)
library(partykit)
library(epitools)

result <- rpart(target~., data = whiskies2)
plot(as.party(result))

library(rpart)

library(partykit)

library(epitools)

result <- rpart(target~., data = whiskies2)

plot(as.party(result))

コンジョイント分析

conjointパッケージなるものがあるようです。こちらの参考文献を元に紹介します。（ Rでコンジョイント分析）

まずは直交表を作ってみます。

library(conjoint)

experiment <- expand.grid(
  麺の量   =c("マシ","普通","少なめ"),
  ニンニク =c("マシ","普通","少なめ"),
  野菜     =c("マシ","普通","少なめ"),
  あぶら   =c("マシ","普通","少なめ")
)

design.ort <- caFactorialDesign(data = experiment,type = "orthogonal")
design.ort

caEncodedDesign(design.ort)
cor(caEncodedDesign(design.ort))

library(conjoint)

experiment <- expand.grid(

麺の量 =c("マシ","普通","少なめ"),

ニンニク =c("マシ","普通","少なめ"),

野菜 =c("マシ","普通","少なめ"),

あぶら =c("マシ","普通","少なめ")

)

design.ort <- caFactorialDesign(data = experiment,type = "orthogonal")

design.ort

caEncodedDesign(design.ort)

cor(caEncodedDesign(design.ort))

> design.ort
   麺の量 ニンニク   野菜 あぶら
5    普通     普通   マシ   マシ
10   マシ     マシ   普通   マシ
27 少なめ   少なめ 少なめ   マシ
34   マシ   少なめ   マシ   普通
42 少なめ     普通   普通   普通
47   普通     マシ 少なめ   普通
57 少なめ     マシ   マシ 少なめ
71   普通   少なめ   普通 少なめ
76   マシ     普通 少なめ 少なめ

> caEncodedDesign(design.ort)
   麺の量 ニンニク 野菜 あぶら
5       2        2    1      1
10      1        1    2      1
27      3        3    3      1
34      1        3    1      2
42      3        2    2      2
47      2        1    3      2
57      3        1    1      3
71      2        3    2      3
76      1        2    3      3
> cor(caEncodedDesign(design.ort))
         麺の量 ニンニク 野菜 あぶら
麺の量        1        0    0      0
ニンニク      0        1    0      0
野菜          0        0    1      0
あぶら        0        0    0      1

> design.ort

麺の量ニンニク野菜あぶら

5 普通普通マシマシ

10 マシマシ普通マシ

27 少なめ少なめ少なめマシ

34 マシ少なめマシ普通

42 少なめ普通普通普通

47 普通マシ少なめ普通

57 少なめマシマシ少なめ

71 普通少なめ普通少なめ

76 マシ普通少なめ少なめ

> caEncodedDesign(design.ort)

麺の量ニンニク野菜あぶら

5 2 2 1 1

10 1 1 2 1

27 3 3 3 1

34 1 3 1 2

42 3 2 2 2

47 2 1 3 2

57 3 1 1 3

71 2 3 2 3

76 1 2 3 3

> cor(caEncodedDesign(design.ort))

麺の量ニンニク野菜あぶら

麺の量 1 0 0 0

ニンニク 0 1 0 0

野菜 0 0 1 0

あぶら 0 0 0 1

残念ながら、面白そうなデータがないので、サンプルについているお茶のデータを使ってみます。

library(conjoint)

data(tea)

print(tprof)
print(tlevn)
print(tprefm)

Conjoint(tprefm,tprof,tlevn)

library(conjoint)

data(tea)

print(tprof)

print(tlevn)

print(tprefm)

Conjoint(tprefm,tprof,tlevn)

価格が効用に与える影響の可視化です。

RFM分析

ほくそ笑むの親分がeasyRFMパッケージを作っていたようです。（ RFM 分析を簡単に実行できる R パッケージ easyRFM を作った）都合良く取引データがなかったので、kaggleの掲示板で落ちていたデータを使いました。（ Sample of transaction data ）

library(easyRFM)

transaction <- read.csv("transactions-sample.csv")
transaction <- transaction %>% mutate(payment=purchasequantity*purchaseamount)
transaction <- transaction[,which(colnames(transaction) %in% c("id","payment","date"))]
transaction$date <- as.character(transaction$date)

result <- rfm_auto(transaction, breaks=3)
head(result$rfm)
result$classes

library(easyRFM)

transaction <- read.csv("transactions-sample.csv")

transaction <- transaction %>% mutate(payment=purchasequantity*purchaseamount)

transaction <- transaction[,which(colnames(transaction) %in% c("id","payment","date"))]

transaction$date <- as.character(transaction$date)

result <- rfm_auto(transaction, breaks=3)

head(result$rfm)

result$classes

結果は以下のとおりです。

> head(result$rfm)
        id    Recency Frequency  Monetary RecencyClass FrequencyClass MonetaryClass
1    86246 2013-04-23     12609  89333.08            1              3             3
2    86252 2013-03-26     12087 104436.35            1              3             3
3 12262064 2013-06-25      1101   5116.08            3              1             1
4 12277270 2013-06-19      1361  10820.49            3              2             2
5 12332190 2013-06-13       684   8709.90            2              1             1
> result$classes
$recency_class
[1] "2013-03-26 00:00:00 to 2013-05-11" "2013-05-11 00:00:01 to 2013-06-18"
[3] "2013-06-18 00:00:01 to 2013-06-26"

$recency_class_days
[1] "92 to 46" "45 to 8"  "7 to 0"  

$frequency_class
[1] "681 to 1200"   "1201 to 8600"  "8601 to 13000"

$monetary_class
[1] "5100 to 9500"    "9501 to 64000"   "64001 to 110000"

> head(result$rfm)

id Recency Frequency Monetary RecencyClass FrequencyClass MonetaryClass

1 86246 2013-04-23 12609 89333.08 1 3 3

2 86252 2013-03-26 12087 104436.35 1 3 3

3 12262064 2013-06-25 1101 5116.08 3 1 1

4 12277270 2013-06-19 1361 10820.49 3 2 2

5 12332190 2013-06-13 684 8709.90 2 1 1

> result$classes

$recency_class

[1] "2013-03-26 00:00:00 to 2013-05-11" "2013-05-11 00:00:01 to 2013-06-18"

[3] "2013-06-18 00:00:01 to 2013-06-26"

$recency_class_days

[1] "92 to 46" "45 to 8" "7 to 0"

$frequency_class

[1] "681 to 1200" "1201 to 8600" "8601 to 13000"

$monetary_class

[1] "5100 to 9500" "9501 to 64000" "64001 to 110000"

共分散構造分析

semパッケージで実行可能です。こちらの参考文献のデータを用います。（ Rによるパス解析）データはこちらにあります。（練習用データ）

library(sem)

data01	<-	read.csv(file="regpath.csv",header=T,
                   fileEncoding="Shift-JIS")

model01	<-	specifyModel()
意欲	->	学力,	b1,	NA
意欲	->	適応,	b2,	NA
適応	->	学力,	b3, NA
学力	->	自尊,	b4,	NA
学力	<->	学力,	v1,	NA
自尊	<->	自尊,	v2,	NA
適応	<->	適応,	v3, NA
意欲	<->	意欲,	v4,	NA

cov01	<-	cov(data01[,2:5])
cov01

fit01	<-	sem(model=model01,	S=cov01,	N=nrow(data01))

summary(fit01,rsquare=T,fit.indices=c("GFI","AGFI","SRMR","RMSEA","AIC","BIC"))
stdCoef(fit01)

library(sem)

data01 <- read.csv(file="regpath.csv",header=T,

fileEncoding="Shift-JIS")

model01 <- specifyModel()

意欲 -> 学力, b1, NA

意欲 -> 適応, b2, NA

適応 -> 学力, b3, NA

学力 -> 自尊, b4, NA

学力 <-> 学力, v1, NA

自尊 <-> 自尊, v2, NA

適応 <-> 適応, v3, NA

意欲 <-> 意欲, v4, NA

cov01 <- cov(data01[,2:5])

cov01

fit01 <- sem(model=model01, S=cov01, N=nrow(data01))

summary(fit01,rsquare=T,fit.indices=c("GFI","AGFI","SRMR","RMSEA","AIC","BIC"))

stdCoef(fit01)

Rを使った分析（SEM）
こちらの方がパスの図も出力できるので、良いかもしれません。

参考文献

マーケティング・サイエンス入門 — 市場対応の科学的マネジメント新版 (有斐閣アルマ)

Cloud Vision APIをRで動かしてみる

Google Cloud Vision APIを使ってAKBのデータを把握する
こちらの所沢義男さん（偽名）のブログにカジュアルにVison APIを使う方法が載せられていたので、早速私も触ってみようと思います。（Vision APIのアカウント作成には住所や氏名や電話番号、クレカ情報が求められます。）

ラーメン二郎のメニュー表

紹介されていたコードを実行してみます。

CROWD_VISION_KEY <- "YOUR KEY"

library("httr")

f <- "http://livedoor.blogimg.jp/musashinonikki/imgs/9/b/9b62aa78.jpg"
tmp <- tempfile()
download.file(f, tmp)
img <- readBin(tmp, "raw", file.info(tmp)[1, "size"])
u <- paste0("https://vision.googleapis.com/v1/images:annotate?key=", CROWD_VISION_KEY)
body <- list(requests = list(image=list(content=img), 
                             features=list(type="TEXT_DETECTION")
                             )
             )
res <- POST(url=u,
            encode="json",
            body=body,
            content_type_json()
)
res <- content(res)
cat(res$responses[[1]]$textAnnotations[[1]]$description[1])

CROWD_VISION_KEY <- "YOUR KEY"

library("httr")

f <- "http://livedoor.blogimg.jp/musashinonikki/imgs/9/b/9b62aa78.jpg"

tmp <- tempfile()

download.file(f, tmp)

img <- readBin(tmp, "raw", file.info(tmp)[1, "size"])

u <- paste0("https://vision.googleapis.com/v1/images:annotate?key=", CROWD_VISION_KEY)

body <- list(requests = list(image=list(content=img),

features=list(type="TEXT_DETECTION")

)

res <- POST(url=u,

encode="json",

body=body,

content_type_json()

)

res <- content(res)

cat(res$responses[[1]]$textAnnotations[[1]]$description[1])

メニュー
券売機番号
券売機番号
11小ラーメン
730円14大ラーメン
830円
12小ラーメン豚入り880円15大ラーメン誘入り 980円
17小つけ麺
880円18大つけ麺
980円
1 19小つけ麺豚入り 1030円20大つけ麺豚入り 1130円
13ミニラーメン 680円
21ビール
500円24黒ウーロン茶
200円
22うずら卵(5個) 100円23かつお魚粉
100円
02刻み生姜
50円
初めてのお客様は
かミニラーメンがお勧めです。
0000
3883
! 1
11ナナ
1つ1
売4580
4-3
券1112
22
00000000
38838005
7880651
睡! | |
| |
券11
1191220

メニュー

券売機番号

11小ラーメン

730円14大ラーメン

830円

12小ラーメン豚入り880円15大ラーメン誘入り 980円

17小つけ麺

880円18大つけ麺

980円

1 19小つけ麺豚入り 1030円20大つけ麺豚入り 1130円

13ミニラーメン 680円

21ビール

500円24黒ウーロン茶

200円

22うずら卵(5個) 100円23かつお魚粉

100円

02刻み生姜

50円

初めてのお客様は

かミニラーメンがお勧めです。

0000

3883

! 1

11ナナ

1つ1

売4580

4-3

券1112

00000000

38838005

7880651

睡! | |

| |

券11

1191220

以下の表は実際の表記との比較ですが、光の反射で見にくくなっているところ以外は正しい表現を抽出できているようです。

大学におけるStanの講義資料などを探してみた

ブログよりは大学の講義の方が体系立てて学べるのではないか、効率的に勉強できるのではないかと思い、各大学が公開しているStanに関するサイトを調べてみました。調べ方は非常に簡単で、Google検索で「stan site:大学のドメイン名」でヒットした上位を基本的に見ています。

東京大学、京都大学、東北大学、大阪大学、慶應大学、早稲田大学、名古屋大学、同志社大学、etc…などを見ましたがなかなかweb上で公開されている資料が見つかりませんでした。非公開か大学ドメイン下での公開がされていないのだと思われます。そのため、検索対象を海外にまで広げてみました。（次回は教員のwebサイトを漁ってみようと思います。）

結論として、学ぶのにちょうどよいと思えるのは、神戸大学の資料とStanford大学の資料でした。今後はこの二つの資料も使って学習を進めていこうと思います。

神戸大学

政治学方法論 II (Research Methods in Political Science II)
ベイズ統計学の授業の内容が公開されています。教科書は「Bayesian Data Analysis, 3rd Edition. CRC Press.」です。ちなみに、こちらは無料のPDFが公開されています。（Bayesian Data Analysis, Third Edition(PDF)）

階層モデルとStan によるベイズ推定
階層ベイズモデルの説明とstanのコードが記されており、学習が捗ります。

講義のスライドはこちらにあります。（ yukiyanai/rm2-Bayes ）

東京工業大学

勉強用(STAN)
stanのコードが載っていました。
正規分布、線形回帰モデル、混合正規分布、ニューラルネットワーク、多種粒子Totally Asymmetric Simple Exclusion Process、混合正規分布でのクラスタリング、ロジスティック回帰などのコードがあるようです。

こちらはstanの説明用の資料です。（Stanによるハミルトニアンモンテカルロ法を用いたサンプリングについて）

Stanford University

Statistical Rethinking A Bayesian Course with Examples in R and Stan
youtubeで2015年の講義が見れるようです。( Statistical Rethinking Winter 2015 )
講義のスライドも公開されています。( Talks by Richard McElreath )
ゴーレムをモデルの引き合いに出して紹介しているのを見て、ユーモアセンスあるなぁと思いました。2016年版の資料も今後アップされると思うので、見逃せないですね。

Colombia University

Home page for the book, “Bayesian Data Analysis”
stanの開発チームの方がコロンビア大学の研究者なので、絶対にあるだろうと思いましたが、スライドとかは特にありませんでした。学生の講義ノートは筆記体で画像になっているので、あまり読むことはお勧めはしません。
私として嬉しいのは。Rstanで教科書のコードを実行するためのスクリプトがGitHubで公開されていることでした。（ avehtari/BDA_R_demos/demos_rstan/ ）

ExploratoryをさわってRでGoogleSearchConsoleのデータを集計・可視化

先日知り合った方から、Exploratoryの紹介をしていただき、Rをカジュアルに利用できる可視化ツールとのことで関心を持っていたのですが、さっそくこのブログのアクセス関連のデータを集計・可視化してみようと思います。

【目次】
・Exploratoryについて
・ダウンロードとインストール
・データの読み込み（Built in）と可視化
・データの読み込み（With R）と可視化
・dplyrの活用
・使ってみて思うところ

Exploratoryについて

ExploratoryはRを用いた、前処理・可視化・分析のためのデスクトップアプリケーションです。Amazon Redshift、Google BigQuery、PostgreSQL、MySQL、MongoDB、Presto、Google Analytics、Google Spreadsheet、Twitter、Web Page Scraping、CSV、Excel、JSON、Statistical files、R scriptなどの豊富なデータソースを利用できるようです。特に、Rのscriptでデータを取得できるのは魅力的だなぁと思いました。加えて、前処理に適したRの関数（dplyrなど）をカジュアルに扱えるようなコマンドラインもありますので、前処理しながら可視化するという作業がしやすいと思われます。分析に関しては、クラスター分析、分類、回帰、アソシエーションルール、相関、類似度、自然言語処理、文書感情分析などができるようです。まぁ、やりたい分析に関するパッケージは新しくRにインストールしてしまえばいいので、R使いにとって自由度の高い環境と言えます。

ダウンロードとインストール

こちらのExploratory公式サイト（ https://exploratory.io/download ）でemailアドレスなどのアカウント情報を入力してダウンロードできます。あとはインストールしてしまえば良いです。Rに関してもExploratoryに対応したバージョンをインストールしてくれます。後は、アプリケーションを開いて、Create Newのボタンを押して分析スタートです。

データの読み込み（Built in）と可視化

私はGoogle Analytics(GA)のデータをよく使うので、さっそくExploratoryで読み込んでみました。

自分のGAのプロパティやビューを選んで、任意の期間を選び、好きなDimensionやMetricsを選択して抽出されたデータを確認してOKすれば読み込みは完了です。非常に簡単でしたし、データを確認しながらDimensionやMetricsを選べるのは試行錯誤しやすく便利です。

ピボットテーブルが非常に使いやすく、デバイス別・ページ別の直帰率の集計も楽にできました。

データの読み込み（With R）と可視化

私が魅力に感じたのはデータ読み込みの際に、Rのscriptをそのまま使えることです。RMeCabを使えばテキストデータを形態素解析して可視化することもできますし、RFinanceなどで取得した株価のデータとかも使えるわけです。これまでに使ってたコードをそのまま使えるのは嬉しい限りです。今回は恥ずかしながら私のブログのGoogle検索結果のクエリをGoogle Search Console APIを使って取得し、可視化にチャレンジしてみます。

Google Search Consoleは以下のコマンドでAPIからデータ取得が可能です。

library("googleAuthR")
library("searchConsoleR")

gar_auth()
sc_websites <- list_websites()

own_site <- "http://kamonohashiperry.com/"
time_period_start <- "2016-07-01" ###任意の開始日を記入
time_period_end   <- "2016-09-20" ###任意の終了日を記入
filter_url <- c("page~~/")

dataset <- search_analytics(siteURL = sc_websites[match(own_site,sc_websites[,1]),1],
                 startDate = as.character(as.Date(time_period_start) ),
                 endDate = as.character(as.Date(time_period_end) ),
                 dimensions = c("query","page","device","date"),
                 dimensionFilterExp = filter_url,
                 searchType = "web", rowLimit = 5000)

library("googleAuthR")

library("searchConsoleR")

gar_auth()

sc_websites <- list_websites()

own_site <- "http://kamonohashiperry.com/"

time_period_start <- "2016-07-01" ###任意の開始日を記入

time_period_end <- "2016-09-20" ###任意の終了日を記入

filter_url <- c("page~~/")

dataset <- search_analytics(siteURL = sc_websites[match(own_site,sc_websites[,1]),1],

startDate = as.character(as.Date(time_period_start) ),

endDate = as.character(as.Date(time_period_end) ),

dimensions = c("query","page","device","date"),

dimensionFilterExp = filter_url,

searchType = "web", rowLimit = 5000)

Google検索クエリの検索結果上での私のブログのデバイス別、順位・CTR・インプレッション数・クリック数を3D散布図に描いてみました。

こんな感じで結果のグラフなどをシェアできるようです。触れば動きます。

時系列データに関してもサクサク作れたので、私のブログのデバイス別の平均順位を時系列で描写してみました。

カーソルを合わせたらデータの詳細が表示されます。

dplyrの活用

細かい前処理がしたい場合は、画面上にあるcommandのウインドウでdplyrのコマンドを叩けば痒いところに手が届かせることも可能です。

使ってみて思うところ

RやPythonを使っている我々にとって、Tableauを使うモチベーションが非常に低かったんですが、Rとつないで好き放題できる環境というのは非常に魅力的だと思いました。私は個人的にWebスクレイピングを結構しているので、テキストデータの可視化などにもチャレンジしてみたいです。数をこなして可視化のセンスを磨きたい。

scikit-learnのモジュールのGitHubでの利用頻度を調べてみた

『Python機械学習プログラミング』を読んで、scikit-learnのモジュールは充実しているなぁと感じたんですが、実際にWebサイトでUser Guide( http://scikit-learn.org/stable/user_guide.html )を見た所、この本に載り切らないような数多くの機械学習手法に応じたモジュールが用意されていました。そこで、世のデータサイエンティストはどのモジュールを良く使っているのだろうと気になったので、GitHubのSearchでヒットしたCodeの数を各モジュール単位で集計してみました。検索クエリは「scikit-learn + モジュール名」なので、正確なものではないのですが、相対的な利用頻度を見るぶんには使えるのではないかと思われます。

データ集計方法

・User Guideに登場するscikit-learnのモジュール名を集めています。
・教師付き学習か教師無し学習かどうかの判断は、User Guideで紹介されているモジュールかどうかで判断しています。
・GitHubのSearchで「scikit-learn + モジュール名」でヒットした件数をそのまま使っています。(2016年9月22日時点)

可視化コード

Jupyterで実行しています。

%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import pylab

module_df = pd.read_csv("scikit_learn_modulelist.csv")

pylab.figure(figsize=(5, 25))
sns.set_style("whitegrid")
ax = sns.barplot(x="Hitscode", y="Module", data=module_df)

%matplotlib inline

import matplotlib.pyplot as plt

import seaborn as sns

import pandas as pd

import pylab

module_df = pd.read_csv("scikit_learn_modulelist.csv")

pylab.figure(figsize=(5, 25))

sns.set_style("whitegrid")

ax = sns.barplot(x="Hitscode", y="Module", data=module_df)

教師付き学習編

圧倒的に多いのがSVM(Support Vector Machine)を扱っているSVCモジュールで、続いて定番のロジスティック回帰やRandom Forestが使われているようです。統計解析ではメジャーなはずの線形回帰が5位なのは、初歩的なのであまりコードがアップされていないのかもしれません。GBDTのモジュールももう少し上位にくるかと思ったんですが15位でした。DMLCのXGBoostモジュールを使っているのかもしれませんね。私も実際のところXGBoostを使ってますし。

教師無し学習編

主成分分析やK-mean法など因子分解などのモジュールが上位を占めています。LDA(Latent Dirichlet Allocation)がもっと上位に来ると思ったんですが、思えばGensimの方が充実しているなぁと思うので、このランキングは妥当なのかもしれません。私もLDAなどはGensimを使っていますし。

収集を終えて

・社内だとデータサイエンティストの方がいないので、scikit-learnのモジュールの利用状況を知れてマニアックな共感をすることができた。
・SVMは実践例が豊富そうなので分析事例を探せば良い発見があるかもしれない。
・scikit-learnのUser Guideは充実していたので、時間を作って向き合ってみたいと思った。

統計的因果推論に関するスライドとRのサンプルコード

導入・入門から実践までのスライドと実践するためのコードを載せています。因果推論を実務でカジュアルに使えるまで上達させたいものです。

以下、
・統計的因果推論に関するスライド
・資料
・用語整理
・統計的因果推論に関するRの実行例
・書籍情報
からなります。

統計的因果推論に関するスライド

統計的因果推論の学習を進める上でのヒントが記されていると思われるスライドです。カジュアルに回帰分析を行うことに関する注意も書かれています。

統計的因果推論勉強会　第1回 from Hikaru GOTO

因果推論の歴史や利用することのモチベーションが非常にわかりやすく書かれています。
多重共線性を引き合いに出していたりしているのも理解が捗ります。

相関と因果について考える：統計的因果推論、その(不)可能性の中心 from takehikoihayashi

以下の2つは星野先生の『調査観察データの統計科学』通称赤本の1~3章に関するスライドです。

星野「調査観察データの統計科学」第1＆2章 from Shuyo Nakatani

星野「調査観察データの統計科学」第3章 from Shuyo Nakatani

Rによる実践に関してのスライドです。コードが書かれているので実行してみると良いでしょう。

傾向スコア：その概念とRによる実装 from takehikoihayashi

傾向スコアの計算をRで実践しているスライドです。

Rで学ぶ傾向スコア解析入門 – 無作為割り当てが出来ない時の因果効果推定 – from Yohei Sato

資料

こちらの資料は定義が書かれていて助かります。
2016/01/23 (Sat) 第 3 回因果推論を学ぶ会

こちらは論文や教科書の紹介もしています。
プロペンシティスコア（Propensity score; PS）（１）－PSの正しい使い方

カーネギーメロン大学の機械学習における因果推論に関しての講義資料です。
Lecture Notes 17 Causal Inference

Googleのハル・バリアンがまとめたペーパーも良いです。
Causal Inference in Economics and Marketing

用語整理

ATE（Average Treatment Effect：平均処置効果（因果効果））

例で述べるとするならば、母集団全てのユーザーにバナーを見せた場合のブランド名検索による訪問数の平均と、全てのユーザーにバナーを見せなかった場合のブランド名検索による訪問数の平均の差として表されます。バナーを見てしまったユーザーにとって、バナーを見なかったら、という反実仮想なデータは当然存在しないので、そのままでは計算できません。ただし、RCT（Randomized Control Trial）、無作為化比較対照実験、の状況ではバイアスなく推定できるとされています。RCTなケースは限られていると思いますが。

ATT（Average Treatment Effect on the Treated：処置群における平均処置効果）

バナーを見せたユーザーにおける、バナーを見せた場合と見せなかった場合の差の期待値。マーケティングにおける施策のROIを計算する際に使うことが望ましいとされています。ATEと同じくRCTにおいてバイアスなく推定できるとされています。

ATU（Average Treatment Effect on the Untreated：対照群における平均処置効果）

バナーを見せていないユーザーにおける、バナーを見せた場合と見せなかった場合の差の期待値。マーケティング施策を拡大させるか否かを判断する際に使うことができます。ATEと同じくRCTにおいてバイアスなく推定できるとされています。

強い意味での無視可能性

共変量に対し求める強い仮定のことで、「バナーを見たか見てないかのバイナリーな変数」や「ブランド名検索での訪問数」などに影響を与えるような共変量に対し、共変量自体で条件をつけて期待値をとると、「バナーを見たか見てないかのバイナリーな変数」と「潜在的なブランド名検索での訪問数」が独立するような特徴が求められています。「バナーを見たか見てないかのバイナリーな変数」が「過去のサイト訪問数（共変量）」や「特定ページへの接触（共変量）」で、配信対象を割り振られている場合は、そのバイナリーな変数は「潜在的なブランド名検索での訪問数」に影響を与えないとされています。

マッチング

バナーを見せられたユーザーの持つ、共変量（サイトへの訪問数や、見たページのカテゴリなど）の値と同じ（完全マッチング）、あるいは近い（距離を使ったマッチング）共変量を持っているが、バナーを見せられていない他のユーザーを「同じ人」と見なして、「バナーを見た・見てない」の与える「ブランド名検索での訪問数」への因果効果を推定します。

傾向スコア（Propensity score）

処置への割り当ての確率。つまり、上述の例でいうところの、バナーを見せられる確率。確率なので、当然0〜1の間の値をとります。推定には2項ロジットモデルが使われているようです。真の傾向スコアを推定できれば、ATE・ATT・ATUを計算することが可能になるそうです。この理屈はベイズの定理より導くことができるようです。詳しくは資料の”第 3 回因果推論を学ぶ会”を見てみてください。

統計的因果推論に関するRの実行例

“Rで学ぶ傾向スコア解析入門 – 無作為割り当てが出来ない時の因果効果推定”で紹介されていたコードを以下に掲載します。

#Matchingパッケージを呼び出す
library(Matching)

#サンプルデータを呼び出す
data(lalonde)

#データの確認
head(lalonde)

# age educ black hisp married nodegr re74 re75     re78 u74 u75 treat
# 1  37   11     1    0       1      1    0    0  9930.05   1   1     1
# 2  22    9     0    1       0      1    0    0  3595.89   1   1     1
# 3  30   12     1    0       0      0    0    0 24909.50   1   1     1
# 4  27   11     1    0       0      1    0    0  7506.15   1   1     1
# 5  33    8     1    0       0      1    0    0   289.79   1   1     1
# 6  22    9     1    0       0      1    0    0  4056.49   1   1     1

#ロジスティック回帰モデルを用いて傾向スコアを推定する。
#re78に対する効果を調べたいので、説明変数から除外しておく。
#treatはCMなどを見せたかどうかのバイナリーデータ。
logi <- glm(treat~., data=lalonde[, -9], family = binomial)

#因果効果を計算するために、マッチングを推定する。
#今回の効果を見るためのre78をYに指定して、CMなどを見せたかどうかのtreatをTrに指定する。
#共変量としてロジスティック回帰の予測値logi$fittedを用いている。
nsw1 <- Match(Y = lalonde$re78, Tr=lalonde$treat, X =logi$fitted)

summary(nsw1)

# Estimate...  2138.6 
# AI SE......  797.76 
# T-stat.....  2.6807 
# p.val......  0.0073468 
# 
# Original number of observations..............  445 
# Original number of treated obs...............  185 
# Matched number of observations...............  185 
# Matched number of observations  (unweighted).  322

#マッチングのペアを確認する。
lalonde2 <- lalonde
lalonde2$id <- 1:nrow(lalonde2)
lalonde2$score <- logi$fitted

pair.df <- cbind(lalonde2[nsw1$index.treated, c("id","score")],
                 lalonde2[nsw1$index.control, c("id", "score")])

names(pair.df) <- c("t.id", "t.score", "c.id", "c.score")

head(pair.df)

# t.id   t.score c.id   c.score
# 1      1 0.3927536  357 0.3935865
# 2      2 0.2271642  231 0.2242716
# 3      3 0.5313484  261 0.5313484
# 4      4 0.3285956  254 0.3285956
# 4.1    4 0.3285956  328 0.3286097
# 4.2    4 0.3285956  333 0.3286097
#確かに近そうだ。

#キャリパーマッチング（ペアが特定の距離以上になる時はマッチングしないマッチング）
nsw2 <- Match(Y =lalonde2$re78, Tr=lalonde2$treat, X =logi$fitted, caliper = T)
summary(nsw2)

# Estimate...  2138.6 
# AI SE......  797.76 
# T-stat.....  2.6807 
# p.val......  0.0073468 
# 
# Original number of observations..............  445 
# Original number of treated obs...............  185 
# Matched number of observations...............  185 
# Matched number of observations  (unweighted).  322 
# 
# Caliper (SDs)........................................   TRUE 
# Number of obs dropped by 'exact' or 'caliper'  0

#IPW(Inverse Probability Weight)の計算
ivec1 <- lalonde$treat
ivec2 <- rep(1, nrow(lalonde)) - ivec1

ivec <- cbind(ivec1, ivec2)
head(ivec)

# ivec1 ivec2
# [1,]     1     0
# [2,]     1     0
# [3,]     1     0
# [4,]     1     0
# [5,]     1     0
# [6,]     1     0

iestp1 <- (ivec1/logi$fitted)*(length(ivec1)/sum(ivec1))
iestp2 <- (ivec2/logi$fitted)*(length(ivec2)/sum(ivec2))

iestp <- iestp1 + iestp2
head(iestp)

# 1         2         3         4         5         6 
# 6.124465 10.588840  4.526984  7.320260  6.053291  6.745951 

ipwe <- lm(re78~ivec - 1, weights=iestp, data = lalonde)
summary(ipwe)

# Call:
#   lm(formula = re78 ~ ivec - 1, data = lalonde, weights = iestp)
# 
# Weighted Residuals:
#   Min     1Q Median     3Q    Max 
# -17241  -9985  -4048   6201 129799 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(&gt;|t|)    
# ivecivec1   6213.0      462.2   13.44   &lt;2e-16 ***
#   ivecivec2   4589.4      436.4   10.52   &lt;2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 15090 on 443 degrees of freedom
# Multiple R-squared:  0.3967,  Adjusted R-squared:  0.3939 
# F-statistic: 145.6 on 2 and 443 DF,  p-value: &lt; 2.2e-16

#因果効果
causal_effect <- ipwe$coefficients[1]- ipwe$coefficients[2]
causal_effect

# ivecivec1 
# 1623.581

#標準誤差
std_causal_effect <- sqrt(summary(ipwe)$coefficients[3]^2 + summary(ipwe)$coefficients[4]^2)
std_causal_effect

# [1] 635.7069

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

#Matchingパッケージを呼び出す

library(Matching)

#サンプルデータを呼び出す

data(lalonde)

#データの確認

head(lalonde)

# age educ black hisp married nodegr re74 re75 re78 u74 u75 treat

# 1 37 11 1 0 1 1 0 0 9930.05 1 1 1

# 2 22 9 0 1 0 1 0 0 3595.89 1 1 1

# 3 30 12 1 0 0 0 0 0 24909.50 1 1 1

# 4 27 11 1 0 0 1 0 0 7506.15 1 1 1

# 5 33 8 1 0 0 1 0 0 289.79 1 1 1

# 6 22 9 1 0 0 1 0 0 4056.49 1 1 1

#ロジスティック回帰モデルを用いて傾向スコアを推定する。

#re78に対する効果を調べたいので、説明変数から除外しておく。

#treatはCMなどを見せたかどうかのバイナリーデータ。

logi <- glm(treat~., data=lalonde[, -9], family = binomial)

#因果効果を計算するために、マッチングを推定する。

#今回の効果を見るためのre78をYに指定して、CMなどを見せたかどうかのtreatをTrに指定する。

#共変量としてロジスティック回帰の予測値logi$fittedを用いている。

nsw1 <- Match(Y = lalonde$re78, Tr=lalonde$treat, X =logi$fitted)

summary(nsw1)

# Estimate... 2138.6

# AI SE...... 797.76

# T-stat..... 2.6807

# p.val...... 0.0073468

# Original number of observations.............. 445

# Original number of treated obs............... 185

# Matched number of observations............... 185

# Matched number of observations (unweighted). 322

#マッチングのペアを確認する。

lalonde2 <- lalonde

lalonde2$id <- 1:nrow(lalonde2)

lalonde2$score <- logi$fitted

pair.df <- cbind(lalonde2[nsw1$index.treated, c("id","score")],

lalonde2[nsw1$index.control, c("id", "score")])

names(pair.df) <- c("t.id", "t.score", "c.id", "c.score")

head(pair.df)

# t.id t.score c.id c.score

# 1 1 0.3927536 357 0.3935865

# 2 2 0.2271642 231 0.2242716

# 3 3 0.5313484 261 0.5313484

# 4 4 0.3285956 254 0.3285956

# 4.1 4 0.3285956 328 0.3286097

# 4.2 4 0.3285956 333 0.3286097

#確かに近そうだ。

#キャリパーマッチング（ペアが特定の距離以上になる時はマッチングしないマッチング）

nsw2 <- Match(Y =lalonde2$re78, Tr=lalonde2$treat, X =logi$fitted, caliper = T)

summary(nsw2)

# Estimate... 2138.6

# AI SE...... 797.76

# T-stat..... 2.6807

# p.val...... 0.0073468

# Original number of observations.............. 445

# Original number of treated obs............... 185

# Matched number of observations............... 185

# Matched number of observations (unweighted). 322

# Caliper (SDs)........................................ TRUE

# Number of obs dropped by 'exact' or 'caliper' 0

#IPW(Inverse Probability Weight)の計算

ivec1 <- lalonde$treat

ivec2 <- rep(1, nrow(lalonde)) - ivec1

ivec <- cbind(ivec1, ivec2)

head(ivec)

# ivec1 ivec2

# [1,] 1 0

# [2,] 1 0

# [3,] 1 0

# [4,] 1 0

# [5,] 1 0

# [6,] 1 0

iestp1 <- (ivec1/logi$fitted)*(length(ivec1)/sum(ivec1))

iestp2 <- (ivec2/logi$fitted)*(length(ivec2)/sum(ivec2))

iestp <- iestp1 + iestp2

head(iestp)

# 1 2 3 4 5 6

# 6.124465 10.588840 4.526984 7.320260 6.053291 6.745951

ipwe <- lm(re78~ivec - 1, weights=iestp, data = lalonde)

summary(ipwe)

# Call:

# lm(formula = re78 ~ ivec - 1, data = lalonde, weights = iestp)

# Weighted Residuals:

# Min 1Q Median 3Q Max

# -17241 -9985 -4048 6201 129799

# Coefficients:

# Estimate Std. Error t value Pr(>|t|)

# ivecivec1 6213.0 462.2 13.44 <2e-16 ***

# ivecivec2 4589.4 436.4 10.52 <2e-16 ***

# ---

# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Residual standard error: 15090 on 443 degrees of freedom

# Multiple R-squared: 0.3967, Adjusted R-squared: 0.3939

# F-statistic: 145.6 on 2 and 443 DF, p-value: < 2.2e-16

#因果効果

causal_effect <- ipwe$coefficients[1]- ipwe$coefficients[2]

causal_effect

# ivecivec1

# 1623.581

#標準誤差

std_causal_effect <- sqrt(summary(ipwe)$coefficients[3]^2 + summary(ipwe)$coefficients[4]^2)

std_causal_effect

# [1] 635.7069

書籍情報

データ分析の力　因果関係に迫る思考法 (光文社新書)
「原因と結果」の経済学―――データから真実を見抜く思考法
 岩波データサイエンス Vol.3
調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)

XGBoostのパラメータチューニング実践 with Python

以前の投稿で紹介したXGBoostのパラメータチューニング方法ですが、実際のデータセットに対して実行するためのプログラムを実践してみようと思います。プログラム自体はAnalytics_Vidhya/Articles/Parameter_Tuning_XGBoost_with_Example/XGBoost models.ipynbに載っているのですが、データセットがついていません。そこで、前回の投稿(不均衡なデータの分類問題について with Python)で赤ワインのデータセットを手に入れているので、こちらのデータセットを用います。誤植なのかところどころ、うまく回らなかったところがあったので、手直しをしています。

以下の工程に従って進みます。結構長いですが、辛抱強く実践してみて下さい。
・ライブラリの読み込み
・データの読み込み
・前処理
・学習用データとテスト用データの作成
・XGBoostの予測結果をもとに、AUCの数値を返すための関数の定義
・モデルの実行
・チューニング

ライブラリの読み込み

import pandas as pd
import numpy as np
import xgboost as xgb
from xgboost.sklearn import XGBClassifier
from sklearn import cross_validation, metrics
from sklearn.grid_search import GridSearchCV

import matplotlib.pylab as plt
%matplotlib inline
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 12, 4

import pandas as pd

import numpy as np

import xgboost as xgb

from xgboost.sklearn import XGBClassifier

from sklearn import cross_validation, metrics

from sklearn.grid_search import GridSearchCV

import matplotlib.pylab as plt

%matplotlib inline

from matplotlib.pylab import rcParams

rcParams['figure.figsize'] = 12, 4

データの読み込み

#importing the red wine data
wine_df = pd.read_csv("https://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-red.csv", sep=";")

1 2	#importing the red wine data wine_df = pd.read_csv("https://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-red.csv", sep=";")

前処理

#ユニークIDを行ごとに割り当てる。
wine_df['ID'] = range(1, len(wine_df) + 1)

#ワインの質に関するデータを0-1データに置換する。qualityが6よりも小さかったら0、それ以外は1とする。
Y = wine_df.quality.values
wine_df.quality = np.asarray([1 if  i>=6 else 0 for i in Y])
wine_df.head(10)

#ユニークIDを行ごとに割り当てる。

wine_df['ID'] = range(1, len(wine_df) + 1)

#ワインの質に関するデータを0-1データに置換する。qualityが6よりも小さかったら0、それ以外は1とする。

Y = wine_df.quality.values

wine_df.quality = np.asarray([1 if i>=6 else 0 for i in Y])

wine_df.head(10)

学習用データとテスト用データの作成

#学習用データとテスト用データの作成
msk = np.random.rand(len(wine_df)) < 0.8 #乱数を発生させて0.8よりも小さいデータを選ぶ
train = wine_df[msk]
test = wine_df[~msk]

train.shape, test.shape
((1236, 13), (363, 13))

target='quality'
IDcol = 'ID'

#訓練データの目的変数の確認
train[target].value_counts()
1    659
0    577
Name: quality, dtype: int64

#学習用データとテスト用データの作成

msk = np.random.rand(len(wine_df)) < 0.8 #乱数を発生させて0.8よりも小さいデータを選ぶ

train = wine_df[msk]

test = wine_df[~msk]

train.shape, test.shape

((1236, 13), (363, 13))

target='quality'

IDcol = 'ID'

#訓練データの目的変数の確認

train[target].value_counts()

1 659

0 577

Name: quality, dtype: int64

XGBoostの予測結果をもとに、AUCの数値を返すための関数の定義

XGBoostの予測結果から、AUCの数値を返し、特徴量に応じた重要度を出力するためのプログラムです。

#テスト結果を格納するデータフレームの生成
test_results = pd.DataFrame(data=test.ID)

#関数の定義
def modelfit(alg, dtrain, dtest, predictors,useTrainCV=True, cv_folds=5, early_stopping_rounds=50):
    
    if useTrainCV:
        xgb_param = alg.get_xgb_params()
        xgtrain = xgb.DMatrix(dtrain[predictors].values, label=dtrain[target].values)
        xgtest = xgb.DMatrix(dtest[predictors].values)
        cvresult = xgb.cv(xgb_param, xgtrain, num_boost_round=alg.get_params()['n_estimators'], nfold=cv_folds,
            metrics=['auc'], early_stopping_rounds=early_stopping_rounds, show_progress=False)
        alg.set_params(n_estimators=cvresult.shape[0])
    
    #Fit the algorithm on the data
    alg.fit(dtrain[predictors], dtrain[target],eval_metric=['auc'])
        
    #Predict training set:
    dtrain_predictions = alg.predict(dtrain[predictors])
    dtrain_predprob = alg.predict_proba(dtrain[predictors])[:,1]
        
    #Print model report:
    print "\nModel Report"
    print "Accuracy : %.4g" % metrics.accuracy_score(dtrain[target].values, dtrain_predictions)
    print "AUC Score (Train): %f" % metrics.roc_auc_score(dtrain[target], dtrain_predprob)
    
    # Predict on testing data:
    dtest['predprob'] = alg.predict_proba(dtest[predictors])[:,1]
    #results = test_results.merge(dtest[['ID','predprob']], on='ID')
    print 'AUC Score (Test): %f' % metrics.roc_auc_score(dtest[target], dtest['predprob'])
                
    feat_imp = pd.Series(alg.booster().get_fscore()).sort_values(ascending=False)
    feat_imp.plot(kind='bar', title='Feature Importances')
    plt.ylabel('Feature Importance Score')

#テスト結果を格納するデータフレームの生成

test_results = pd.DataFrame(data=test.ID)

#関数の定義

def modelfit(alg, dtrain, dtest, predictors,useTrainCV=True, cv_folds=5, early_stopping_rounds=50):

if useTrainCV:

xgb_param = alg.get_xgb_params()

xgtrain = xgb.DMatrix(dtrain[predictors].values, label=dtrain[target].values)

xgtest = xgb.DMatrix(dtest[predictors].values)

cvresult = xgb.cv(xgb_param, xgtrain, num_boost_round=alg.get_params()['n_estimators'], nfold=cv_folds,

metrics=['auc'], early_stopping_rounds=early_stopping_rounds, show_progress=False)

alg.set_params(n_estimators=cvresult.shape[0])

#Fit the algorithm on the data

alg.fit(dtrain[predictors], dtrain[target],eval_metric=['auc'])

#Predict training set:

dtrain_predictions = alg.predict(dtrain[predictors])

dtrain_predprob = alg.predict_proba(dtrain[predictors])[:,1]

#Print model report:

print "\nModel Report"

print "Accuracy : %.4g" % metrics.accuracy_score(dtrain[target].values, dtrain_predictions)

print "AUC Score (Train): %f" % metrics.roc_auc_score(dtrain[target], dtrain_predprob)

# Predict on testing data:

dtest['predprob'] = alg.predict_proba(dtest[predictors])[:,1]

#results = test_results.merge(dtest[['ID','predprob']], on='ID')

print 'AUC Score (Test): %f' % metrics.roc_auc_score(dtest[target], dtest['predprob'])

feat_imp = pd.Series(alg.booster().get_fscore()).sort_values(ascending=False)

feat_imp.plot(kind='bar', title='Feature Importances')

plt.ylabel('Feature Importance Score')

モデルの実行

predictors = [x for x in train.columns if x not in [target, IDcol]]
xgb1 = XGBClassifier(
        learning_rate =0.1,
        n_estimators=1000,
        max_depth=5,
        min_child_weight=1,
        gamma=0,
        subsample=0.8,
        colsample_bytree=0.8,
        objective= 'binary:logistic',
        nthread=4,
        scale_pos_weight=1,
        seed=27)

modelfit(xgb1, train, test, predictors)

Will train until cv error hasn't decreased in 50 rounds.
Stopping. Best iteration: 237

Model Report
Accuracy : 1
AUC Score (Train): 1.000000
AUC Score (Test): 0.875199

predictors = [x for x in train.columns if x not in [target, IDcol]]

xgb1 = XGBClassifier(

learning_rate =0.1,

n_estimators=1000,

max_depth=5,

min_child_weight=1,

gamma=0,

subsample=0.8,

colsample_bytree=0.8,

objective= 'binary:logistic',

nthread=4,

scale_pos_weight=1,

seed=27)

modelfit(xgb1, train, test, predictors)

Will train until cv error hasn't decreased in 50 rounds.

Stopping. Best iteration: 237

Model Report

Accuracy : 1

AUC Score (Train): 1.000000

AUC Score (Test): 0.875199

チューニング

max_depthとmin_child_weightの数値をチューニングするためのプログラムです。

#Grid seach on subsample and max_features
#Choose all predictors except target & IDcols
param_test1 = {
    'max_depth':range(3,10,2),
    'min_child_weight':range(1,6,2)
}
gsearch1 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=1000, max_depth=5,
                                        min_child_weight=1, gamma=0, subsample=0.8, colsample_bytree=0.8,
                                        objective= 'binary:logistic', nthread=4, scale_pos_weight=1, seed=27), 
                       param_grid = param_test1, scoring='roc_auc',n_jobs=4,iid=False, cv=5)
gsearch1.fit(train[predictors],train[target])

#Grid seach on subsample and max_features

#Choose all predictors except target & IDcols

param_test1 = {

'max_depth':range(3,10,2),

'min_child_weight':range(1,6,2)

}

gsearch1 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=1000, max_depth=5,

min_child_weight=1, gamma=0, subsample=0.8, colsample_bytree=0.8,

objective= 'binary:logistic', nthread=4, scale_pos_weight=1, seed=27),

param_grid = param_test1, scoring='roc_auc',n_jobs=4,iid=False, cv=5)

gsearch1.fit(train[predictors],train[target])

gsearch1.grid_scores_, gsearch1.best_params_, gsearch1.best_score_

([mean: 0.76728, std: 0.03045, params: {'max_depth': 3, 'min_child_weight': 1},
  mean: 0.76649, std: 0.03378, params: {'max_depth': 3, 'min_child_weight': 3},
  mean: 0.76540, std: 0.03620, params: {'max_depth': 3, 'min_child_weight': 5},
  mean: 0.76509, std: 0.03183, params: {'max_depth': 5, 'min_child_weight': 1},
  mean: 0.76430, std: 0.02988, params: {'max_depth': 5, 'min_child_weight': 3},
  mean: 0.76221, std: 0.03336, params: {'max_depth': 5, 'min_child_weight': 5},
  mean: 0.77162, std: 0.03335, params: {'max_depth': 7, 'min_child_weight': 1},
  mean: 0.76575, std: 0.03585, params: {'max_depth': 7, 'min_child_weight': 3},
  mean: 0.76277, std: 0.03511, params: {'max_depth': 7, 'min_child_weight': 5},
  mean: 0.77235, std: 0.03283, params: {'max_depth': 9, 'min_child_weight': 1},
  mean: 0.76452, std: 0.03414, params: {'max_depth': 9, 'min_child_weight': 3},
  mean: 0.76114, std: 0.03561, params: {'max_depth': 9, 'min_child_weight': 5}],
 {'max_depth': 9, 'min_child_weight': 1},
 0.77235073909956886)

gsearch1.grid_scores_, gsearch1.best_params_, gsearch1.best_score_

([mean: 0.76728, std: 0.03045, params: {'max_depth': 3, 'min_child_weight': 1},

mean: 0.76649, std: 0.03378, params: {'max_depth': 3, 'min_child_weight': 3},

mean: 0.76540, std: 0.03620, params: {'max_depth': 3, 'min_child_weight': 5},

mean: 0.76509, std: 0.03183, params: {'max_depth': 5, 'min_child_weight': 1},

mean: 0.76430, std: 0.02988, params: {'max_depth': 5, 'min_child_weight': 3},

mean: 0.76221, std: 0.03336, params: {'max_depth': 5, 'min_child_weight': 5},

mean: 0.77162, std: 0.03335, params: {'max_depth': 7, 'min_child_weight': 1},

mean: 0.76575, std: 0.03585, params: {'max_depth': 7, 'min_child_weight': 3},

mean: 0.76277, std: 0.03511, params: {'max_depth': 7, 'min_child_weight': 5},

mean: 0.77235, std: 0.03283, params: {'max_depth': 9, 'min_child_weight': 1},

mean: 0.76452, std: 0.03414, params: {'max_depth': 9, 'min_child_weight': 3},

mean: 0.76114, std: 0.03561, params: {'max_depth': 9, 'min_child_weight': 5}],

{'max_depth': 9, 'min_child_weight': 1},

0.77235073909956886)

より細かい数値で再度最適なパラメータを探します。

#Grid seach on subsample and max_features
#Choose all predictors except target & IDcols
param_test2 = {
    'max_depth':[4,5,6,7,8,9],
    'min_child_weight':[1,2,3,4,5,6]
}
gsearch2 = GridSearchCV(estimator = XGBClassifier( learning_rate=0.1, n_estimators=1000, max_depth=5,
                                        min_child_weight=1, gamma=0, subsample=0.8, colsample_bytree=0.8,
                                        objective= 'binary:logistic', nthread=4, scale_pos_weight=1,seed=27), 
                       param_grid = param_test2, scoring='roc_auc',n_jobs=4,iid=False, cv=5)
gsearch2.fit(train[predictors],train[target])

#Grid seach on subsample and max_features

#Choose all predictors except target & IDcols

param_test2 = {

'max_depth':[4,5,6,7,8,9],

'min_child_weight':[1,2,3,4,5,6]

}

gsearch2 = GridSearchCV(estimator = XGBClassifier( learning_rate=0.1, n_estimators=1000, max_depth=5,

min_child_weight=1, gamma=0, subsample=0.8, colsample_bytree=0.8,

objective= 'binary:logistic', nthread=4, scale_pos_weight=1,seed=27),

param_grid = param_test2, scoring='roc_auc',n_jobs=4,iid=False, cv=5)

gsearch2.fit(train[predictors],train[target])

gsearch2.grid_scores_, gsearch2.best_params_, gsearch2.best_score_

([mean: 0.76820, std: 0.03336, params: {'max_depth': 4, 'min_child_weight': 1},
  mean: 0.76431, std: 0.02792, params: {'max_depth': 4, 'min_child_weight': 2},
  mean: 0.76171, std: 0.03356, params: {'max_depth': 4, 'min_child_weight': 3},
  mean: 0.76257, std: 0.03277, params: {'max_depth': 4, 'min_child_weight': 4},
  mean: 0.76128, std: 0.03661, params: {'max_depth': 4, 'min_child_weight': 5},
  mean: 0.75902, std: 0.03280, params: {'max_depth': 4, 'min_child_weight': 6},
  mean: 0.76509, std: 0.03183, params: {'max_depth': 5, 'min_child_weight': 1},
  mean: 0.76426, std: 0.02974, params: {'max_depth': 5, 'min_child_weight': 2},
  mean: 0.76430, std: 0.02988, params: {'max_depth': 5, 'min_child_weight': 3},
  mean: 0.76262, std: 0.02992, params: {'max_depth': 5, 'min_child_weight': 4},
  mean: 0.76221, std: 0.03336, params: {'max_depth': 5, 'min_child_weight': 5},
  mean: 0.76655, std: 0.03397, params: {'max_depth': 5, 'min_child_weight': 6},
  mean: 0.77066, std: 0.02936, params: {'max_depth': 6, 'min_child_weight': 1},
  mean: 0.76422, std: 0.03038, params: {'max_depth': 6, 'min_child_weight': 2},
  mean: 0.76126, std: 0.03021, params: {'max_depth': 6, 'min_child_weight': 3},
  mean: 0.76334, std: 0.03176, params: {'max_depth': 6, 'min_child_weight': 4},
  mean: 0.76347, std: 0.03245, params: {'max_depth': 6, 'min_child_weight': 5},
  mean: 0.76437, std: 0.03546, params: {'max_depth': 6, 'min_child_weight': 6},
  mean: 0.77162, std: 0.03335, params: {'max_depth': 7, 'min_child_weight': 1},
  mean: 0.76140, std: 0.03245, params: {'max_depth': 7, 'min_child_weight': 2},
  mean: 0.76575, std: 0.03585, params: {'max_depth': 7, 'min_child_weight': 3},
  mean: 0.76345, std: 0.03518, params: {'max_depth': 7, 'min_child_weight': 4},
  mean: 0.76277, std: 0.03511, params: {'max_depth': 7, 'min_child_weight': 5},
  mean: 0.75858, std: 0.03375, params: {'max_depth': 7, 'min_child_weight': 6},
  mean: 0.77487, std: 0.03377, params: {'max_depth': 8, 'min_child_weight': 1},
  mean: 0.76740, std: 0.03472, params: {'max_depth': 8, 'min_child_weight': 2},
  mean: 0.76048, std: 0.03267, params: {'max_depth': 8, 'min_child_weight': 3},
  mean: 0.76288, std: 0.03773, params: {'max_depth': 8, 'min_child_weight': 4},
  mean: 0.76138, std: 0.03045, params: {'max_depth': 8, 'min_child_weight': 5},
  mean: 0.76233, std: 0.03652, params: {'max_depth': 8, 'min_child_weight': 6},
  mean: 0.77235, std: 0.03283, params: {'max_depth': 9, 'min_child_weight': 1},
  mean: 0.76929, std: 0.03267, params: {'max_depth': 9, 'min_child_weight': 2},
  mean: 0.76452, std: 0.03414, params: {'max_depth': 9, 'min_child_weight': 3},
  mean: 0.76152, std: 0.03731, params: {'max_depth': 9, 'min_child_weight': 4},
  mean: 0.76114, std: 0.03561, params: {'max_depth': 9, 'min_child_weight': 5},
  mean: 0.76551, std: 0.03394, params: {'max_depth': 9, 'min_child_weight': 6}],
 {'max_depth': 8, 'min_child_weight': 1},
 0.77486987248915451)

gsearch2.grid_scores_, gsearch2.best_params_, gsearch2.best_score_

([mean: 0.76820, std: 0.03336, params: {'max_depth': 4, 'min_child_weight': 1},

mean: 0.76431, std: 0.02792, params: {'max_depth': 4, 'min_child_weight': 2},

mean: 0.76171, std: 0.03356, params: {'max_depth': 4, 'min_child_weight': 3},

mean: 0.76257, std: 0.03277, params: {'max_depth': 4, 'min_child_weight': 4},

mean: 0.76128, std: 0.03661, params: {'max_depth': 4, 'min_child_weight': 5},

mean: 0.75902, std: 0.03280, params: {'max_depth': 4, 'min_child_weight': 6},

mean: 0.76509, std: 0.03183, params: {'max_depth': 5, 'min_child_weight': 1},

mean: 0.76426, std: 0.02974, params: {'max_depth': 5, 'min_child_weight': 2},

mean: 0.76430, std: 0.02988, params: {'max_depth': 5, 'min_child_weight': 3},

mean: 0.76262, std: 0.02992, params: {'max_depth': 5, 'min_child_weight': 4},

mean: 0.76221, std: 0.03336, params: {'max_depth': 5, 'min_child_weight': 5},

mean: 0.76655, std: 0.03397, params: {'max_depth': 5, 'min_child_weight': 6},

mean: 0.77066, std: 0.02936, params: {'max_depth': 6, 'min_child_weight': 1},

mean: 0.76422, std: 0.03038, params: {'max_depth': 6, 'min_child_weight': 2},

mean: 0.76126, std: 0.03021, params: {'max_depth': 6, 'min_child_weight': 3},

mean: 0.76334, std: 0.03176, params: {'max_depth': 6, 'min_child_weight': 4},

mean: 0.76347, std: 0.03245, params: {'max_depth': 6, 'min_child_weight': 5},

mean: 0.76437, std: 0.03546, params: {'max_depth': 6, 'min_child_weight': 6},

mean: 0.77162, std: 0.03335, params: {'max_depth': 7, 'min_child_weight': 1},

mean: 0.76140, std: 0.03245, params: {'max_depth': 7, 'min_child_weight': 2},

mean: 0.76575, std: 0.03585, params: {'max_depth': 7, 'min_child_weight': 3},

mean: 0.76345, std: 0.03518, params: {'max_depth': 7, 'min_child_weight': 4},

mean: 0.76277, std: 0.03511, params: {'max_depth': 7, 'min_child_weight': 5},

mean: 0.75858, std: 0.03375, params: {'max_depth': 7, 'min_child_weight': 6},

mean: 0.77487, std: 0.03377, params: {'max_depth': 8, 'min_child_weight': 1},

mean: 0.76740, std: 0.03472, params: {'max_depth': 8, 'min_child_weight': 2},

mean: 0.76048, std: 0.03267, params: {'max_depth': 8, 'min_child_weight': 3},

mean: 0.76288, std: 0.03773, params: {'max_depth': 8, 'min_child_weight': 4},

mean: 0.76138, std: 0.03045, params: {'max_depth': 8, 'min_child_weight': 5},

mean: 0.76233, std: 0.03652, params: {'max_depth': 8, 'min_child_weight': 6},

mean: 0.77235, std: 0.03283, params: {'max_depth': 9, 'min_child_weight': 1},

mean: 0.76929, std: 0.03267, params: {'max_depth': 9, 'min_child_weight': 2},

mean: 0.76452, std: 0.03414, params: {'max_depth': 9, 'min_child_weight': 3},

mean: 0.76152, std: 0.03731, params: {'max_depth': 9, 'min_child_weight': 4},

mean: 0.76114, std: 0.03561, params: {'max_depth': 9, 'min_child_weight': 5},

mean: 0.76551, std: 0.03394, params: {'max_depth': 9, 'min_child_weight': 6}],

{'max_depth': 8, 'min_child_weight': 1},

0.77486987248915451)

max_depthを8、min_child_weightを1として、他のパラメータチューニングに移ります。
続いて、gammaのチューニングを行います。

#Grid seach on subsample and max_features
#Choose all predictors except target & IDcols
param_test3 = {
    'gamma':[i/10.0 for i in range(0,5)]
}
gsearch3 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=1000, max_depth=8,
                                        min_child_weight=1, gamma=0, subsample=0.8, colsample_bytree=0.8,
                                        objective= 'binary:logistic', nthread=4, scale_pos_weight=1,seed=27), 
                       param_grid = param_test3, scoring='roc_auc',n_jobs=4,iid=False, cv=5)
gsearch3.fit(train[predictors],train[target])

#Grid seach on subsample and max_features

#Choose all predictors except target & IDcols

param_test3 = {

'gamma':[i/10.0 for i in range(0,5)]

}

gsearch3 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=1000, max_depth=8,

min_child_weight=1, gamma=0, subsample=0.8, colsample_bytree=0.8,

objective= 'binary:logistic', nthread=4, scale_pos_weight=1,seed=27),

param_grid = param_test3, scoring='roc_auc',n_jobs=4,iid=False, cv=5)

gsearch3.fit(train[predictors],train[target])

gsearch3.grid_scores_, gsearch3.best_params_, gsearch3.best_score_

([mean: 0.77487, std: 0.03377, params: {'gamma': 0.0},
  mean: 0.77689, std: 0.03298, params: {'gamma': 0.1},
  mean: 0.77735, std: 0.03117, params: {'gamma': 0.2},
  mean: 0.78163, std: 0.03076, params: {'gamma': 0.3},
  mean: 0.78790, std: 0.03328, params: {'gamma': 0.4}],
 {'gamma': 0.4},
 0.78789976715320331)

gsearch3.grid_scores_, gsearch3.best_params_, gsearch3.best_score_

([mean: 0.77487, std: 0.03377, params: {'gamma': 0.0},

mean: 0.77689, std: 0.03298, params: {'gamma': 0.1},

mean: 0.77735, std: 0.03117, params: {'gamma': 0.2},

mean: 0.78163, std: 0.03076, params: {'gamma': 0.3},

mean: 0.78790, std: 0.03328, params: {'gamma': 0.4}],

{'gamma': 0.4},

0.78789976715320331)

gammaを0.4と置きます。
ここで、いままでにチューニングしたパラメータを用いて再度推定を行います。先ほどの0.875よりも高くなっています。

predictors = [x for x in train.columns if x not in [target, IDcol]]
xgb2 = XGBClassifier(
        learning_rate =0.1,
        n_estimators=1000,
        max_depth=8,
        min_child_weight=1,
        gamma=0.4,
        subsample=0.8,
        colsample_bytree=0.8,
        objective= 'binary:logistic',
        nthread=4,
        scale_pos_weight=1,
        seed=27)
modelfit(xgb2, train, test, predictors)


Will train until cv error hasn't decreased in 50 rounds.
Stopping. Best iteration: 120

Model Report
Accuracy : 1
AUC Score (Train): 1.000000
AUC Score (Test): 0.884028

predictors = [x for x in train.columns if x not in [target, IDcol]]

xgb2 = XGBClassifier(

learning_rate =0.1,

n_estimators=1000,

max_depth=8,

min_child_weight=1,

gamma=0.4,

subsample=0.8,

colsample_bytree=0.8,

objective= 'binary:logistic',

nthread=4,

scale_pos_weight=1,

seed=27)

modelfit(xgb2, train, test, predictors)

Will train until cv error hasn't decreased in 50 rounds.

Stopping. Best iteration: 120

Model Report

Accuracy : 1

AUC Score (Train): 1.000000

AUC Score (Test): 0.884028

続いて、subsampleとcolsample_bytreeのチューニングを行います。

#Grid seach on subsample and max_features
#Choose all predictors except target & IDcols
param_test4 = {
    'subsample':[i/10.0 for i in range(6,10)],
    'colsample_bytree':[i/10.0 for i in range(6,10)]
}
gsearch4 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=177, max_depth=8,
                                        min_child_weight=1, gamma=0.4, subsample=0.8, colsample_bytree=0.8,
                                        objective= 'binary:logistic', nthread=4, scale_pos_weight=1,seed=27), 
                       param_grid = param_test4, scoring='roc_auc',n_jobs=4,iid=False, cv=5)
gsearch4.fit(train[predictors],train[target])

#Grid seach on subsample and max_features

#Choose all predictors except target & IDcols

param_test4 = {

'subsample':[i/10.0 for i in range(6,10)],

'colsample_bytree':[i/10.0 for i in range(6,10)]

}

gsearch4 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=177, max_depth=8,

min_child_weight=1, gamma=0.4, subsample=0.8, colsample_bytree=0.8,

objective= 'binary:logistic', nthread=4, scale_pos_weight=1,seed=27),

param_grid = param_test4, scoring='roc_auc',n_jobs=4,iid=False, cv=5)

gsearch4.fit(train[predictors],train[target])

gsearch4.grid_scores_, gsearch4.best_params_, gsearch4.best_score_

([mean: 0.78994, std: 0.02779, params: {'subsample': 0.6, 'colsample_bytree': 0.6},
  mean: 0.78900, std: 0.03519, params: {'subsample': 0.7, 'colsample_bytree': 0.6},
  mean: 0.78509, std: 0.03202, params: {'subsample': 0.8, 'colsample_bytree': 0.6},
  mean: 0.78706, std: 0.02848, params: {'subsample': 0.9, 'colsample_bytree': 0.6},
  mean: 0.78511, std: 0.03140, params: {'subsample': 0.6, 'colsample_bytree': 0.7},
  mean: 0.78343, std: 0.03336, params: {'subsample': 0.7, 'colsample_bytree': 0.7},
  mean: 0.78939, std: 0.03203, params: {'subsample': 0.8, 'colsample_bytree': 0.7},
  mean: 0.78646, std: 0.04090, params: {'subsample': 0.9, 'colsample_bytree': 0.7},
  mean: 0.77809, std: 0.03452, params: {'subsample': 0.6, 'colsample_bytree': 0.8},
  mean: 0.78994, std: 0.03483, params: {'subsample': 0.7, 'colsample_bytree': 0.8},
  mean: 0.79369, std: 0.03232, params: {'subsample': 0.8, 'colsample_bytree': 0.8},
  mean: 0.79207, std: 0.03057, params: {'subsample': 0.9, 'colsample_bytree': 0.8},
  mean: 0.78466, std: 0.02672, params: {'subsample': 0.6, 'colsample_bytree': 0.9},
  mean: 0.78863, std: 0.03289, params: {'subsample': 0.7, 'colsample_bytree': 0.9},
  mean: 0.78905, std: 0.02660, params: {'subsample': 0.8, 'colsample_bytree': 0.9},
  mean: 0.78501, std: 0.03666, params: {'subsample': 0.9, 'colsample_bytree': 0.9}],
 {'colsample_bytree': 0.8, 'subsample': 0.8},
 0.79369231068019075)

gsearch4.grid_scores_, gsearch4.best_params_, gsearch4.best_score_

([mean: 0.78994, std: 0.02779, params: {'subsample': 0.6, 'colsample_bytree': 0.6},

mean: 0.78900, std: 0.03519, params: {'subsample': 0.7, 'colsample_bytree': 0.6},

mean: 0.78509, std: 0.03202, params: {'subsample': 0.8, 'colsample_bytree': 0.6},

mean: 0.78706, std: 0.02848, params: {'subsample': 0.9, 'colsample_bytree': 0.6},

mean: 0.78511, std: 0.03140, params: {'subsample': 0.6, 'colsample_bytree': 0.7},

mean: 0.78343, std: 0.03336, params: {'subsample': 0.7, 'colsample_bytree': 0.7},

mean: 0.78939, std: 0.03203, params: {'subsample': 0.8, 'colsample_bytree': 0.7},

mean: 0.78646, std: 0.04090, params: {'subsample': 0.9, 'colsample_bytree': 0.7},

mean: 0.77809, std: 0.03452, params: {'subsample': 0.6, 'colsample_bytree': 0.8},

mean: 0.78994, std: 0.03483, params: {'subsample': 0.7, 'colsample_bytree': 0.8},

mean: 0.79369, std: 0.03232, params: {'subsample': 0.8, 'colsample_bytree': 0.8},

mean: 0.79207, std: 0.03057, params: {'subsample': 0.9, 'colsample_bytree': 0.8},

mean: 0.78466, std: 0.02672, params: {'subsample': 0.6, 'colsample_bytree': 0.9},

mean: 0.78863, std: 0.03289, params: {'subsample': 0.7, 'colsample_bytree': 0.9},

mean: 0.78905, std: 0.02660, params: {'subsample': 0.8, 'colsample_bytree': 0.9},

mean: 0.78501, std: 0.03666, params: {'subsample': 0.9, 'colsample_bytree': 0.9}],

{'colsample_bytree': 0.8, 'subsample': 0.8},

0.79369231068019075)

より細かい範囲で再度パラメータをチューニングします。

#Grid seach on subsample and max_features
#Choose all predictors except target & IDcols
param_test5 = {
    'subsample':[i/100.0 for i in range(75,90,5)],
    'colsample_bytree':[i/100.0 for i in range(75,90,5)]
}
gsearch5 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=177, max_depth=8,
                                        min_child_weight=1, gamma=0.4, subsample=0.8, colsample_bytree=0.8,
                                        objective= 'binary:logistic', nthread=4, scale_pos_weight=1,seed=27), 
                       param_grid = param_test5, scoring='roc_auc',n_jobs=4,iid=False, cv=5)
gsearch5.fit(train[predictors],train[target])

#Grid seach on subsample and max_features

#Choose all predictors except target & IDcols

param_test5 = {

'subsample':[i/100.0 for i in range(75,90,5)],

'colsample_bytree':[i/100.0 for i in range(75,90,5)]

}

gsearch5 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=177, max_depth=8,

min_child_weight=1, gamma=0.4, subsample=0.8, colsample_bytree=0.8,

objective= 'binary:logistic', nthread=4, scale_pos_weight=1,seed=27),

param_grid = param_test5, scoring='roc_auc',n_jobs=4,iid=False, cv=5)

gsearch5.fit(train[predictors],train[target])

gsearch5.grid_scores_, gsearch5.best_params_, gsearch5.best_score_

([mean: 0.78890, std: 0.03171, params: {'subsample': 0.75, 'colsample_bytree': 0.75},
  mean: 0.79369, std: 0.03232, params: {'subsample': 0.8, 'colsample_bytree': 0.75},
  mean: 0.79374, std: 0.03061, params: {'subsample': 0.85, 'colsample_bytree': 0.75},
  mean: 0.78890, std: 0.03171, params: {'subsample': 0.75, 'colsample_bytree': 0.8},
  mean: 0.79369, std: 0.03232, params: {'subsample': 0.8, 'colsample_bytree': 0.8},
  mean: 0.79374, std: 0.03061, params: {'subsample': 0.85, 'colsample_bytree': 0.8},
  mean: 0.78418, std: 0.03232, params: {'subsample': 0.75, 'colsample_bytree': 0.85},
  mean: 0.78905, std: 0.02660, params: {'subsample': 0.8, 'colsample_bytree': 0.85},
  mean: 0.78367, std: 0.03582, params: {'subsample': 0.85, 'colsample_bytree': 0.85}],
 {'colsample_bytree': 0.75, 'subsample': 0.85},
 0.79374219292158221)

gsearch5.grid_scores_, gsearch5.best_params_, gsearch5.best_score_

([mean: 0.78890, std: 0.03171, params: {'subsample': 0.75, 'colsample_bytree': 0.75},

mean: 0.79369, std: 0.03232, params: {'subsample': 0.8, 'colsample_bytree': 0.75},

mean: 0.79374, std: 0.03061, params: {'subsample': 0.85, 'colsample_bytree': 0.75},

mean: 0.78890, std: 0.03171, params: {'subsample': 0.75, 'colsample_bytree': 0.8},

mean: 0.79369, std: 0.03232, params: {'subsample': 0.8, 'colsample_bytree': 0.8},

mean: 0.79374, std: 0.03061, params: {'subsample': 0.85, 'colsample_bytree': 0.8},

mean: 0.78418, std: 0.03232, params: {'subsample': 0.75, 'colsample_bytree': 0.85},

mean: 0.78905, std: 0.02660, params: {'subsample': 0.8, 'colsample_bytree': 0.85},

mean: 0.78367, std: 0.03582, params: {'subsample': 0.85, 'colsample_bytree': 0.85}],

{'colsample_bytree': 0.75, 'subsample': 0.85},

0.79374219292158221)

続いて、reg_alphaをチューニングします。

#Grid seach on subsample and max_features
#Choose all predictors except target & IDcols
param_test6 = {
    'reg_alpha':[1e-5, 1e-2, 0.1, 1, 100]
}
gsearch6 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=177, max_depth=8,
                                        min_child_weight=1, gamma=0.4, subsample=0.85, colsample_bytree=0.75,
                                        objective= 'binary:logistic', nthread=4, scale_pos_weight=1,seed=27), 
                       param_grid = param_test6, scoring='roc_auc',n_jobs=4,iid=False, cv=5)
gsearch6.fit(train[predictors],train[target])

#Grid seach on subsample and max_features

#Choose all predictors except target & IDcols

param_test6 = {

'reg_alpha':[1e-5, 1e-2, 0.1, 1, 100]

}

gsearch6 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=177, max_depth=8,

min_child_weight=1, gamma=0.4, subsample=0.85, colsample_bytree=0.75,

objective= 'binary:logistic', nthread=4, scale_pos_weight=1,seed=27),

param_grid = param_test6, scoring='roc_auc',n_jobs=4,iid=False, cv=5)

gsearch6.fit(train[predictors],train[target])

gsearch6.grid_scores_, gsearch6.best_params_, gsearch6.best_score_

([mean: 0.79377, std: 0.03058, params: {'reg_alpha': 1e-05},
  mean: 0.79068, std: 0.02953, params: {'reg_alpha': 0.01},
  mean: 0.79298, std: 0.03268, params: {'reg_alpha': 0.1},
  mean: 0.78731, std: 0.03270, params: {'reg_alpha': 1},
  mean: 0.72370, std: 0.03333, params: {'reg_alpha': 100}],
 {'reg_alpha': 1e-05},
 0.79376831622356758)

gsearch6.grid_scores_, gsearch6.best_params_, gsearch6.best_score_

([mean: 0.79377, std: 0.03058, params: {'reg_alpha': 1e-05},

mean: 0.79068, std: 0.02953, params: {'reg_alpha': 0.01},

mean: 0.79298, std: 0.03268, params: {'reg_alpha': 0.1},

mean: 0.78731, std: 0.03270, params: {'reg_alpha': 1},

mean: 0.72370, std: 0.03333, params: {'reg_alpha': 100}],

{'reg_alpha': 1e-05},

0.79376831622356758)

範囲が粗かったので、より細かくパラメータをチューニングします。

#Grid seach on subsample and max_features
#Choose all predictors except target & IDcols
param_test7 = {
    'reg_alpha':[0, 0.001, 0.005, 0.01, 0.05]
}
gsearch7 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=177, max_depth=8,
                                        min_child_weight=1, gamma=0.4, subsample=0.85, colsample_bytree=0.75,
                                        objective= 'binary:logistic', nthread=4, scale_pos_weight=1,seed=27), 
                       param_grid = param_test7, scoring='roc_auc',n_jobs=4,iid=False, cv=5)
gsearch7.fit(train[predictors],train[target])

#Grid seach on subsample and max_features

#Choose all predictors except target & IDcols

param_test7 = {

'reg_alpha':[0, 0.001, 0.005, 0.01, 0.05]

}

gsearch7 = GridSearchCV(estimator = XGBClassifier( learning_rate =0.1, n_estimators=177, max_depth=8,

min_child_weight=1, gamma=0.4, subsample=0.85, colsample_bytree=0.75,

objective= 'binary:logistic', nthread=4, scale_pos_weight=1,seed=27),

param_grid = param_test7, scoring='roc_auc',n_jobs=4,iid=False, cv=5)

gsearch7.fit(train[predictors],train[target])

gsearch7.grid_scores_, gsearch7.best_params_, gsearch7.best_score_

([mean: 0.79374, std: 0.03061, params: {'reg_alpha': 0},
  mean: 0.79433, std: 0.03076, params: {'reg_alpha': 0.001},
  mean: 0.79099, std: 0.02989, params: {'reg_alpha': 0.005},
  mean: 0.79068, std: 0.02953, params: {'reg_alpha': 0.01},
  mean: 0.79160, std: 0.02950, params: {'reg_alpha': 0.05}],
 {'reg_alpha': 0.001},
 0.79432567460197734)

gsearch7.grid_scores_, gsearch7.best_params_, gsearch7.best_score_

([mean: 0.79374, std: 0.03061, params: {'reg_alpha': 0},

mean: 0.79433, std: 0.03076, params: {'reg_alpha': 0.001},

mean: 0.79099, std: 0.02989, params: {'reg_alpha': 0.005},

mean: 0.79068, std: 0.02953, params: {'reg_alpha': 0.01},

mean: 0.79160, std: 0.02950, params: {'reg_alpha': 0.05}],

{'reg_alpha': 0.001},

0.79432567460197734)

これまでにチューニングしてきたパラメータを用いて再度推定を行います。

xgb3 = XGBClassifier(
        learning_rate =0.1,
        n_estimators=1000,
        max_depth=8,
        min_child_weight=1,
        gamma=0.4,
        subsample=0.85,
        colsample_bytree=0.75,
        reg_alpha=0.001,
        objective= 'binary:logistic',
        nthread=4,
        scale_pos_weight=1,
        seed=27)
modelfit(xgb3, train, test, predictors)

Will train until cv error hasn't decreased in 50 rounds.
Stopping. Best iteration: 153

Model Report
Accuracy : 1
AUC Score (Train): 1.000000
AUC Score (Test): 0.880331

xgb3 = XGBClassifier(

learning_rate =0.1,

n_estimators=1000,

max_depth=8,

min_child_weight=1,

gamma=0.4,

subsample=0.85,

colsample_bytree=0.75,

reg_alpha=0.001,

objective= 'binary:logistic',

nthread=4,

scale_pos_weight=1,

seed=27)

modelfit(xgb3, train, test, predictors)

Will train until cv error hasn't decreased in 50 rounds.

Stopping. Best iteration: 153

Model Report

Accuracy : 1

AUC Score (Train): 1.000000

AUC Score (Test): 0.880331

ブログであるように試行回数を1,000回から5,000回まで増やしてみます。

xgb4 = XGBClassifier(
        learning_rate =0.01,
        n_estimators=5000,
        max_depth=8,
        min_child_weight=1,
        gamma=0.4,
        subsample=0.85,
        colsample_bytree=0.75,
        reg_alpha=0.001,
        objective= 'binary:logistic',
        nthread=4,
        scale_pos_weight=1,
        seed=27)
modelfit(xgb4, train, test, predictors)

Will train until cv error hasn't decreased in 50 rounds.
Stopping. Best iteration: 604

Model Report
Accuracy : 0.9951
AUC Score (Train): 0.999955
AUC Score (Test): 0.888000

xgb4 = XGBClassifier(

learning_rate =0.01,

n_estimators=5000,

max_depth=8,

min_child_weight=1,

gamma=0.4,

subsample=0.85,

colsample_bytree=0.75,

reg_alpha=0.001,

objective= 'binary:logistic',

nthread=4,

scale_pos_weight=1,

seed=27)

modelfit(xgb4, train, test, predictors)

Will train until cv error hasn't decreased in 50 rounds.

Stopping. Best iteration: 604

Model Report

Accuracy : 0.9951

AUC Score (Train): 0.999955

AUC Score (Test): 0.888000

88.8%まで向上しました。色々と数値いじっても、1%高めるだけにとどまってしまうのですね。

とにかく、XGBoostをPythonで実行してパラメータチューニングするという一連の試行がこのコードでできそうなので、今後も使いまわしてみようと思います。

不均衡なデータの分類問題について with Python

データマイニング界隈で人気のKDnuggetsで紹介されていた、”Dealing with Unbalanced Classes, SVMs, Random Forests, and Decision Trees in Python“のプログラムが残念なことに画像だったので、写経しました。せっかくなので、紹介させていただきます。内容としては不均衡データに対する処方の紹介で、プログラムはPythonで書かれています。ライブラリさえインストールできれば皆さんもすぐに実行できるので、是非チャレンジしてみて下さい。

まずはもろもろライブラリを呼び出します。

%matplotlib inline
import numpy as np
import scipy as sp
import pandas as pd
import sklearn
import seaborn as sns
from matplotlib import  pyplot as plt

import sklearn.cross_validation

%matplotlib inline

import numpy as np

import scipy as sp

import pandas as pd

import sklearn

import seaborn as sns

from matplotlib import pyplot as plt

import sklearn.cross_validation

CSV形式のデータセットをWebサイトから取得します。ワインの評価と、ワインに関した特徴量からなるデータセットです。

wine_df = pd.read_csv("https://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-red.csv", sep=";")
wine_df.head()

1 2	wine_df = pd.read_csv("https://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-red.csv", sep=";") wine_df.head()

分類のための目的変数を作成します。

#ワインの質に関する数値
Y = wine_df.quality.values
#質に関するデータを落としている。
wine_df = wine_df.drop('quality', axis =1)
#7よりも小さかったら0、それ以外は1とする。
Y = np.asarray([1 if  i>=7 else 0 for i in Y])
wine_df.head()

#ワインの質に関する数値

Y = wine_df.quality.values

#質に関するデータを落としている。

wine_df = wine_df.drop('quality', axis =1)

#7よりも小さかったら0、それ以外は1とする。

Y = np.asarray([1 if i>=7 else 0 for i in Y])

wine_df.head()

X =  wine_df.as_matrix()

1	X = wine_df.as_matrix()

ランダムフォレストを実行します。

from sklearn.ensemble import RandomForestClassifier
from sklearn.cross_validation import  cross_val_score

scores =[]

#1~41までの木の数のランダムフォレストを実行する。
for val in range(1,41):
    clf = RandomForestClassifier(n_estimators  =val)
    validated = cross_val_score(clf, X, Y, cv =10)
    scores.append(validated)

from sklearn.ensemble import RandomForestClassifier

from sklearn.cross_validation import cross_val_score

scores =[]

#1~41までの木の数のランダムフォレストを実行する。

for val in range(1,41):

clf = RandomForestClassifier(n_estimators =val)

validated = cross_val_score(clf, X, Y, cv =10)

scores.append(validated)

どんな結果が返ってくるのか、試しに一つだけツリーの数を2にして実行してみます。10回分のクロスバリデーションを行った推定結果が出力されています。これは、いわゆる正解率のことを指します。

#木の数が2のランダムフォレストの結果を返す
clf1 = RandomForestClassifier(n_estimators = 2)
validated = cross_val_score(clf1, X, Y, cv=10)
validated

array([ 0.8757764 ,  0.86956522,  0.8625    ,  0.84375   ,  0.8875    ,
        0.85      ,  0.8375    ,  0.85534591,  0.85534591,  0.88679245])

#木の数が2のランダムフォレストの結果を返す

clf1 = RandomForestClassifier(n_estimators = 2)

validated = cross_val_score(clf1, X, Y, cv=10)

validated

array([ 0.8757764 , 0.86956522, 0.8625 , 0.84375 , 0.8875 ,

0.85 , 0.8375 , 0.85534591, 0.85534591, 0.88679245])

ツリーの数に応じた正解率を可視化します。

sns.boxplot(data=scores)
plt.xlabel('number of trees')
plt.ylabel('Classification scores')
plt.title('Classification score for number of trees')
plt.show()

sns.boxplot(data=scores)

plt.xlabel('number of trees')

plt.ylabel('Classification scores')

plt.title('Classification score for number of trees')

plt.show()

正解率はツリーの数を増やすことで増すようです。

しかしながら、正解率は誤解されやすい指標です。不均衡データでは偏りのある方ばかりを当てていても、正解率は増してしまいます。当たりだけを予測できて、ハズレを予測できないというのは分類器として使いみちが限られると思います。そこで、悪いワインの割合を直線で引いてみます。

len_y = len(Y)
temp = [i for i in Y if i ==0]
temp_1 = temp.count(0)

percentage = float(temp_1)/float(len_y)

print(float(temp_1)/float(len_y)*100)

sns.boxplot(data=scores)
plt.axhline(y = percentage, ls = '--')
plt.xlabel('number of trees')
plt.ylabel('Classification Scores')
plt.title('Classification scores of  for trees')
plt.show()

len_y = len(Y)

temp = [i for i in Y if i ==0]

temp_1 = temp.count(0)

percentage = float(temp_1)/float(len_y)

print(float(temp_1)/float(len_y)*100)

sns.boxplot(data=scores)

plt.axhline(y = percentage, ls = '--')

plt.xlabel('number of trees')

plt.ylabel('Classification Scores')

plt.title('Classification scores of for trees')

plt.show()

悪いワインの割合がそもそも多いので、悪いワインと判定しまくっていても、正解率は高いわけです。

そこで、機械学習における予測精度の評価指標とされているF値を使います。
ツリーの数を増やしても、F値は良くなっていないようです。

scores = []

for val in range(1, 41):
    cfl = RandomForestClassifier(n_estimators = val)
    validated = cross_val_score(clf, X, Y, cv=10, scoring = 'f1')
    scores.append(validated)

scores = []

for val in range(1, 41):

cfl = RandomForestClassifier(n_estimators = val)

validated = cross_val_score(clf, X, Y, cv=10, scoring = 'f1')

scores.append(validated)

sns.boxplot( data=scores)
plt.xlabel('number of trees')
plt.ylabel('F1 Scores')
plt.title('F1 scores as a function of the number of trees')
plt.show()

sns.boxplot( data=scores)

plt.xlabel('number of trees')

plt.ylabel('F1 Scores')

plt.title('F1 scores as a function of the number of trees')

plt.show()

ここでは、0.5よりも大きいとする予測になる特徴量のデータを切り捨てます。その切り捨てる割合がどこが望ましいのかを以下で探していきます。

clf = RandomForestClassifier(n_estimators= 15)
clf.fit(X, Y)

(clf.predict_proba(X)[:,1] > 0.5).astype(int)

clf = RandomForestClassifier(n_estimators= 15)

clf.fit(X, Y)

(clf.predict_proba(X)[:,1] > 0.5).astype(int)

def cutoff_predict(clf, X,  cutoff):
    return (clf.predict_proba(X)[:,1] > cutoff).astype(int)

scores = []

def custom_f1(cutoff):
    def f1_cutoff(clf, X, Y):
        ypred = cutoff_predict(clf, X,  cutoff)
        return sklearn.metrics.f1_score(Y, ypred)
    
    return f1_cutoff

for cutoff in np.arange(0.1,  0.9, 0.1):
    clf = RandomForestClassifier(n_estimators=15)
    validated = cross_val_score(clf,  X, Y, cv=10, scoring=custom_f1(cutoff)) 
    scores.append(validated)

def cutoff_predict(clf, X, cutoff):

return (clf.predict_proba(X)[:,1] > cutoff).astype(int)

scores = []

def custom_f1(cutoff):

def f1_cutoff(clf, X, Y):

ypred = cutoff_predict(clf, X, cutoff)

return sklearn.metrics.f1_score(Y, ypred)

return f1_cutoff

for cutoff in np.arange(0.1, 0.9, 0.1):

clf = RandomForestClassifier(n_estimators=15)

validated = cross_val_score(clf, X, Y, cv=10, scoring=custom_f1(cutoff))

scores.append(validated)

sns.boxplot(data=scores, names= np.arange(0.1, 0.9, 0.1))
plt.xlabel('each cut off value')
plt.ylabel('F1 Scores')
plt.title('custom F scores')
plt.show()

sns.boxplot(data=scores, names= np.arange(0.1, 0.9, 0.1))

plt.xlabel('each cut off value')

plt.ylabel('F1 Scores')

plt.title('custom F scores')

plt.show()

どうやら、階級値が2～4、つまり割合にして0.3～0.5のカットオフ値が望ましい水準のようです。

以下では、決定境界の可視化を行います。しかしながら、二次元の可視化となると、複数あるデータの中から特徴量を二つだけ選ばなければなりません。その変数を決めるに際して、変数の重要度を用います。変数の重要度はランダムフォレストで計算可能です。

clf = RandomForestClassifier(n_estimators=15)
clf.fit(X, Y)

imp = clf.feature_importances_
names = wine_df.columns

imp, names = zip(*sorted(zip(imp, names)))

plt.barh(range(len(names)), imp, align='center')
plt.yticks(range(len(names)), names)

plt.xlabel('Importance of features')
plt.ylabel('Features')
plt.title('Importance of each feature')
plt.show()

clf = RandomForestClassifier(n_estimators=15)

clf.fit(X, Y)

imp = clf.feature_importances_

names = wine_df.columns

imp, names = zip(*sorted(zip(imp, names)))

plt.barh(range(len(names)), imp, align='center')

plt.yticks(range(len(names)), names)

plt.xlabel('Importance of features')

plt.ylabel('Features')

plt.title('Importance of each feature')

plt.show()

from sklearn.tree import DecisionTreeClassifier
import sklearn.linear_model
import sklearn.svm

def plot_decision_surface(clf, X_train, Y_train):
    plot_step=0.1
    
    if X_train.shape[1] != 2:
        raise ValueError("X_train should have exactly 2 columns!")
        
    x_min, x_max = X_train[:, 0].min() - plot_step, X_train[:, 0].max() + plot_step
    y_min, y_max = X_train[:, 1].min() - plot_step, X_train[:, 1].max() + plot_step
    xx, yy = np.meshgrid(np.arange(x_min, x_max, plot_step),
                        np.arange(y_min, y_max, plot_step))
    
    clf.fit(X_train, Y_train)
    if hasattr(clf, 'predict_proba'):
        Z = clf.predict_proba(np.c_[xx.ravel(), yy.ravel()])[:,1]
    else:
        Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])
    Z = Z.reshape(xx.shape)
    cs = plt.contourf(xx, yy, Z, cmap = plt.cm.Reds)
    plt.scatter(X_train[:,0], X_train[:,1], c=Y_train, cmap=plt.cm.Paired)
    plt.show()

from sklearn.tree import DecisionTreeClassifier

import sklearn.linear_model

import sklearn.svm

def plot_decision_surface(clf, X_train, Y_train):

plot_step=0.1

if X_train.shape[1] != 2:

raise ValueError("X_train should have exactly 2 columns!")

x_min, x_max = X_train[:, 0].min() - plot_step, X_train[:, 0].max() + plot_step

y_min, y_max = X_train[:, 1].min() - plot_step, X_train[:, 1].max() + plot_step

xx, yy = np.meshgrid(np.arange(x_min, x_max, plot_step),

np.arange(y_min, y_max, plot_step))

clf.fit(X_train, Y_train)

if hasattr(clf, 'predict_proba'):

Z = clf.predict_proba(np.c_[xx.ravel(), yy.ravel()])[:,1]

else:

Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])

Z = Z.reshape(xx.shape)

cs = plt.contourf(xx, yy, Z, cmap = plt.cm.Reds)

plt.scatter(X_train[:,0], X_train[:,1], c=Y_train, cmap=plt.cm.Paired)

plt.show()

こちらで、重要度が上位のものに絞って決定境界を可視化します。ここでは、ランダムフォレストのみならず、SVMや決定木も実行されています。

imp_fe = np.argsort(imp)[::-1][0:2]
X_imp = X[:, imp_fe]

algorithms = [DecisionTreeClassifier(),
             RandomForestClassifier(),
             sklearn.svm.SVC(C = 100.0, gamma = 1)]

title = ['Decision Tree Classifier', 'Random Forest Classifier',
        'Support Vector Maachine']

for i in xrange(3):
    plt.title(title[i])
    plt.xlabel('Feature1')
    plt.ylabel('Feature2')
    plot_decision_surface(algorithms[i], X_imp, Y)

imp_fe = np.argsort(imp)[::-1][0:2]

X_imp = X[:, imp_fe]

algorithms = [DecisionTreeClassifier(),

RandomForestClassifier(),

sklearn.svm.SVC(C = 100.0, gamma = 1)]

title = ['Decision Tree Classifier', 'Random Forest Classifier',

'Support Vector Maachine']

for i in xrange(3):

plt.title(title[i])

plt.xlabel('Feature1')

plt.ylabel('Feature2')

plot_decision_surface(algorithms[i], X_imp, Y)

sklearnのSVMはデフォルトではクラスごとの重み付けを行わないが、自動で重み付けを行うことが出来る。以下の例では、C=1、gamma=1でクラスごとの重み付けを行う・行わないでの決定境界を描いている。重み付けを行うことで、赤色の少ない方のデータの識別が比較的できていることが伺えるが、他方で、多くの青を誤判定している。さらなる改善にはパラメータチューニングが必要となります。

svm = [sklearn.svm.SVC(C = 1.0, gamma = 1.0, class_weight=None),
      sklearn.svm.SVC(C = 1.0, gamma = 1.0, class_weight='auto')]

title = ['Svm without class weight', 'Svm with class weight']

for i in xrange(2):
    plt.title(title[i])
    plt.xlabel('Feature1')
    plt.ylabel('Feature2')
    
    plot_decision_surface(svm[i], X_imp, Y)

svm = [sklearn.svm.SVC(C = 1.0, gamma = 1.0, class_weight=None),

sklearn.svm.SVC(C = 1.0, gamma = 1.0, class_weight='auto')]

title = ['Svm without class weight', 'Svm with class weight']

for i in xrange(2):

plt.title(title[i])

plt.xlabel('Feature1')

plt.ylabel('Feature2')

plot_decision_surface(svm[i], X_imp, Y)

不均衡データに対するアプローチや、Pythonによる機械学習を学ぶ良い機会になりました。KDnuggetsは非常に勉強になりますね。

Tokyo.R#53で得たパッケージ情報とその実践

第53回のTokyo.Rで気になったパッケージの情報と実行例をいくつかあげました。スライドなどもろもろの発表はこちらの方のブログ「第53回R勉強会@東京で発表してきた」が非常に詳しく書かれています。

【目次】
・ggradarパッケージ
・proxyパッケージ
・因果推論（CBPSパッケージ）
・MXNetパッケージ
・missForestパッケージ
・RFinanceパッケージ

ggradarパッケージ

簡単にレーダーチャートを作れるパッケージです。こちらのブログを参考にしています。

install.packages("devtools")
devtools::install_github("ricardo-bion/ggradar")

1 2	install.packages("devtools") devtools::install_github("ricardo-bion/ggradar")

企業の職場環境に関してまとめられた某口コミサイトから4個ほどデータを拝借してきました。

> CompanyVoiceData
  company growth stability salary rewarding idea difficulty welfare education
1  google    5.0       5.0    4.9       5.0  4.3        5.0     5.0       4.6
2   yahoo    3.9       5.0    3.2       3.8  3.7        3.9     3.1       3.3
3 recruit    4.4       4.8    5.0       5.0  5.0        5.0     4.0       5.0
4  amazon    5.0       5.0    4.2       4.0  4.2        5.0     3.6       3.3

> CompanyVoiceData

company growth stability salary rewarding idea difficulty welfare education

1 google 5.0 5.0 4.9 5.0 4.3 5.0 5.0 4.6

2 yahoo 3.9 5.0 3.2 3.8 3.7 3.9 3.1 3.3

3 recruit 4.4 4.8 5.0 5.0 5.0 5.0 4.0 5.0

4 amazon 5.0 5.0 4.2 4.0 4.2 5.0 3.6 3.3

ggradarをそのまま使おうとすると、Circular Air Lightというフォントが必要だと怒られるので、参考のブログにある通り、OSXの場合はこちらをダブルクリックでインストールして再起動します。

先ほどのデータに対して、以下のコードを実行すれば非常に簡単にレーダーチャートが作れました。

library("ggradar")
CompanyVoiceData <- data.frame(read.csv(file ="company_voice.csv",header = TRUE))

ggradar(CompanyVoiceData, 
        grid.max = max(CompanyVoiceData[, 2:ncol(CompanyVoiceData)]),
        background.circle.colour = "#ffdd99", #背景色の指定
        background.circle.transparency = 1, #背景色の透明度を指定
        group.line.width = 2, #線の太さの指定
        group.point.size = 6, #シンボルの大きさの指定
        axis.label.size = 5, #軸ラベルサイズの指定
        gridline.min.colour = "#4b61ba", #最小円の線色の指定
        gridline.mid.colour = "#a87963", #中円の線色の指定
        gridline.max.colour = "#e1e6ea", #最大円の線色の指定
        grid.line.width = 1.5, #各円の線の太さの指定
        gridline.min.linetype = "longdash", #線種の指定
        gridline.mid.linetype = "longdash", #線種の指定
        gridline.max.linetype = "longdash") #線種の指定

library("ggradar")

CompanyVoiceData <- data.frame(read.csv(file ="company_voice.csv",header = TRUE))

ggradar(CompanyVoiceData,

grid.max = max(CompanyVoiceData[, 2:ncol(CompanyVoiceData)]),

background.circle.colour = "#ffdd99", #背景色の指定

background.circle.transparency = 1, #背景色の透明度を指定

group.line.width = 2, #線の太さの指定

group.point.size = 6, #シンボルの大きさの指定

axis.label.size = 5, #軸ラベルサイズの指定

gridline.min.colour = "#4b61ba", #最小円の線色の指定

gridline.mid.colour = "#a87963", #中円の線色の指定

gridline.max.colour = "#e1e6ea", #最大円の線色の指定

grid.line.width = 1.5, #各円の線の太さの指定

gridline.min.linetype = "longdash", #線種の指定

gridline.mid.linetype = "longdash", #線種の指定

gridline.max.linetype = "longdash") #線種の指定

proxyパッケージ

距離や類似度を計算するパッケージです。
先ほどのデータに対して類似度と距離を計算してみます。

library(proxy)
> simil(CompanyVoiceData[,-1])
          1         2         3
2 0.2286639                    
3 0.6373648 0.1339713          
4 0.6499133 0.5787506 0.4188571
> dist(CompanyVoiceData[,-1])
         1        2        3
2 3.522783                  
3 1.435270 3.401470         
4 2.269361 1.989975 2.393742

library(proxy)

> simil(CompanyVoiceData[,-1])

1 2 3

2 0.2286639

3 0.6373648 0.1339713

4 0.6499133 0.5787506 0.4188571

> dist(CompanyVoiceData[,-1])

1 2 3

2 3.522783

3 1.435270 3.401470

4 2.269361 1.989975 2.393742

こんな感じで、類似度や距離の計算ができます。

因果推論

こちらはパッケージとかそういうものではなく、既存の関数などで計算できるようです。
こちらのブログ、「調査観察データにおける因果推論(3) – Rによる傾向スコア，IPW推定量，二重にロバストな推定量の算出」に詳しく書かれています。
・glm関数での傾向スコアの算出
・傾向スコアを共変量としてlm関数で回帰分析
・コードを愚直に書いてIPW推定量の算出
・期待値の標準誤差を出すための関数を作成
・DR推定量の算出をするための関数を作成
などで、推定自体は実現できるようです。

ただし、CBPS(Covariate Balancing Propensity Score)というパッケージがあるらしく、このパッケージを用いれば因果推論の計算を行えるようです。

Package ‘CBPS’
以下のようなExampleコードが載っていたので、実行してみましたが、なかなか結果が返ってこなかったので不安になりました。計算が終わるまで10分以上はかかったと思います。

library(CBPS)
data(Blackwell)

form1<-"d.gone.neg ~ d.gone.neg.l1 + d.gone.neg.l2 + d.neg.frac.l3 + camp.length + camp.length +
deminc + base.poll + year.2002 + year.2004 + year.2006 + base.und + office"

##Fitting the models in Imai and Ratkovic (2014)
##Warning: may take a few mintues; setting time.vary to FALSE
##Results in a quicker fit but with poorer balance
fit1 <- CBMSM(formula = form1, time=Blackwell$time,id=Blackwell$demName,data=Blackwell, type="MSM",
            iterations = NULL, twostep = TRUE, msm.variance = "full", time.vary = TRUE)
fit2 <- CBMSM(formula = form1, time=Blackwell$time,id=Blackwell$demName,data=Blackwell, type="MSM",
            iterations = NULL, twostep = TRUE, msm.variance = "approx", time.vary = TRUE)

##Assessing balance
bal1 <- balance(fit1)
bal2 <- balance(fit2)

##Effect estimation: Replicating Effect Estimates in
##Table 3 of Imai and Ratkovic (2014)
lm1 <- lm(demprcnt[time==1]~fit1$treat.hist,data=Blackwell,weights=fit1$glm.weights)
lm2 <- lm(demprcnt[time==1]~fit1$treat.hist,data=Blackwell,weights=fit1$weights)
lm3 <- lm(demprcnt[time==1]~fit1$treat.hist,data=Blackwell,weights=fit2$weights)
lm4 <- lm(demprcnt[time==1]~fit1$treat.cum,data=Blackwell,weights=fit1$glm.weights)
lm5 <- lm(demprcnt[time==1]~fit1$treat.cum,data=Blackwell,weights=fit1$weights)
lm6 <- lm(demprcnt[time==1]~fit1$treat.cum,data=Blackwell,weights=fit2$weights)

library(CBPS)

data(Blackwell)

form1<-"d.gone.neg ~ d.gone.neg.l1 + d.gone.neg.l2 + d.neg.frac.l3 + camp.length + camp.length +

deminc + base.poll + year.2002 + year.2004 + year.2006 + base.und + office"

##Fitting the models in Imai and Ratkovic (2014)

##Warning: may take a few mintues; setting time.vary to FALSE

##Results in a quicker fit but with poorer balance

fit1 <- CBMSM(formula = form1, time=Blackwell$time,id=Blackwell$demName,data=Blackwell, type="MSM",

iterations = NULL, twostep = TRUE, msm.variance = "full", time.vary = TRUE)

fit2 <- CBMSM(formula = form1, time=Blackwell$time,id=Blackwell$demName,data=Blackwell, type="MSM",

iterations = NULL, twostep = TRUE, msm.variance = "approx", time.vary = TRUE)

##Assessing balance

bal1 <- balance(fit1)

bal2 <- balance(fit2)

##Effect estimation: Replicating Effect Estimates in

##Table 3 of Imai and Ratkovic (2014)

lm1 <- lm(demprcnt[time==1]~fit1$treat.hist,data=Blackwell,weights=fit1$glm.weights)

lm2 <- lm(demprcnt[time==1]~fit1$treat.hist,data=Blackwell,weights=fit1$weights)

lm3 <- lm(demprcnt[time==1]~fit1$treat.hist,data=Blackwell,weights=fit2$weights)

lm4 <- lm(demprcnt[time==1]~fit1$treat.cum,data=Blackwell,weights=fit1$glm.weights)

lm5 <- lm(demprcnt[time==1]~fit1$treat.cum,data=Blackwell,weights=fit1$weights)

lm6 <- lm(demprcnt[time==1]~fit1$treat.cum,data=Blackwell,weights=fit2$weights)

MXNet

XGBoostのパッケージを作ったチームが手がけているパッケージで、深層学習を実行できます。

インストール方法はここに書かれています。
Deep Learning for R

install.packages("drat", repos="https://cran.rstudio.com")
drat:::addRepo("dmlc")
install.packages("mxnet")

install.packages("drat", repos="https://cran.rstudio.com")

drat:::addRepo("dmlc")

install.packages("mxnet")

あれ、OSXではエラーが返ってきてライブラリが読み込めないですね。どうやら私のためにあるようなブログ「Installing mxnet for R on Yosemite」があったので、時間を見つけてチャレンジしてみようと思います。

ディープラーニングを用いた回帰分析については、Neural Network with MXNet in Five Minutesにコードがもろもろ載っていますので、チャレンジしてみると良いと思います。

リンク先に載っているのですが、一応コードを以下に記しておきます。

data(BostonHousing, package="mlbench")

train.ind = seq(1, 506, 3)
train.x = data.matrix(BostonHousing[train.ind, -14])
train.y = BostonHousing[train.ind, 14]
test.x = data.matrix(BostonHousing[-train.ind, -14])
test.y = BostonHousing[-train.ind, 14]

# Define the input data
data <- mx.symbol.Variable("data")
# A fully connected hidden layer
# data: input source
# num_hidden: number of neurons in this hidden layer
fc1 <- mx.symbol.FullyConnected(data, num_hidden=1)

# Use linear regression for the output layer
lro <- mx.symbol.LinearRegressionOutput(fc1)

preds = predict(model, test.x)

## Auto detect layout of input matrix, use rowmajor..
sqrt(mean((preds-test.y)^2))

demo.metric.mae <- mx.metric.custom("mae", function(label, pred) {
  res <- mean(abs(label-pred))
  return(res)
})

mx.set.seed(0)
model <- mx.model.FeedForward.create(lro, X=train.x, y=train.y,
                                     ctx=mx.cpu(), num.round=50, array.batch.size=20,
                                     learning.rate=2e-6, momentum=0.9, eval.metric=demo.metric.mae)

data(BostonHousing, package="mlbench")

train.ind = seq(1, 506, 3)

train.x = data.matrix(BostonHousing[train.ind, -14])

train.y = BostonHousing[train.ind, 14]

test.x = data.matrix(BostonHousing[-train.ind, -14])

test.y = BostonHousing[-train.ind, 14]

# Define the input data

data <- mx.symbol.Variable("data")

# A fully connected hidden layer

# data: input source

# num_hidden: number of neurons in this hidden layer

fc1 <- mx.symbol.FullyConnected(data, num_hidden=1)

# Use linear regression for the output layer

lro <- mx.symbol.LinearRegressionOutput(fc1)

preds = predict(model, test.x)

## Auto detect layout of input matrix, use rowmajor..

sqrt(mean((preds-test.y)^2))

demo.metric.mae <- mx.metric.custom("mae", function(label, pred) {

res <- mean(abs(label-pred))

return(res)

})

mx.set.seed(0)

model <- mx.model.FeedForward.create(lro, X=train.x, y=train.y,

ctx=mx.cpu(), num.round=50, array.batch.size=20,

learning.rate=2e-6, momentum=0.9, eval.metric=demo.metric.mae)

missForest

ランダムフォレストを用いて、欠損値補完を行うためのパッケージです。目的変数が欠損していても適用できるようです。
詳しくは、スライドを見ていただいた方がいいですが、以下のプログラムで実行できました。ちなみにスライドはこちら、「Imputation of Missing Values using Random Forest」

library(missForest)
library(dplyr)

#ggplot2のデータセットを読み込む
data(diamonds, package = "ggplot2")
dia.sample <- sample_n(diamonds, size=2000)
dia.sample <- as.data.frame(dia.sample)

#既存データセットに5%の欠損を与える
dia.mis <- prodNA(dia.sample, noNA=0.05)

#補完の実行
dia.imp <- missForest(dia.mis, verbose=TRUE)
dia.imp %>% str(max.level=1)

#補完精度の推定
dia.imp$OOBerror
dia.imp <- missForest(dia.mis, verbose=TRUE, variablewise=TRUE)

#補完精度の検証
mixError(ximp = dia.imp$ximp,
         xmis = dia.mis,
         xtrue = dia.sample)

library(missForest)

library(dplyr)

#ggplot2のデータセットを読み込む

data(diamonds, package = "ggplot2")

dia.sample <- sample_n(diamonds, size=2000)

dia.sample <- as.data.frame(dia.sample)

#既存データセットに5%の欠損を与える

dia.mis <- prodNA(dia.sample, noNA=0.05)

#補完の実行

dia.imp <- missForest(dia.mis, verbose=TRUE)

dia.imp %>% str(max.level=1)

#補完精度の推定

dia.imp$OOBerror

dia.imp <- missForest(dia.mis, verbose=TRUE, variablewise=TRUE)

#補完精度の検証

mixError(ximp = dia.imp$ximp,

xmis = dia.mis,

xtrue = dia.sample)

RFinanceYJ

Yohei Sato, Nobuaki Oshiro, Shinichi Takayanagiさんたちが作った、Yahoo!ファイナンスの株価データを取得できるパッケージです。だいぶ前からあったようですが、使って分析している人は初めて見ました。どうやらYahoo!ファイナンスの仕様によって書き換えていかないといけないようです。「2015-01-20 Rでチャートを書いてみる(9)」のブログに実行可能なプログラムがあります。以下、実行可能なコードを転載いたします。

library(RFinanceYJ)

#API
quoteStockTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')
{
  time.interval <- substr(time.interval,1,1)
  function.stock <- function(quote.table.item){
    if( xmlSize(quote.table.item) < 5) return(NULL) 
    d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)
    o <- as.number(xmlValue(quote.table.item[[2]]))
    h <- as.number(xmlValue(quote.table.item[[3]]))
    l <- as.number(xmlValue(quote.table.item[[4]]))
    c <- as.number(xmlValue(quote.table.item[[5]]))
    v <- ifelse(xmlSize(quote.table.item) >= 6,as.number(xmlValue(quote.table.item[[6]])),0)
    a <- ifelse(xmlSize(quote.table.item) >= 7,as.number(xmlValue(quote.table.item[[7]])),0)
    return(data.frame(date=d,open=o,high=h,low=l,close=c,volume=v, adj_close=a))
  }
  return(quoteTsData(x,function.stock,since,start.num,date.end,time.interval,type="stock"))
}
quoteFundTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')
{
  time.interval <- substr(time.interval,1,1)
  function.fund <- function(quote.table.item){
    d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)
    if(time.interval=='monthly'){
      d <- endOfMonth(d)
    }
    c <- as.number(xmlValue(quote.table.item[[2]]))
    v <- as.number(xmlValue(quote.table.item[[3]]))
    return(data.frame(date=d,constant.value=c,NAV=v))
  }
  return(quoteTsData(x,function.fund,since,start.num,date.end,time.interval,type="fund"))
}
quoteFXTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')
{
  time.interval <- substr(time.interval,1,1)
  function.fx <- function(quote.table.item){
    d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)
    o <- as.number(xmlValue(quote.table.item[[2]]))
    h <- as.number(xmlValue(quote.table.item[[3]]))
    l <- as.number(xmlValue(quote.table.item[[4]]))
    c <- as.number(xmlValue(quote.table.item[[5]]))
    return(data.frame(date=d,open=o,high=h,low=l,close=c))
  }
  return(quoteTsData(x,function.fx,since,start.num,date.end,time.interval,type="fx"))
}
######  private functions  #####
#get time series data from Yahoo! Finance.
quoteTsData <- function(x,function.financialproduct,since,start.num,date.end,time.interval,type="stock"){
  r <- NULL
  result.num <- 51
  financial.data <- data.frame(NULL)
  #start <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&c=\\1&a=\\2&b=\\3",since))
  #end   <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&f=\\1&d=\\2&e=\\3",date.end))
  start <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&sy=\\1&sm=\\2&sd=\\3",since))
  end   <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&ey=\\1&em=\\2&ed=\\3",date.end))
  
  if(!any(time.interval==c('d','w','m'))) stop("Invalid time.interval value")
  
  extractQuoteTable <- function(r,type){
    if(type %in% c("fund","fx")){
      tbl <- r[[2]][[2]][[7]][[3]][[3]][[9]][[2]]
    }
    else{
      tbl <- r[[2]][[2]][[7]][[3]][[3]][[10]][[2]]
    }
    return(tbl)
  }
  
  #while( result.num >= 51 ){
  while(1){
    start.num <- start.num + 1
    quote.table <- NULL
    quote.url <- paste('http://info.finance.yahoo.co.jp/history/?code=',x,start,end,'&p=',start.num,'&tm=',substr(time.interval,1,1),sep="")
    #cat(quote.url)
    #try( r <- xmlRoot(htmlTreeParse(quote.url,error=xmlErrorCumulator(immediate=F))), TRUE)  # これだと取得時にエラーが出た。。
    try(r<-htmlParse(quote.url))
    if( is.null(r) ) stop(paste("Can not access :", quote.url))
    
    #try( quote.table <- r[[2]][[1]][[1]][[16]][[1]][[1]][[1]][[4]][[1]][[1]][[1]], TRUE )
    #try( quote.table <- extractQuoteTable(r,type), TRUE )
    try( quote.table <- xpathApply(r,"//table")[[2]], TRUE )
    
    quote.size<-xmlSize(quote.table)
    #cat(paste("size:",quote.size))
    if(xmlSize(quote.table)<=1){
      return (financial.data)
    }
    if( is.null(quote.table) ){
      if( is.null(financial.data) ){
        stop(paste("Can not quote :", x))
      }else{
        financial.data <- financial.data[order(financial.data$date),]
        return(financial.data)
      }
    }
    
    size <- xmlSize(quote.table)
    for(i in 2:size){
      financial.data <- rbind(financial.data,function.financialproduct(quote.table[[i]]))
    }
    
    #result.num <- xmlSize(quote.table)
    Sys.sleep(1)
  }
  financial.data <- financial.data[order(financial.data$date),]
  return(financial.data)  
}
#convert string formart date to POSIXct object
convertToDate <- function(date.string,time.interval)
{
  #data format is different between monthly and dialy or weekly
  if(any(time.interval==c('d','w'))){
    result <- gsub("^([0-9]{4})([^0-9]+)([0-9]{1,2})([^0-9]+)([0-9]{1,2})([^0-9]+)","\\1-\\3-\\5",date.string)
  }else if(time.interval=='m'){
    result <- gsub("^([0-9]{4})([^0-9]+)([0-9]{1,2})([^0-9]+)","\\1-\\3-01",date.string)
  }
  return(as.POSIXct(result))
}
#convert string to number.
as.number <- function(string)
{
  return(as.double(as.character(gsub("[^0-9.]", "",string))))
}
#return end of month date.
endOfMonth <- function(date.obj)
{
  startOfMonth     <- as.Date(format(date.obj,"%Y%m01"),"%Y%m%d")
  startOfNextMonth <- as.Date(format(startOfMonth+31,"%Y%m01"),"%Y%m%d")
  return(startOfNextMonth-1)
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

library(RFinanceYJ)

#API

quoteStockTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')

{

time.interval <- substr(time.interval,1,1)

function.stock <- function(quote.table.item){

if( xmlSize(quote.table.item) < 5) return(NULL)

d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)

o <- as.number(xmlValue(quote.table.item[[2]]))

h <- as.number(xmlValue(quote.table.item[[3]]))

l <- as.number(xmlValue(quote.table.item[[4]]))

c <- as.number(xmlValue(quote.table.item[[5]]))

v <- ifelse(xmlSize(quote.table.item) >= 6,as.number(xmlValue(quote.table.item[[6]])),0)

a <- ifelse(xmlSize(quote.table.item) >= 7,as.number(xmlValue(quote.table.item[[7]])),0)

return(data.frame(date=d,open=o,high=h,low=l,close=c,volume=v, adj_close=a))

}

return(quoteTsData(x,function.stock,since,start.num,date.end,time.interval,type="stock"))

}

quoteFundTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')

{

time.interval <- substr(time.interval,1,1)

function.fund <- function(quote.table.item){

d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)

if(time.interval=='monthly'){

d <- endOfMonth(d)

}

c <- as.number(xmlValue(quote.table.item[[2]]))

v <- as.number(xmlValue(quote.table.item[[3]]))

return(data.frame(date=d,constant.value=c,NAV=v))

}

return(quoteTsData(x,function.fund,since,start.num,date.end,time.interval,type="fund"))

}

quoteFXTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')

{

time.interval <- substr(time.interval,1,1)

function.fx <- function(quote.table.item){

d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)

o <- as.number(xmlValue(quote.table.item[[2]]))

h <- as.number(xmlValue(quote.table.item[[3]]))

l <- as.number(xmlValue(quote.table.item[[4]]))

c <- as.number(xmlValue(quote.table.item[[5]]))

return(data.frame(date=d,open=o,high=h,low=l,close=c))

}

return(quoteTsData(x,function.fx,since,start.num,date.end,time.interval,type="fx"))

}

###### private functions #####

#get time series data from Yahoo! Finance.

quoteTsData <- function(x,function.financialproduct,since,start.num,date.end,time.interval,type="stock"){

r <- NULL

result.num <- 51

financial.data <- data.frame(NULL)

#start <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&c=\\1&a=\\2&b=\\3",since))

#end <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&f=\\1&d=\\2&e=\\3",date.end))

start <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&sy=\\1&sm=\\2&sd=\\3",since))

end <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&ey=\\1&em=\\2&ed=\\3",date.end))

if(!any(time.interval==c('d','w','m'))) stop("Invalid time.interval value")

extractQuoteTable <- function(r,type){

if(type %in% c("fund","fx")){

tbl <- r[[2]][[2]][[7]][[3]][[3]][[9]][[2]]

}

else{

tbl <- r[[2]][[2]][[7]][[3]][[3]][[10]][[2]]

}

return(tbl)

}

#while( result.num >= 51 ){

while(1){

start.num <- start.num + 1

quote.table <- NULL

quote.url <- paste('http://info.finance.yahoo.co.jp/history/?code=',x,start,end,'&p=',start.num,'&tm=',substr(time.interval,1,1),sep="")

#cat(quote.url)

#try( r <- xmlRoot(htmlTreeParse(quote.url,error=xmlErrorCumulator(immediate=F))), TRUE) # これだと取得時にエラーが出た。。

try(r<-htmlParse(quote.url))

if( is.null(r) ) stop(paste("Can not access :", quote.url))

#try( quote.table <- r[[2]][[1]][[1]][[16]][[1]][[1]][[1]][[4]][[1]][[1]][[1]], TRUE )

#try( quote.table <- extractQuoteTable(r,type), TRUE )

try( quote.table <- xpathApply(r,"//table")[[2]], TRUE )

quote.size<-xmlSize(quote.table)

#cat(paste("size:",quote.size))

if(xmlSize(quote.table)<=1){

return (financial.data)

}

if( is.null(quote.table) ){

if( is.null(financial.data) ){

stop(paste("Can not quote :", x))

}else{

financial.data <- financial.data[order(financial.data$date),]

return(financial.data)

}

size <- xmlSize(quote.table)

for(i in 2:size){

financial.data <- rbind(financial.data,function.financialproduct(quote.table[[i]]))

}

#result.num <- xmlSize(quote.table)

Sys.sleep(1)

}

financial.data <- financial.data[order(financial.data$date),]

return(financial.data)

}

#convert string formart date to POSIXct object

convertToDate <- function(date.string,time.interval)

{

#data format is different between monthly and dialy or weekly

if(any(time.interval==c('d','w'))){

result <- gsub("^([0-9]{4})([^0-9]+)([0-9]{1,2})([^0-9]+)([0-9]{1,2})([^0-9]+)","\\1-\\3-\\5",date.string)

}else if(time.interval=='m'){

result <- gsub("^([0-9]{4})([^0-9]+)([0-9]{1,2})([^0-9]+)","\\1-\\3-01",date.string)

}

return(as.POSIXct(result))

}

#convert string to number.

as.number <- function(string)

{

return(as.double(as.character(gsub("[^0-9.]", "",string))))

}

#return end of month date.

endOfMonth <- function(date.obj)

{

startOfMonth <- as.Date(format(date.obj,"%Y%m01"),"%Y%m%d")

startOfNextMonth <- as.Date(format(startOfMonth+31,"%Y%m01"),"%Y%m%d")

return(startOfNextMonth-1)

}

このコードでYahoo!ジャパンの株価を見てみましょう。ちなみに番号は4689です。どうやら上手く取れているようです。

> quoteStockTsData("4689.t",since="2016-01-01")
         date open high low close   volume adj_close
1  2016-05-02  476  483 475   478 18498100       478
2  2016-04-28  504  508 493   496 11966300       496
3  2016-04-27  505  511 495   497 12973800       497
4  2016-04-26  507  508 495   500  7712600       500
5  2016-04-25  513  515 506   509  7350600       509
6  2016-04-22  515  517 509   514  8908900       514
7  2016-04-21  512  517 506   514 13249900       514
8  2016-04-20  511  515 493   506 14455700       506
9  2016-04-19  516  523 511   516 13345800       516
10 2016-04-18  503  509 499   503 10275900       503
11 2016-04-15  504  519 504   513 16962900       513

> quoteStockTsData("4689.t",since="2016-01-01")

date open high low close volume adj_close

1 2016-05-02 476 483 475 478 18498100 478

2 2016-04-28 504 508 493 496 11966300 496

3 2016-04-27 505 511 495 497 12973800 497

4 2016-04-26 507 508 495 500 7712600 500

5 2016-04-25 513 515 506 509 7350600 509

6 2016-04-22 515 517 509 514 8908900 514

7 2016-04-21 512 517 506 514 13249900 514

8 2016-04-20 511 515 493 506 14455700 506

9 2016-04-19 516 523 511 516 13345800 516

10 2016-04-18 503 509 499 503 10275900 503

11 2016-04-15 504 519 504 513 16962900 513