組み合わせの数4 ーモンテカルロシミュレーション版

ちょうしにのってモンテカルロシミュレーション

ある事象が確率pで起きるとき
n回のうち、k回の事象が起きる確率は
{}_n C_k p^k (1-p)^{n-k}
となる。

そこで、「n回のうちk回成功」という確率をシミュレーションで求めて
その後にコンビネーションを逆算する。

シミュレーションの試行回数は10万回とします。

combi n k = do let t = 100000 --試行回数
               let p = 0.5 --成功確率
               nums <- sequence(map sequence (replicate t (replicate n randomFIO)))
               putStrLn $ show $ (fromIntegral(length(filter (== k) $
         map (length.filter (> p)) nums)))/ fromIntegral(t) / (sub p n k)
            where
               sub p n k= p ^ k *  (1-p)^(n-k)
               randomFIO :: IO Float
               randomFIO = randomIO

モナドを理解してないので、すっきりしないコードです。
(副作用ある上に返り値はモナドだし)

*Main> combi 7 3
34.98752
*Main> combi 7 3
35.26272
*Main> combi 12 4
493.52704
*Main> combi 6 1
5.91104
*Main> c 7 3
35
*Main> c 12 4
495
*Main> c 6 1
6

シミュレーションなので誤差があります。
数が大きくなればなるほど目立つのかな。
ちなみに計算に結構時間がかかりますので注意。