[Python] Pi Estimation Using Monte Carlo Method

 

“So here we have pi squared, which an engineer would call 10.”

Frank King

 

There are several methods for Pi estimation and this one uses Monte Carlo method in doing so. If you need to know more about the theory of what is done here, you can read this post.

What we are going to do in a nutshell, is that we want to produce random numbers, and check whether they will fall inside an imaginary circle with a radius of 1. This circle is inscribed in a square, like this:

So, the numbers we will produce will either fall inside the circle or they will be inside the square (and not in the circle). From simple mathematics, we can calculate Pi number using the following euqation:

So, let’s jump to the code:

"""Econowmics.com"""

#First I begin by importing the modules we need
import random
from math import sqrt

#Setting the seed
random.seed(1991)    #Guess my birth.year

#This function does all the calculations for us
def PiEstimation(n):
    #I define a variable to store the number of points inside the circle
    circle = 0

    #loop to calculate the points we need for our estimation
    for i in range(1,n+1):
        x = random.random()
        y = random.random()
        
        #Checking to see if the produced number falls into the circle
        if sqrt(x**2 + y**2) <= 1:
            circle += 1

    #I calculate Pi for this iteration
    Pi = 4 * (circle/n)

    #returning the result
    return Pi


#n
n=[10,100,1000,10000,100000,10000000,100000000,200000000]

#I define an empty dictionary to store the results of iterations to it
Pi_storage ={}

#Running a loop to estimate Pi with different number of random points
for j in n:
    Pi_storage[j]=PiEstimation(j)

#Printing the results to the console
for k in n:
    print ("The result of pi estimation after " + str(k) + " attempts is equal to: " + str(Pi_storage[k]))

 

Let’s quickly check the code before we jump to the results:

1. First, I have started by importing the modules that are needed for this estimation, and then the seed is set. Setting a set for your random number generation helps to produce the same results whenever the experiment is done.

2. Then a function is created for estimating the Pi number, given number of attempts (n). Inside the function, first a variable is created so that later the number of points inside the circle can be stored to it. Then, and by using a for loop, n instances are created for x and y. Each instance is a random number between 0 and 1. The function then checks if the produced number falls into our circle or not. After all the points have been created, an estimation of Pi is created using the formula described above.

3.  At the next step, I have first created a list containing different values of n that I want to do the Monte Carlo experiment with. Moreover, I create an empty dictionary so that I can store the results of estimation with different values of n inside it. Subsequently, using a for loop and the function created in the previous step, the results of estimation is computed and stored in the dictionary.

4. At the final step, and by using one last loop, the results are printed out to the console in the following format:

 

As it is evident from the results, the larger the number of our attempts (n) is, the closer we get to the value of Pi. However, even after 20 million attempts, just the first two digits of Pi have been estimated correctly.