読者です 読者をやめる 読者になる 読者になる

ビッグデータ時代のマーケティングの事例をRでやってみた(1)

R 統計

ビッグデータ時代のマーケティングは名著なのですが、事例に関するコードが公開されていないのが難点です。 逆に一つずつ自分でコードを書くと勉強になりそうなのでひとつやってみました。

2.4.「マーケティングにおけるモデリング例」です。

まず、P46の「データ生成の手順」に従ってデータフレームを作ります。

library(dplyr)
#データ生成
err <- rnorm(100, mean=0, sd=0.3)
Z1 <- runif(100, min=0.5, max=1)
u <- runif(100)
Z2 <- (u <= 0.3)*(Z1 <= 0.8)
y <- exp(0.5 - 2.5 * log(Z1) + 0.8 * Z2 + err)  %>% floor()
data <- data.frame (y, Z1, Z2)

散布図は次のようになり、だいたい、本と近くなっています。 f:id:datascienceconsultant:20161227102812p:plain

モデルが10個あるのでリストを作って格納します。
#モデルリスト
mf <- list(y~1, y~Z1, y~Z1+Z2, y~log(Z1), y~log(Z1)+Z2,
           log(y)~1, log(y)~Z1, log(y)~Z1+Z2, log(y)~log(Z1), log(y)~log(Z1)+Z2)

AICヤコビアンを計算します。

#AIC格納用リスト
AICb <- vector("double", 10)
AICa <- vector("double", 10)
#補正前AIC
for(i in 1:10){
  lm(mf[[i]], data=data) %>%
    AIC() -> AICb[i]
}
#ヤコビアン
jacob <- sum(log(numDeriv::grad(log, x=data$y)))
#補正後AIC
AICa <- AICb + c(rep(0, 5), rep(-2 * jacob, 5))

data.frame(AICb, AICa) %>%
  set_rownames(paste0("model", 1:10)) %>%
  knitr::kable(format="markdown")
AICb AICa
model1 621.08099 621.0810
model2 532.42325 532.4232
model3 502.95307 502.9531
model4 525.10030 525.1003
model5 493.74816 493.7482
model6 264.85676 553.5136
model7 136.97222 425.6290
model8 90.41715 379.0740
model9 139.38691 428.0437
model10 89.20195 377.8588

確かにmodel1~5では、1>2>4>3>5

model6~10では、6>9>7>8>10

model5と10を比べると5>10となり、model10が最適という結果が得られました。

最後に現況再現性の確認をします。

model[[10]] %>% predict() %>% plot(type="l", ylim=c(-1, 3.5))
par(new=T) 
log(y) %>% plot(type="l", lty="dashed", ylim=c(-1, 3.5))
par(new=T) 
(predict(model[[10]]) - log(y)) %>%
  plot(type="l", lwd=2, ylim=c(-1, 3.5))

f:id:datascienceconsultant:20161227104203p:plain

麻美子数

統計 アニメ

ベーコン数をご存知でしょうか。 ケビン・ベーコンという俳優と競演している俳優はベーコン数:1 直接競演はしていないが、ベーコン数1の俳優と競演している俳優はベーコン数:2 という風に数値をつけます。

これの声優版は無いものでしょうか? .lainで調べると出演件数の多い声優さんは

1.能登麻美子 - 533作品

2.子安武人 - 426作品

3.石田彰 - 399作品

4.置鮎龍太郎 - 351作品

5.三木眞一郎 - 341作品

というわけで、能登麻美子さんを起点とした「麻美子数」を定義したい。

The Oracle of Baconのように測定できると良いな……

と思ったら、MamikoNotoさんがIMDbに登録されているので、 The Oracle of Baconそのもので、MamikoNumberの測定ができます!

西明日香:1(モーレツ宇宙海賊映画版)

洲崎綾:2(タマ子ラブストーリー→日高里奈→劇場版とある魔術の禁書目録

大橋彩香:2(映画ドキドキプリキュア→寿美奈子→劇場版とある魔術の禁書目録

巽悠衣子:1(モーレツ宇宙海賊映画版)

うーん。IMDbにはTVアニメは登録されていないっぽい

Rの乱数生成法による違い

R 乱数

blogを始めようと思い立った日に、

R Advent Calendar 2016 - QiitaRの正規乱数生成法(既定)を振り返るという記事が気になりました。

 

Rの乱数生成法には以下の7つがあるのはよく判ったのですが、

Wichmann-Hill, Marsaglia-Multicarry, Super-Duper, Mersenne-Twister,Knuth-TAOCP, Knuth-TAOCP-2002, L'Ecuyer-CMRG

実用上どれくらいの違いが有るのでしょうか?

 

岡田謙介「Rを利用したモンテカルロ法による統計量の評価」(専修大学 心理学研究センター年報 第1号 2012年3月)にある、モンテカルロ法で試してみました。

library(magrittr)
library(dplyr)
library(moments)

#関数の定義
f <- function(x){(cos(50*x)+sin(30*x))^2}
#正解値
corAns <- integrate(f,0,1)$value
#乱数のパターン
rndPtn <- c("Wichmann-Hill", "Marsaglia-Multicarry", "Super-Duper",
"Mersenne-Twister","Knuth-TAOCP", "Knuth-TAOCP-2002", "L'Ecuyer-CMRG")
#結果格納用
result <- matrix(rep(NA, 7*10^3), ncol=7)
result_time <- rep(NA,7)
#メインルーチン
for(j in 1:7){
t<-proc.time()
for(i in 1:(10^3)){
10^3 %>% runif() %>% f() %>% 
sum() %>% divide_by(10^3) -> result[i,j]
}
#実行時間計測
proc.time() - t -> result_time[j]
}
#結果をデータフレームに
result %<>%
as.data.frame() %>%
set_colnames(rndPtn)
#ヒストグラム描画
par(mfrow = c(2, 4))
for(j in 1:7){
hist(result[,j] - corAns, xlab="正解との差", 
main=paste0(rndPtn[j], "\n(time:", round(result_time[j],2), "sec)"))
}
par(mfrow = c(1, 1))

実行時間だとSuper-Duperが頭一つ抜けて速いようです。

#平均,標準偏差,尖度,歪度
colMeans(result) - corAns -> ave
apply(result, 2, sd) -> sde
apply(result, 2, skewness) -> skw
apply(result, 2, kurtosis) -3  -> krt
rbind(ave, sde, skw, krt) %>%
  multiply_by(1000) %>%
  round(1) %>% knitr::kable(format="markdown")
Wichmann-Hill Marsaglia-Multicarry Super-Duper Mersenne-Twister Knuth-TAOCP Knuth-TAOCP-2002 L'Ecuyer-CMRG
ave -0.4 0.6 0.0 2.5 -1.0 1.5 -0.6
sde 35.9 34.3 34.6 34.0 34.5 33.5 35.0
skw -37.3 39.9 -50.5 71.6 42.3 31.7 0.6
krt -138.4 -183.7 -179.2 -127.9 85.5 -98.3 107.2

結果は1000倍しています(歪度は-3した後1000倍)。

試行回数が1000回と少ないからかもしれませんが、これくらいだとMersenne-Twisterのありがたみが薄い、というか他を選びたくなります。

というわけで、目的によってはMT以外を選ぶ意味があるのではないか、と思いました(特に実行時間的な意味で)