module Main where

import Data.Array (Array, accumArray, elems, Ix, inRange)
import System.Random (RandomGen, newStdGen, random)
import Control.Monad.State (State, state, evalState)

histogram :: Ix i => (i, i) -> [i] -> Array i Int
histogram bounds bins = accumArray (+) 0 bounds [(x, 1) | x <- bins, inRange bounds x]

bernoulli :: RandomGen gen => Double -> State gen Int
bernoulli p =
  do x <- state random
     return $ if x < p then 1 else 0

binomial :: RandomGen gen => Int -> Double -> State gen Int
binomial count p =
  do ns <- sequence $ replicate count $ bernoulli p
     return $ sum ns

simulate :: RandomGen gen => gen -> Int -> State gen a -> [a]
simulate gen tries m = evalState (sequence $ replicate tries m) gen

test :: Double -> Int -> Int -> IO ()
test bias nails tries =
  do gen <- newStdGen
     let bins = simulate gen tries $ binomial nails bias
     let h = histogram (0, nails) bins
     print $ elems h

main :: IO ()
main = test 0.5 100 10000