シミュレーションの仕方

「楽しい　統計学セミナー　- 証拠を科学する方法 -」の「第1講」 において，シミュレーションプログラムが掲示されているが

system.time({
set.seed(12345) # 乱数の初期値を設定しておこう
n <- 10000      # プログラムの中で何回も出てくる数値は変数に代入しておいて使おう（一箇所で変更できる）

# for (i in 1:n) { # 等差数列を作るというのだが...
#  x[i] <- 6*i
# }
x <- 1:n*6         # ベクトル演算はこんなときにも使える (1:n)*6 と書くかは，好みの問題

y <- numeric(n)    # １～n回の試行を格納するための変数を確保．

for (i in 1:n) {
  a <-sample(x=1:6, size=x[i], replace=TRUE, prob=c(rep(1/6, 6)))
# y[i] <- (sum(a[a == 1])) / x[i] # 足してデータの個数で割るというのは...
  y[i] <- mean(a == 1)            # 平均値を求めるということ
}                                 # ちなみに，(a[a == 1]) は冗長

plot(y, type="l")
abline(a = 1 / 6, b = 0, col="red")
})

このプログラムでは，「6 回サイコロを振る」ということを，毎回 1～x[i] 繰り返す。
もう一つの考え方は，1 ~ n 回繰り返し，その途中経過を見るということでもよい。もっとも，こちらの考え方では k 回の時点での結果は k-1 回までの結果に従属であるということではあるが（前述の方法では，各回の結果は独立）。

もう一つのポイントは，1の目が出るということに注目するのであるから，2～6の目が出ることをシミュレーションする必要はないということ。つまり 1 とそれ以外の目が 1/6, 5/6 で出るというシミュレーションをすればよいということ。TRUE/FALSE や，成功回数を数えるときに sum が直接使えるということを考えれば，1 の目がでることを 1，それ以外の目が出ることを 0 で表すとよい。つまり n 回の試行結果は x <- sample(0:1, n, replace=TRUE, prob=c(5/6, 1/6)) で得られ，k 回目までの成功数は cumsum(x) で得られ，成功率は cumsum(x)/1:n で得られることになる。ということで，纏めると

system.time({
set.seed(12345)
n <- 10000
x <- matrix(sample(0:1, 6*n, replace=TRUE, prob=c(5/6, 1/6)), 6)
y <- cumsum(colSums(x))/6/1:n
plot(y, type="l")
abline(a = 1 / 6, b = 0, col="red")
})

前のプログラムの実行時間は
   ユーザ   システム       経過  
    15.418      2.327     17.636

後ろのプログラムの実行時間は
   ユーザ   システム       経過  
     0.136      0.004      0.155