ダメ出し：コンパイラなんかに頼らない

「TokyoR25」にて，library(compiler) の威力を見せつけようと？？しているようだが，compile したときとしないときの比較ではなくて，単に標準誤差はサンプルサイズの平方根に反比例するという統計学上の当たり前の話をなぞっているだけのような。

ここでも，もう何回となく取り上げてきて，うんざりということだが，プログラムを書く一人一人についてはうんざりでもなんでもなく，え～そんなことがあるの？ということなんだろうねえ。

元のプログラムは，シミュレーションごとに標本を生成している。それは，ダメだっっていってるでしょ。まとめて生成するの。必然的に for も使わなくて済むから。

> library(ggplot2)
> library(compiler)
>
> age <- data.frame(age = c(23, 25, 27, 27, 29, 26, 31, 31, 28, 29, 23, 26, 29,
+     30, 30, 27, 26))
>
> getMean <- cmpfun(function(times, s = 10) {
+     sample_mean <- numeric(times)
+     for (i in 1:times) {
+         res <- sample(age$age, size = s, replace = TRUE)
+         sample_mean[i] <- mean(res)
+     }
+     return(sample_mean)
+ })

だからね，13秒もかかっちゃうわけよ。

> system.time(res1 <- getMean(1e+06))
   ユーザ   システム       経過  
    13.384      0.143     13.650

あいかわらず，ggplot の図も「きたない」し。

> ggplot(data.frame(res1), aes(x = res1)) + geom_histogram(binwidth = 0.1, colour = "white") +
+     geom_vline(x = mean(res1), color = "red") + annotate(x = 26, y = 40000,
+     geom = "text", label = paste(sep = ":", "mean", round(mean(res1), 1))) +
+     annotate(x = 26, y = 30000, geom = "text", label = paste(sep = ":", "sd",
+         round(sd(res1), 3))) + xlim(24, 32)

まとめて標本を生成するとはこういうこと。

> getMean2 <- cmpfun(function(times, s = 10) {
+     res <- matrix(sample(age$age, size = s*times, replace = TRUE), times) # はい，ここですよ！
+     return(rowMeans(res)) # シミュレーション結果を返すのも一発で
+ })

で，こうなるわけさ。

> system.time(res2 <- getMean2(1e+06))
   ユーザ   システム       経過  
     0.231      0.057      0.287

> ggplot(data.frame(res2), aes(x = res2)) + geom_histogram(binwidth = 0.1, colour = "white") +
+     geom_vline(x = mean(res2), color = "red") + annotate(x = 26, y = 40000,
+     geom = "text", label = paste(sep = ":", "mean", round(mean(res2), 1))) +
+     annotate(x = 26, y = 30000, geom = "text", label = paste(sep = ":", "sd",
+         round(sd(res2), 3))) + xlim(24, 32)

定石なんだから，ああだこうだいわずに，いつでも採用すべき方策。