28人中にB型が3名以下になる確率

OOPL


Rに苦戦中にて候 - ときどきの雑記帖 i戦士篇 を読んで R の練習にちょうど良さそうだな…とチャレンジしてみようと思ったのですが、その前にとりあえず Squeak Smalltalk を使ってオーソドックスに。

| 血液型群 Ｂ型が三名以下 試行回数 bag |
血液型群 := RunArray runs: #(4 3 2 1) values: #(A O B AB).
Ｂ型が三名以下 := 0.
試行回数 := 100000.
bag := (1 to: 28) asBag.
試行回数 timesRepeat: [
    ((bag collect: [:each | 血液型群 atRandom]) occurrencesOf: #B) <= 3
        ifTrue: [Ｂ型が三名以下 := Ｂ型が三名以下 + 1]].
Ｂ型が三名以下 / 試行回数 * 100.0   "=> 16.062 % "


あくまで候補のＢ型の割合がおかしなことになっていないという前提ですが。^^;


直接求める方法はないものかと、Wikipedia を眺めていたらこんなのを見つけました。

全住民の5%がある感染症に罹患しており、その中から無作為に500人を抽出する。このとき、抽出された集団の中に罹患者が30人以上いる確率はどれくらいか。

二項分布 - Wikipedia


これを適当に置き換えて、

28人を無作為に抽出したとき、B型の人が3名以下になる確率を求めよ。


同じく Squeak Smalltalk でこの確率を計算してみると約16％と似たような値が出たのであっているっぽいですね。

(3 to: 0 by: -1) inject: 0 into: [:sum :each |
    sum + ((28 take: each) * (0.2 raisedTo: each) * (1 - 0.2 raisedTo: 28 - each) asFloat)]
"=> 0.1601826710989386 "


では、R でもシミュレーションを…と思ったのですが時間切れなので後ほど。


追記：
…ということで R で書いてみました。新しく runif、table、factor を覚えました。

xs <- array(runif(28))
p2b <- function(x) if (x>=0.6) "A" else if (x>=0.3) "O" else if (x>=0.1) "B" else "AB"
bs <- function(x) table(factor(mapply(p2b,runif(28)),levels=c("A","O","B","AB")))["B"]
n <- 10000
table(mapply(bs, 1:n) <= 3)["TRUE"] / n  #=> 0.1642


ついでに Smalltalk のを Ruby に書き換えた版も。

btypes = [4,3,2,1].zip(["A","O","B","AB"]).collect{ |n,t| [t]*n }.flatten
count = 0
n = 10000
n.times{
  count += 1 if ((1..28).collect{ |e| btypes.sample }).count("B") <= 3
}
count.to_f / n   #=> 0.1577