No.13062 Re: 多項の超幾何分布について 【青木繁伸】 2010/07/15(Thu) 11:00
単なる二項分布？
問題を言い換えると（シミュレーションの手順でもありますが）以下のようになりますか？
出
現比率 (0.1, 0.2, 0.3, 0.1, 0.16, 0.14) の事象を考える。それぞれの事象を数字 1 〜 6 であらわす。1 〜 6
の数字を出現比率（母比率）で含む，例えば 30000 個の数字の母集団を考える。母集団から 3 個の数字を取り出して，その数字が 1, 2, 3
であれば正解。求めたいのは，3 つの数字のうち 1, 2, 3 のいずれかである個数。
つまり，1 ~ 6 の数字のうち 1, 2, 3 の数字になる確率。1，2，3 いずれかの数字になる確率は 0.1+0.2+0.3 = 0.6。つまり，試行数 3，母比率 0.6 の二項分布。> n <- 1600000 # シミュレーション試行回数
> p <- c(0.1, 0.2, 0.3, 0.1, 0.16, 0.14) # 母比率
> d <- matrix(sample(6, n*3, prob=p, replace=TRUE), n) # 母集団作製。n行3列。行が 1 シミュレーション（3 つの数字の組）
> head(d) # 最初の方の様子
[,1] [,2] [,3]
[1,] 3 3 1 # 正解が 3 つ
[2,] 2 3 2 # 正解が 3 つ
[3,] 2 3 1 # 正解が 3 つ
[4,] 3 5 5 # 正解が 1 つ
[5,] 4 3 3 # 正解が 2 つ
[6,] 3 4 3 # 正解が 2 つ
： 以下略
> table(d)/(n*3) # n*3 個の数字の集計をして，n*3 で割って母比率を確認
d
1 2 3 4 5 6
0.1000371 0.1999710 0.3000129 0.1000596 0.1599306 0.1399887 # 十分な母集団
> # 何をやるかを明らかにするため，あえて for を使う
> tbl <- numeric(4) # 正解の個数の集計に使う
> for (i in 1:n) { # n 回のシミュレーションの集計
+ s <- 0 # 3 つの数字のうち 3 より小さいもの（正解）の数を数える
+ for (j in 1:3) { # 3 つの数字について
+ if (d[i,j] <= 3) s <- s+1 # 3 以下（正解）ならカウント
+ }
+ tbl[s+1] <- tbl[s+1]+1 # 正解の数によって，頻度を蓄積（集計）
+ }
> print(tbl) # 集計表を出力
[1] 101990 460538 692853 344619 # 順に，正解が 0, 1, 2, 3 となった回数
> # 本当は以下の一行でできる
> print(table(apply(d, 1, function(x) sum(x <= 3))))
0 1 2 3
101990 460538 692853 344619 # 上の集計と同じになる（あたりまえ）
> # で。二項分布ならどうなるか。
> # 二項分布
> p <- dbinom(0:3, 3, p=0.6)
> p
[1] 0.064 0.288 0.432 0.216 # 理論比
> p*n
[1] 102400 460800 691200 345600 # 理論度数
No.13073 Re: 多項の超幾何分布について 【ひの】 2010/07/15(Thu) 18:34
先の考え方で計算してみました。処理系は仮称10進BASIC（http://hp.vector.co.jp/authors/VA008683/）です。DIM X(1 to 6) ! 事象Xi の確率を保持する配列
DIM N(0 to 3) ! 正解数i個の確率を集計していく配列
LET X(1) = 0.1
LET X(2) = 0.2
LET X(3) = 0.3
LET X(4) = 0.1
LET X(5) = 0.16
LET X(6) = 0.14
LET N(0) = 0
LET N(1) = 0
LET N(2) = 0
LET N(3) = 0
for i = 1 to 6
for j = 1 to 6
for k = 1 to 6
if i<>j and j<>k and k<>i then ! 同じものが重複するケースを除外
! 確率の計算
LET P = X(i)*X(j)*X(k) / (1 - X(i)) / (1 - X(i) - X(j))
! 正解数(A)のカウント
let A = 0
if i <=3 then let A = A + 1
if j <=3 then let A = A + 1
if k <=3 then let A = A + 1
! A個の正解数の確率をを加算
let N(A) = N(A) + P
end IF
next k
next j
next i
! 各正解数の確率を表示
for i = 0 to 3
print using$("----%",i), using$("---%.##########", N(i))
next i
! 念のため確率の合計を表示
print "Total" ; using$("----------------------%.##########" , N(0)+N(1)+N(2)+N(3))
END
====================出力結果=========================
0 0.0211994474
1 0.3519786887
2 0.5490440861
3 0.0777777778
Total 1.0000000000
13073k and k