ビッグデータ時代のマーケティングの事例をRでやってみた(1)
ビッグデータ時代のマーケティングは名著なのですが、事例に関するコードが公開されていないのが難点です。 逆に一つずつ自分でコードを書くと勉強になりそうなのでひとつやってみました。
まず、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格納用リスト 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))
麻美子数
ベーコン数をご存知でしょうか。 ケビン・ベーコンという俳優と競演している俳優はベーコン数: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の測定ができます!
洲崎綾:2(タマ子ラブストーリー→日高里奈→劇場版とある魔術の禁書目録)
大橋彩香:2(映画ドキドキプリキュア→寿美奈子→劇場版とある魔術の禁書目録)
うーん。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以外を選ぶ意味があるのではないか、と思いました(特に実行時間的な意味で)