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

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

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

まず、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)

散布図は次のようになり、だいたい、本と近くなっています。

モデルが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))

麻美子数

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

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

1.能登麻美子 - 533作品

2.子安武人 - 426作品

3.石田彰 - 399作品

4.置鮎龍太郎 - 351作品

5.三木眞一郎 - 341作品

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

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

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

西明日香：１（モーレツ宇宙海賊映画版）

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

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

巽悠衣子：１（モーレツ宇宙海賊映画版）

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

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

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

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

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以外を選ぶ意味があるのではないか、と思いました（特に実行時間的な意味で）