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

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