Bending bits...

Bytes and Words...

Monte Carlo algorithms in OCaml

Monte Carlo Integration

module type MonteCarlo_type=sig
  val initializePRNG:unit->unit
  val getRandomNumber:unit->float
  val approximate:(float->float)->float->float->int->int->float->float  
end

let rec approximate f a b n index value=
    if index=(n-1) then
      ((b-.a)/.(float_of_int index))*.value
    else
      let randVal=getRandomNumber() in
      let inInterval=a+.randVal*.(b-.a) in
      let fVal=f inInterval in
      approximate f a b n (index+1) (value+.fVal)

open Random
open MonteCarlo_type

module MonteCarlo:MonteCarlo_type=struct
  let initializePRNG ()=Random.self_init()
  let getRandomNumber ()=Random.float 1.0
  
  let rec approximate f a b n index value=
    if index=(n-1) then
      ((b-.a)/.(float_of_int index))*.value
    else
      let randVal=getRandomNumber() in
      let inInterval=a+.randVal*.(b-.a) in
      let fVal=f inInterval in
      approximate f a b n (index+1) (value+.fVal)
end

open MonteCarlo

let pi = 4.0 *. atan 1.0;;
let func (x:float)=sin x
let () = 
  MonteCarlo.initializePRNG();
  print_float (MonteCarlo.approximate func 0.0 pi 100000 0 0.0)

Monte Carlo PI Estimation

Formula to determine PI:

\begin{equation} PI=4.0*\frac{hits}{darts thrown} \end{equation}
let rnd=System.Random(System.DateTime.Now.Millisecond)

let genRandomNumbers (count:int) =
    List.init count (fun _ -> rnd.NextDouble ())

let isInside (x:double) (y:double)=(sqrt (x*x+y*y))<1.0

let sum (x:^a list)=
    match (isInside (x.Head*2.0-1.0) (x.Tail.Head*2.0-1.0)) with
        | false->0
        | true-> 1

let rec computePi (numThrows:int) (hits:int) (idx:int)=
    if numThrows=idx then
        4.0*((double)hits/(double)numThrows)
    else
          computePi numThrows (hits+(genRandomNumbers 2 |> sum)) (idx+1)

[<EntryPoint>]
let main argv = 
    let l=computePi 1000000000 0 0
    printfn "%F" l
    0