dice.hs

import Text.Printf
import Control.Monad
import Control.Monad.MC

-- runGame returns True if Bob rolls a 6 on his second turn
-- Sue cannot roll a 6 or the premise of the problem is invalid.
-- Therefore don't even model her rolls.
runGame :: MC Bool
runGame = do
    turns <- (sequence.repeat) runTurn
    return $ elemIndex True turns == Just 1

runTurn :: MC Bool
runTurn = do
    roll <- sample 6 [1..6]
    return $ roll == 6

runGames :: Int -> MC (Double,Double)
runGames n = do
    xs <- replicateM n runGame
    let m  = length $ filter id xs
        p  = toDouble m / toDouble n
        se = sqrt (p * (1 - p) / toDouble n)
    return (4*p, 4*se)
  where
    toDouble = realToFrac . toInteger

main = let
    n       = 100
    seed    = 0
    (mu,se) = evalMC (runGames n) $ mt19937 seed
    delta   = 2.575*se
    (l,u)   = (mu-delta, mu+delta)
  in do
    printf "Estimate:                  %g\n" mu
    printf "99%% Confidence Interval:   (%g,%g)\n" l u