4章 二項分布 大数の法則とシミュレーション


サイコロを6回振っただけでは、1回以上1の目が出る確率は100%にならず66.5%程度でした。

今回は、何回振れば1の目が100%出ると言えるのか、それを試してみたいと思います。

サイコロを振って1の目が出る確率InftyProjectをpとし、その他の目が出る確率InftyProjectをqとします。

n回振って得られる確率は二項分布の公式InftyProjectで求められます。

プログラムでループさせれば簡単に計算できますので、

n回振って1回以上1の目が出る確率を、99回まで試してみます。

 1回目=0.166667

 2回目=0.305556

 3回目=0.421296

 4回目=0.517747

 5回目=0.598122

 6回目=0.665102

 7回目=0.720918

 8回目=0.767432

 9回目=0.806193

10回目=0.838494

11回目=0.865412

12回目=0.887843

13回目=0.906536

14回目=0.922113

15回目=0.935095

16回目=0.945912

17回目=0.954927

18回目=0.962439

19回目=0.968699

20回目=0.973916

21回目=0.978263

22回目=0.981886

23回目=0.984905

24回目=0.987421

25回目=0.989517

26回目=0.991265

27回目=0.992720

28回目=0.993934

29回目=0.994945

30回目=0.995787

31回目=0.996489

32回目=0.997074

33回目=0.997562

34回目=0.997968

35回目=0.998307

36回目=0.998589

37回目=0.998824

38回目=0.999020

39回目=0.999184

40回目=0.999320

41回目=0.999433

42回目=0.999528

43回目=0.999606

44回目=0.999672

45回目=0.999727

46回目=0.999772

47回目=0.999810

48回目=0.999842

49回目=0.999868

50回目=0.999890

51回目=0.999908

52回目=0.999924

53回目=0.999936

54回目=0.999947

55回目=0.999956

56回目=0.999963

57回目=0.999969

58回目=0.999974

59回目=0.999979

60回目=0.999982

61回目=0.999985

62回目=0.999988

63回目=0.999990

64回目=0.999991

65回目=0.999993

66回目=0.999994

67回目=0.999995

68回目=0.999996

69回目=0.999997

70回目=0.999997

71回目=0.999998

72回目=0.999998

73回目=0.999998

74回目=0.999999

75回目=0.999999

76回目=0.999999

77回目=0.999999

78回目=0.999999

79回目=0.999999

80回目=1.000000

81回目=1.000000

82回目=1.000000

83回目=1.000000

84回目=1.000000

85回目=1.000000

86回目=1.000000

87回目=1.000000

88回目=1.000000

89回目=1.000000

90回目=1.000000

91回目=1.000000

92回目=1.000000

93回目=1.000000

94回目=1.000000

95回目=1.000000

96回目=1.000000

97回目=1.000000

98回目=1.000000

99回目=1.000000


この出力を得るのに作ったプログラムはこれです。

##########################################

##二項分布を得る関数 binoD

##   確率分布表 binomial distribution

########################################## 

binoD<-function(n,pp){

 p=pp

 q=(1-p)

 ans=0

 for( i in 0:n){

   ans[i+1]=choose(n,i)*p^i*q^(n-i)

 }

 return(ans)

}

##########################################

#上記関数binoDで1回以上1の目が出る確率を計算

##########################################

n=99

p=1/6

for(i in 1:n){

  a=binoD(i,p)        #関数を呼んでる

  cat(sprintf("%2i回=%f\n",i,sum(a)-a[1]))

}


1行ずつ出力されますが、手作業で上の表に直しています。

この表を見れば、1回以上1の目が出る確率は、17回目で95%を超えて、80回目で100%に達してます。


確率ですので厳密には100%ではありませんが、ほぼ100%とみなすことが出来るます。

80回やって1回以上目が出でる確率を、普通に計算するとInftyProjectですね。


次に進む前に、期待値と平均について触れておきたいと思います。

統計で平均と言っている値は、確率では期待値と呼んでいます。


二項分布にも期待値の公式があり、npで表します。

nは試行回数でpは確率です。

サイコロを3回振った場合、公式を用いると3×InftyProjectとなりInftyProjectが期待値となります。

これを手動で計算すると、p(x=0)+p(x=1)+p(x=2)+p(x=3)を求める事になります。

xは、1の目が出る回数です。3回振るのですから0回~3回の範囲で1の目が出ますね。

p(x=0)なら、 InftyProject

p(x=1)なら、 InftyProject

p(x=2)なら、 InftyProject

p(x=3)なら、 InftyProject

これらを全部足すと0.5になります。計算は同じパターンなので省略しますが6回なら1になります。


二項分布は、成功か失敗かの2つしか結果がありません。

確率に掛けている数字は成功回数そのものになります。

成功とはInftyProjectで起こる事象が生起するという意味で、この場合1の目が出るという事です。

つまり期待値とは理論値と同じになります。

