離散的な確率分布をする確率変数の定義に基づいた生成

２項分布は高校でも扱うので、いくらかイメージを持っている。それ以外の超幾何分布やら何やら読んだだけで分かったつもりになっていたが、問題を解こうとすると理解できていないことに気づく。理解の助けになるだろうかと思い、Ｒで、定義に基づいて確率変数を生成してみた。



１　サイコロ
　d=sample(1:6,100,replace=T)

　サイコロを１００回振って出た目をベクトルdに保存。

２　２項分布
 n=2;p=0.5
 d=replicate(100,sum(runif(n)<p))

１回の試行で事象Ａの起こる確率がpである独立試行（ベルヌーイ試行）をn回繰り返す。事象Ａの起こった回数xの確率分布が２項分布Bi(n,p)である。
（０，１）間の一様乱数がpより小さい確率はpであるので、これを持って事象Ａが起きたと考える。runif(n)<pはTRUEとFALSEからなるベクトルだが、これをn回の試行の結果と考え、このうちのTRUEの数を数えればよい。sumはTRUEを１FALSEを０として足す。これを１００回繰り返したときの結果をベクトルdに保存。

３　幾何分布　１回の試行で事象Ａの起こる確率がpである試行を繰り返しx回目にはじめてＡが起きたとする。xの分布を幾何分布という。q=1-pとして、確率関数pq^xが幾何数列（等比数列）の一般項の形をしているところから幾何分布というと「統計学入門」にあった。

 p=0.5
 v=runif(1000)<p

 i=which(v)
 d=i[-1]-i[-length(i)]

通常のプログラムではforやwhileなどを使って起きるまで繰り返しそれまでの回数を数えると考えれば実に素直に実現できそうであるが、Ｒでfor文は使いたくない。で、上のようにした。runif(1000)<pでこの試行１０００回のシミュレーションをする。

TRUE FALSE FALSE FALSE TURE TRUE FALSE FALSE TRUE ・・・
といった列になる。TRUEが起きたこと、FALSEが起きなかったことを表す。
これをいったんvに代入しておく。
whichはTRUEの位置のインデックスからなるベクトルを返す。上の例では
1 5 6 9 ・・・
となる。これを変数iに代入しておく。
TRUEになってから次のTRUEが何回目に起きているかを調べるには隣同士の差を取ればよい。
4 1 3 ・・・
これが求めるものである。隣同士の差を取るために次のようにした。
i[-1]はベクトルiの先頭を削除するので、5 6 9 ・・・となる。
i[-length(i)]は最後を削除して長さをそろえた。
i[-1]からi[-length(i)]を引けばＯＫである。

スクリプト中のいい加減な定数1000は何とかしたいのだが、忍ぶこととした。
iの先頭（上の例では１）を活かすなら、iの先頭に0を付け加えたものを引くことも考えられる。


４　負の２項分布　「統計学入門」では　Bi(1,p)を繰り返し、k回目の成功を得るまでの失敗の回数x、「自然科学の統計学」では　k回目の成功を得るまでに要した回数xとなっていて、定義がちょっと違う。

 p=0.5 ; k=5
 v=runif(100)<p

 i=which(v)
 l=length(i)
 d=i[-(1:k)]-i[-((l-k+1):l)] #「自然科学・・」の定義
 e=d-k #「入門」の定義

幾何分布ではiを１ずらしたが、kずらした。

５　超幾何分布 M個のAとN個のBからn個取り出した時のAの数xの分布

 M=5;N=3;n=2
 S=c(rep(1,M),rep(0,N)) #1がM個、0がN個
 d=replicate(1000, sum( sample(S,n) ) )

意外に簡単にできた。

いくらかイメージがわいてきた。幾何分布するk個の確率変数X1,X2,・・・Xk　の和が負の２項分布するというのも理解できる。ポアソン分布については上のように定義でパシッとはできないかも。近似的な解決になりそうだが今後の課題。