Chapter 14. Example: Estimating Pi

Generally speaking, the Monte Carlo method refers to a technique that attempts to find an approximate solution to a numerical problem by making use of randomly generated inputs. I would like to present as follows an implementation based on such a technique that approximates the value of the constant pi (the ratio of the circumference of circle over its diameter).

The meaning of the following declarations should be clear:

// #define N 1024 // typedef point = $tup(double, double) // typedef points = list0(point) typedef pointss = list0(points) typedef pointsss = list0(pointss) //

Given a grid of the dimension N by N, there are N2 unit squares in it. For each point on the grid, the area of the unit disk centered at the point equals pi (as the radius of the unit disk is 1). Given a unit contained in the grid, the probability of a randomly chosen point on the grid residing within the disk is clearly pi/N2.

Assume that there are K points on the grid. A randomly chosen point generates one hit if it resides within the unit disk centered at one point among the K points. Accumulatively, a randomly chosen point can generate n hits if it resides within the unit disks centered at n points among the K points.

Let create a matrix of the dimension N by N for storing lists of points:

// val theGrid = matrix0_make_elt<points>(N, N, list0_nil) //

We are to play N2 rounds of experiment. At round K for each K between 1 and N2, there are already K-1 points on the grid chosen randomly in the previous rounds. We choose a new point on the grid randomly for round K and records the number of hits it generates (with respect to the K-1 points on the grid). We then add this new point to the list stored at the matrix cell (i, j), where i and j are the integer parts of the x-coordinate and y-coordinate, respectively.

The following functions are for computing the number of hits generated by a point with respect to a collection of points:

// fun hit_point_point ( p1: point, p2: point ) : int = ( if dist(p1, p2) <= 1.0 then 1 else 0 ) // fun hit_point_points ( p0: point, ps: points ) : int = ( list0_foldleft<int><point> ( ps, 0 , lam(res, p1) => res + hit_point_point(p0, p1) ) (* list0_foldleft *) ) // fun hit_point_pointss ( p0: point, pss: pointss ) : int = ( list0_foldleft<int><points> ( pss, 0 , lam(res, ps) => res + hit_point_points(p0, ps) ) (* list0_foldleft *) ) // fun hit_point_pointsss ( p0: point, psss: pointsss ) : int = ( list0_foldleft<int><pointss> ( psss, 0 , lam(res, ps) => res + hit_point_pointss(p0, ps) ) (* list0_foldleft *) ) //

A unit square in the grid is located at (i, j) if the lower-left corner of the unit square is at the point (i, j). Two unit squares are neighbors if they locate at (i1, j1) and (i2, j2) such that abs(i1-i2) and abs(j1-j2) are at most 1, where abs refers to the absolute value function (on integers). Given a point inside the unit square located at (i, j), the point can only generate hits with respect to those residing inside the neighbors of the unit square. Given a unit square at (i, j), the following function generate a collection of the entire neighbors of the unit square:

// fun theGrid_get_neighbors ( i: int, j: int ) : pointsss = let // fun isvalid_row (i: int): bool = (0 <= i && i < N) fun isvalid_col (j: int): bool = (0 <= j && j < N) // fun fopr ( i: int , j: int ) : points = if ( isvalid_row(i) && isvalid_col(j) ) then theGrid[i, j] else list0_nil() // in // int_list0_map<pointss> ( 3 , lam(i') => int_list0_map<points> ( 3 , lam(j') => fopr(i + i' - 1, j + j' - 1) ) ) (* end of [int_list0_map] *) // end // end of [theGrid_get_neighbors] //

It should be clear that the following function performs one round of the aforementioned experiment:

// fun do_one(): int = let // val x = N * randfloat() val y = N * randfloat() // val p = $tup(x, y) // val i = g0float2int_double_int(x) val j = g0float2int_double_int(y) // val nhit = hit_point_pointsss (p, theGrid_get_neighbors(i, j)) // in // theGrid[i, j] := list0_cons(p, theGrid[i, j]); nhit // end // end of [do_one] //

Note that the function randfloat returns a randomly chosen floating point number between 0 and 1 (based on the uniform distribution). The integer returned by a call to do_one records the number of hits generated by the point chosen at one specific round. In order to perform K rounds of experiment, the following function can be called:

// fun do_all ( K: int ) : int = ( fix loop ( i: int, res: int ) : int =<cloref1> if i < K then loop(i+1, do_one()+res) else res )(0, 0) // end of [do_all] //

The integer returned by a call to do_all(K) is the tally of the numbers of hits recorded during rounds n for n ranging from 1 to K.

It can be proven that do_all(N*N)/(N*N) converges to pi/2 (with probability 1) as N approaches the positive infinity. Please find on-line the entirety of the code used in this chapter.