大数の法則のシミュレーション

『統計学入門』（東京大学出版会）の中で大数の法則を説明している章があり、その中でコンピュータシミュレーションによる実験が紹介されている。

成功確率p=0.4のベルヌーイ試行を20000回行い,100回おきに成功率を計算して、nが増えるに従いp=0.4に近づくかを確かめている。

本にはこの実験を4回行った結果が紹介され、その内3回は成功率が0.4に収束していく様子が図示されており、１回については理論的な成功率である0.4に収束していない様子が示されている。

さて、この20000回の試行で経験的確率が理論的確率に収束しないことはどの程度の頻度で起きるのだろうか。本にはそこまで書いていないので実際に試してみるしかあるまい。

Rでこのシミュレーションを再現する。とりあえず次のようにコードを書いて同じ実験をできるようにする。

n <- 20000　　　   #試行回数
aa <- numeric(length=n)
bb <- numeric(length=n)
cc <- numeric(length=n/100)

for (i in 1:n){
  aa[i] <- runif(1,min = 0, max = 1) 　#0から1の一様乱数を発生
  bb[i] <- ifelse(aa[i]<0.4,1,0)　　　 #0.4より小さい場合に成功(1)
  if(i%%100 == 0) {                    #100回に1回成功率を計算
    cc[i/100] <- sum(bb[1:i])/i
    }
}

plot(cc,type="l",ylim=c(0.35,0.45))
abline(h=0.4,col="blue")

一回試してみたのが次のもの、この場合だとやや0.4からは離れていますね。

何回か試して見ると、0.4に近かったり離れていたりする。

どの程度の確率でずれるのかを見るために、20000回の試行を行って得た成功率を記録し、同じシミュレーションを1000回繰り返すことにする。（つまり乱数は都合2千万回発生することになる）

for (j in 1:){　
　#ここにさっきの関数

　dd[j] <-cc[length(cc)]　 #一番最後の試行を記録
}

そして得た1000回分の成功率をヒストグラムにしたものが次のもの。

0.01よりも大きな差はないようだが、0.005ぐらいの差は結構な頻度で起きるようだ。（0.395を下回るのは82回、0.405を上回るのは63回あった）元の本の図を見る限りだと、0.005ぐらいの差で収束していないと判断しているようなので、その基準に従うなら、今回のシミュレーションでも15％弱は経験的確率が理論的な確率に収束していないと判断できそうである。