¶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