Edited at

乱数(モンテカルロ法)を使って確率を計算しよう

More than 3 years have passed since last update.


乱数を使って確率を計算


  • 0から1の乱数をn個作成 runif(n, 0, 1)

  • 0.5%未満の数だけ取得 which(x<rate)

  • 0.5%未満の数を数える

  • 上の施行を1万回繰り返して、分布を調べる

  • quantileは、下位x%番目の数字を出している

  • => コイン投げの実験結果とほぼ等しい

> n=100;rate=0.5;z<-c();for(i in 1:10000){x=runif(n, 0, 1);y=which(x<rate);z[i]=length(y)/n};summary(z)

Min. 1st Qu. Median Mean 3rd Qu. Max.
0.3000 0.4700 0.5000 0.4996 0.5300 0.6900
> quantile(z, c(.01,.05,.25,.5,.75,.95,.99))
1% 5% 25% 50% 75% 95% 99%
0.39 0.42 0.47 0.50 0.53 0.58 0.61
> n=1000;rate=0.5;z<-c();for(i in 1:10000){x=runif(n, 0, 1);y=which(x<rate);z[i]=length(y)/n};quantile(z, c(.01,.05,.25,.5,.75,.95,.99))
1% 5% 25% 50% 75% 95% 99%
0.463 0.474 0.489 0.500 0.511 0.526 0.536


施行数と精度の関係


  • 施行数を重ねると、だんだん精度が良くなる。

  • 施行数を100倍にすると10倍の精度。
    => 詳しく知りたい人は、「大数の定理」「一様分布」など調べてみて

n=10000;rate=0.05;z<-c();for(i in 1:10000){x=runif(n, 0, 1);y=which(x<rate);z[i]=length(y)/n};quantile(z, c(.01,.05,.25,.5,.75,.95,.99))

> n=100000;rate=0.03;z<-c();for(i in 1:1000){x=runif(n, 0, 1);y=which(x<rate);z[i]=length(y)/n};quantile(z, c(.01,.05,.25,.5,.75,.95,.99))
1% 5% 25% 50% 75% 95% 99%
0.0289197 0.0291700 0.0296400 0.0299700 0.0303525 0.0309700 0.0312300
> n=100000;rate=0.03;z<-c();for(i in 1:4000){x=runif(n, 0, 1);y=which(x<rate);z[i]=length(y)/n};quantile(z, c(.01,.05,.25,.5,.75,.95,.99))
1% 5% 25% 50% 75% 95% 99%
0.02879 0.02911 0.02962 0.02998 0.03036 0.03088 0.03128
=> 1000回と4000回だと4倍施行数が違うので、精度が2倍になる。その差を見て、誤差だと思えるなら、その施行数で問題ない。95%だと0.0309700と0.03088で0.3%くらい誤差がある。


関連