Calculating $\Pi$ by Monte-Carlo Simulation

-- by Charlie Chengrui Zheng 07/10/2021

As we learned more and more math, we found more and more ways to calculate $\pi$. In computational statistics, there is a way to calculate $\pi$ by brute force -- Monte-Carlo Simulation. In this article, I will do a simple Monte-Carlo Simulation on the calculation of $\pi$, or the area of a circle. This method can also be applied to the calculation of any area of geometric shapes.

First, we are going to build a Cartesian coordinates and set both x-axis and y-aixs from -1 to 1.

Second, draw a circle on it by the formula $x^2 + y^2 = 1$.

Then, we can apply Monte-Carlo Simulation by generating random dots on the coordinates. The dot's x and y are random numbers within -1 and 1. if its $x^2 + y^2 < 1$, then this dot lies inside the circle like the yellow dot (0,0) below. if its $x^2 + y^2 > 1$, then this dot lies outside the circle like the red dot (0,0) below.

import matplotlib.pyplot as plt
import random

fig, ax = plt.subplots()
fig.set_size_inches(5, 5)
#build a Cartesian coordinates and set both x-axis and y-aixs from -1 to 1
ax.set_xlim((-1, 1))
ax.set_ylim((-1, 1))
#draw a circle
circle = plt.Circle((0, 0), 1)
ax.add_artist(plt.Circle((0, 0), 1,fill=False))
#generate the yellow dot (0, 0) that lies inside the circle
ax.plot((0), (0), 'o', color='y')
#generate the red dot (0.9, 0.5) that lies outside the circle
ax.plot((0.9), (0.5), 'o', color='r')

[<matplotlib.lines.Line2D at 0x7ffaf79c4040>]

The idea of calculating $\pi$ is to calculate the probability that a randomly generated dot lies inside the circle.

$$p = (numbe\space of\space dots\space inside\space the\space circle)\space/\space (total\space number\space of\space dots)$$

Since the total area of our coordinates is 4, by multiplying 4 to the probability, we can get the estimated $\pi$.

As for the example above, there is 1 dot inside the circle out of 2 dots in totoal, so $$p = 1/2 = 0.5$$ $$E(\pi) = p*4 = 2$$

This is a very sloppy estimated $\pi$ because the sample size is too small. By scaling up the sample size, we can have more accurate estimated $\pi$.

def simulation(num):
    """
    Monte-Carlo Simulation according to the input sample size and prints out the estimated pi
    
    Args:
        num (int): the sample size
    """
    inside = 0
    fig, ax = plt.subplots()
    fig.set_size_inches(5, 5)
    ax.set_xlim((-1, 1))
    ax.set_ylim((-1, 1))
    ax.add_artist(plt.Circle((0, 0), 1,fill=False))
    #generate num random dots
    for i in range(num):
        x = random.uniform(-1, 1)
        y = random.uniform(-1, 1)
        #calculate whether the dot lies in the circle
        if x**2 + y**2 <1:
            inside += 1
            ax.plot((x), (y), 'o', color='y')
        else:
            ax.plot((x), (y), 'o', color='r')
    #calculate estimated pi
    print(f"pi = {4*inside/num}")

Now, let us try a sample size of 10. In this sample, there are 8 dots inside the circle out of 10 dots. Therefore, p = 0.8 and E($\pi$) = 3.2 Note that, we have an estimated $\pi$ rounding to 1 decimal place because it is divided by 10.

simulation(10)

pi = 3.2

Let us try a sample size of 1000. In this sample, E($\pi$) = 3.184. Note that, we have an estimated $\pi$ rounding to 3 decimal places because it is divided by $10^3$. We can see that the circle is roughly filled up and our estimate is getting closer to $\pi$.

simulation(1000)

pi = 3.184

Let us try a sample size of 100,000. In this sample, E($\pi$) = 3.15468. Note that, we have an estimated $\pi$ rounding to 5 decimal places because it is divided by $10^5$. We can see that the circle is almost filled up and our estimate is getting closer to $\pi$. By passsing in greater sample, we can continue to have more accurate estimated $\pi$.

simulation(100000)

pi = 3.15468