InftyProjectの確率で成功するなら、6回やれば1回は成功する」それが期待値です。


もし分散を求めるのであれば、この値を2乗して下さい。

InftyProjectと計算し、0.666667を求めます。

この値から期待値を2乗した値、つまりInftyProjectを引きます。

InftyProjectが分散値になります。

これも2項分布の分散の公式があり、np(1-p)で表します。

これに当て嵌めると、InftyProjectとなり分散値が求まります。


次はサイコロを沢山振った時の話です。

実際には振れないので、代わりに乱数を用いてシミュレーションしてみます。

image001.png  

# サイコロを30,60,600,6000回振るシミュレーション

a <- c(30,60,600,6000)

g <- par(mfrow=c(2,2))

for (i in 1:4) {

  y <- floor(runif(a[i],1,7))

  s <- paste("n=",floor(a[i]))

  barplot(table(y),main=s) 

}

par(g)


乱数でシミュレーションをしていますので、実行する度に結果は違うはずですが、

nが大きくなるほど、各目の出る差が少なくなっていきます。



今度はサイコロではなくてコインで考えてみましょう。

コインを投げて表を1、裏を0とします。

10回投げて理想通り、表が5回、裏が5回出たとします。

こんな感じに裏表が出て、結果={1,1,1,0,1,0,1,0,0,0}だったとします。

結果の総和はInftyProjectですね。

これだと試行が増えるたびに総和も増えていってしますので、試行回数nで割った式InftyProjectで考えてみると、

InftyProjectとなり、コインの裏表が出る期待値になります。


これも数式よりプログラムでシミュレーションをした方が分り易いと思います。


#コイン投げ5000回のシミュレーション

for(i in 1:5000){

  x[i]<-floor(runif(1,0,2))

  p[i]<-sum(x[1:i])/i

  if( i%%500 == 0 || i<10 || (i>2990 && i<=3000) || (i>4990) ){

     cat(sprintf("P(x%i=%f)\n",i, p[i]))

  }

}


上のプログラムで得た確率変数P(1-5000)をプロットしたものです。

image002.png


途中経過です。

P(x1=0.000000)

P(x2=0.500000)

P(x3=0.666667)

P(x4=0.500000)

P(x5=0.400000)

P(x6=0.333333)

P(x7=0.428571)

P(x8=0.375000)

P(x9=0.333333)

P(x500=0.462000)

P(x1000=0.481000)

P(x1500=0.492000)

P(x2000=0.492000)

P(x2500=0.494800)

P(x2991=0.492477)

P(x2992=0.492647)

P(x2993=0.492482)

P(x2994=0.492318)

P(x2995=0.492487)

P(x2996=0.492323)

P(x2997=0.492492)

P(x2998=0.492662)

P(x2999=0.492497)

P(x3000=0.492333)

P(x3500=0.496571)

P(x4000=0.493500)

P(x4500=0.491778)

P(x4991=0.496293)

P(x4992=0.496394)

P(x4993=0.496495)

P(x4994=0.496596)

P(x4995=0.496697)

P(x4996=0.496597)

P(x4997=0.496498)

P(x4998=0.496399)

P(x4999=0.496499)

P(x5000=0.496400)



試行回数が増えるほどバラツキが小さくなり理論値に近づいているのが分ると思います。


これを大数の法則といいます。

この証明にはチェビシェフの不等式が使われます。

母平均InftyProject、母分散InftyProjectをも確率変数xとその母平均InftyProjectとの差について次の不等式が成り立つというのがそれです。

 InftyProject

チェビシェフの不等式は分布が分らない時でも使える利用範囲が広い式です。

分散が±1の時にこの不等式を利用してみます。

InftyProject

±2,±3,±5,±10でもやってみます。

InftyProject

InftyProject

InftyProject

InftyProject

正規分布なら±1σ=69%、±2σ=96%、±3σ=99.7%ですので、比べるとかなり大雑把な感じになりますね。


さて、大数の法則ですが、チェビシェフの不等式における左辺の確率変数xをn個の互いに独立な確率変数{x1,x2,x3...xn}の平均 Xと想定してみます。

メンドクサイ言い回しですが、x1,x2,x3...xnは観測値と考えて下さい。

サイコロなら1回1回振った結果です。

独立とはサイコロ1回1回の結果は、他に影響を与えないという事です。1回目と2回目の結果は無関係ですよね。

n個の観測値の平均Xの分散はInftyProjectです。

その分散InftyProjectを右辺の母分散に当てはめるとInftyProjectとなるので、nが無限大に大きくなると限りなくゼロに近づいていくのが分ります。つまり期待値と観測値の差が無くなってしまいます。




今回は試行回数と確率の関係について紹介しました。

始めの章で、受け手側と説明した確率ですね。

くれぐれも理論値につられてギャンブルなどしないように気をつけて下さい。


次回は、ベイジアンフィルターを取り上げます。


前章へ 次章へ

メインメニュー