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
