Calculating PI using Monte Carlo Simulation

Monte Carlo Simulation to Compute PI

Equations

Code

Follows a sample OCaml code to implement the above idea.

let one_point () = (Random.float (2.) -. 1., Random.float (2.) -. 1.)

let n_points n =
  let rec iter acc = function
      0 -> acc
    | n -> iter ((one_point ()) :: acc) (n - 1)
  in
  iter [] n

(* tail-recursive implementation of map to prevent stack-overflow. *)
let mymap f l =
  let rec iter f acc = function
      [] -> acc
    | h :: t -> iter f ((f h) :: acc) t
  in
  iter f [] l

(* tail-recursive implementation of filter to prevent stack-overflow. *)
let myfilter f l =
  let rec iter f acc = function
      [] -> acc
    | h :: t -> if (f h) then (iter f ((f h) :: acc) t)
                else (iter f acc t)
  in
  iter f [] l

let pi =
  let total = 10000000 in                                         (* number of times the dart is thrown *)
  let points = n_points total in                                  (* list of points all the darts hit *)
  let rs = mymap (fun (x, y) -> x *. x +. y *. y) points in       (* the distance from origin of all points *)
  let inner_points = myfilter (fun r -> r <= 1.) rs in            (* list of all points inside or on the unit circle *)
  let num_inner_points = List.length inner_points in              (* number of points inside or on the unit circle *)
  4. *. ((float_of_int num_inner_points) /. (float_of_int total)) (* value of PI. *)

OCaml Interpretor Interaction

# #use "pi.ml";;
val one_point : unit -> float * float = <fun>
val n_points : int -> (float * float) list = <fun>
val mymap : ('a -> 'b) -> 'a list -> 'b list = <fun>
val myfilter : ('a -> bool) -> 'a list -> bool list = <fun>
val pi : float = 3.1419352

Notes

A large number of points are generated within the bounding square. The probability distribution for the points over this area is uniform, which means that the probability of occurrence of a dot at all points in the given square is equal. With large values of N, the ratio of the number of points inside the unit circle to N starts approaching the ratio of the areas of the unit circle and the square. This idea is used to compute a close approximate to Pi.

As we use OCaml's library module Random to generate the points, the value of Pi computed in every simulation is non-deterministic and is different from the other by a small amount. The time needed to do the above computation was about 10 seconds on my laptop.

Most of the lines in the above code are really there to either enhance the readability or scalability of the program. As List.map and List.filter are non-tail recursive implementations, they lead to stack overflow for very large values of total (>= 1000000). The tail-recursive implementations of map and filter are provided only to prevent stack overflow for very large values of total (>= 1000000).