Monte Carlo Pi

Using Monte Carlo methods, we can create a simulation that approximates pi. In this post, we will go over the math behind the approximation and the code.

The Math

Pi is a mathematical constant consisting of the ratio between the circumfrence of a circle and it’s diameter.

The circumfrence of the circle is defined to be $$ C = 2\pi r$$ while the diameter of the circle is $$d = 2r$$

Take the ratio between the two and you get $$\frac{2\pi r}{2r} = \pi$$

Now let us consider the area of a circle. One can derive the area of a circle by taking the integral of the circumfrence with respect to it’s radius $$ A_{circle} = \int{(2\pi r) dr} = \pi r^2 $$

Let us simplify the formula more by setting the radius equal to one. $$A_{circle} = \pi$$

Now consider only the first quadrant of the circle. Since our circle is centered at the origin and all the points on the circumfrence is equidistant from the center, our area is now $$A_{circle} = \frac{1}{4} \pi$$

And bound the quarter-circle in a 1×1 box with an area of $$A_{square} = 1^2 = 1$$

Notice that the ratio between the circle and the square is a quarter of pi $$\frac{A_{circle}}{A_{square}} = \frac{\frac{1}{4} \pi}{1} = \frac{1}{4} \pi$$

Simulation and Statisitcs

The formula for a circle centered at the origin with radius one is $$x^2 + y^2 = 1$$

Let us focus again on the first quadrent, and do a Monte Carlo simulation to find the area of the quarter-circle

We can do this by what is called the dart board method. We generate a random x and y between 0 and 1. If it satisfies the inequality $$x^2 + y^2 \leq 1$$ then it counts as being inside the circle, if not then it lies outside the circle.

That point will count as an really small area. The point will always be inside the square but may sometimes be inside the circle. Running the simulations a large number of times allows us to add up all the tiny little areas that make up the circle and the square.

To add up these small areas we need to make an assumption. The assumption is that the variance of all the little Monte Carlo trials are the same. Since we are using a psuedo-random number generator, it is safe to assume it is met.

This will allow us to perform a pooled empiricle probability on the simulations to sum up the areas.

Meaning the area of the circle will be the number of times that the inequality was satisfied $$A_{circle} = # Successes$$

And the area of the square will be the number of times the simulation was run, since the random numbers generated will always be between 0 and 1 $A_{square} = # Trials$

Recall that taking the ratio of the area of the circle and the area of the square is a fourth of pi. $$\frac{\frac{1}{4} \pi}{1} = \frac{1}{4} \pi$$

Multiply this number by 4 and you get the value for pi.

This tells us that four times the probability that the randomly generated point is in the circle is equal to pi.

$$\pi = 4 * (Probability\ of\ being\ inside\ circle) = 4 * \frac{# Success}{# Trials} = 4 * \frac{A_{circle}}{A_{square}}$$

Implementation

For the Monte Carlo simulation I used Java. The BigDecimal implementation was used so that there wouldn’t be any issue with integer size limits

/** Calculates Pi
  * @author Brandon Rozek
*/
// Big Integers are used so we don't run into the integer size limit
import java.math.BigInteger;
import java.math.BigDecimal;
class MonteCarloPi {
        public static void main(String[] args) {

                BigInteger successes = BigInteger.ZERO;
                BigInteger trials = BigInteger.ZERO;

For this simulation, we will run 1,000,000,000 trials

BigInteger numTrials = new BigInteger("1000000000");
/*
    Monte Carlo Simulation
        Generate a random point 0 <= x < 1 and 0 <= y < 1
        If the generated point satisfies x^2 + x^2 < 1
           Count as a success
        Keep track of the number of trials and successes
*/
for (; trials.compareTo(numTrials) < 0; trials = trials.add(BigInteger.ONE)) {
    double randomX = Math.random();
    double randomY = Math.random();
    if (Math.pow(randomX, 2) + Math.pow(randomY, 2) < 1) {
        successes = successes.add(BigInteger.ONE);
    }
}

And then we finalize it with a quick calculation of pi

// (Number of successes) / (Number of trials) * 4 gives the approximation for pi
BigDecimal pi = new BigDecimal(successes)
                       .divide(new BigDecimal(trials))
                       .multiply(new BigDecimal("4"));
System.out.println("The calculated value of pi is: " + pi);
}}

Conclusion

We found an approximation of pi using the Monte Carlo methods! I find that really awesome, however, there are some concerns I have with this approach.

1. We don’t keep track of double counting. One possible solution for this is increasing the radius and bounding box appropriately so that the probability of double counting is low.

2. Speed. The more trials you ask it to run, the longer it takes to perform all of the simulations. One possible way around this is to write a parrallel version of this code. That’s possible because of the equal variance that we spoke of earlier. Pooling the successses and trials will still result in a good approximation.