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

乱数を使って確率を計算

• 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%くらい誤差がある。

関連

• 偶然とはなにか。コインを投げてみよう http://qiita.com/reikubonaga/items/12329122772d8161439d
• p値とはなにか http://qiita.com/reikubonaga/items/c4684fcea0810441a955
• A/Bテストは完全ではない http://qiita.com/reikubonaga/items/6f12543aa2478a6ac3b4
• 乱数(モンテカルロ法)を使って確率を計算しよう http://qiita.com/reikubonaga/items/28ee2c2c8074b98b0b04