Example: Estimating the Constant Pi

I present as follows a Monte Carlo approach to estimating the constant Pi, the ratio of the circumference of a circle over its diameter.

Assume that we have a square of the dimension N by N, where N is a relatively large natural number (e.g., 1000), and a disk of radius 1 that is contained in the square. Let N2 stand for N*N, that is, the square of N. If we randomly choose a point inside the square, then the probability for the point to hit the disk is Pi/N2.

The experiment we use to estimate the constant Pi can be described as follows. Given a natural number K, let us randomly choose K points inside the square in K rounds. In each round, we choose exactly one point. If the point chosen in round k hits on the disk centered at a previously chosen point, then we record one hit. Clearly, the expected number of hits recorded in round k is (k-1)*Pi/N2 as k-1 points have already being chosen in the previous rounds. Therefore, in K rounds, the expected total number of hits is (K*(K-1)/2)*Pi/N2. If K is fixed to be N2, then the expected total number of hits is (N2-1)*Pi/2. It can be proven that the total number of hits divided by N2 converges to Pi/2 (with probability 1) as N approaches infinity.

If we implement the above experiment directly based on the given description, the time-complexity of the implementation is evidently proportional to N2*N2 as the time spent in round k is proportional to k, where k ranges from 1 to N2. An implementation as such is simply impractical for handling N around the order 1000 (and thus N2 around the order of 1,000,000). To address the issue, we can impose a grid on the square, dividing it into N2 unit squares (of the dimension 1 by 1). We then associate with each unit square a list of chosen points that are inside it. In each round, we first choose a point randomly inside the original square; we next locate the unit square that contains this point; we then only search the lists associated with the unit square or any of its neighbors to count the number of hits generated by the point chosen in this round as this point cannot hit any disks centered at points that are not on these lists. As each unit square can have at most 8 neighbors and the average length of the list associated with each square is less than 1 during the experiment, the time spent during each round is O(1), that is, bounded by a constant. Hence, the time taken by the entire experiment is O(N2).

An implementation that precisely matches the above description plus some testing code is available on-line.