As one advances in physics, statistics start to become a very useful and intuitive tool that many students will appreciate, especially during their Statistical Mechanics course. So today, we will be exploring a very interesting, and simple for that matter, application of statistics to help us estimate the value of Pi.
The Math For this method, we will be imagining a simple scenario. Imagine for a moment, that we have the unit circle, inside a square.
By having the unit circle, we immediately figure out that the area of the square will be four since the radius of the circle is defined at one, which means that our square will have sides with a value of two. Now here’s where things get interesting, if we take the ratio of both of areas, we end up getting the following
Both of these geometric figures end up having a ratio of pi over four between them, which is an important value for our next step in which we use a bit of imagination. For a moment, image that you have on the ground a circle inside a square, and that you started randomly throwing some sticks into the circle/square diagram. Some of your sticks will most likely land inside the circle and a couple will likely land inside the square but outside the circle. Using this concept is how we will code our pi estimator, by throwing some random numbers into the unit circle equation, as seen below
Furthermore, taking a ratio of throws that landed inside our circle and the total number of throws, we can then formulate the following
And by combining our ratio between the unit circle with the square with this new equation, we can assume the next equation
With this equation, we can finally start coding our estimator and see how close we can get to pi’s actual value.
The code For this code, we will be importing three packages
import random #To use random numbers import numpy as np #For arrays import matplotlib.pyplot as plt #For plotting From here, we start initializing our variables that we will be using later on
trials = 1000 #How many times we are going to "throw" our values inside_the_unit_circle = 0.0 #How many land inside the circle, starting with zero attempts = 1 #We will be starting at one with attempts since our total attempts can't be zero since we can't divide by zero
pi_estimates = np.array([]) #Array to store our pi values for graphing
Now, we just start a while loop that will continue going until we reach the amount of trials we’ve set our code to do
while(attempts <= trials): random_throws = random.random(), random.random() #Setting two random values #By using the unit circle equation, we will be throwing x and y values inside the equation and if they are less than 1, that is, contained within the unit circle, then we will be counting that as a successful throw and we will be adding that estimation to our array if (random_throws[0]**2 + random_throws[1]*_2 < 1): inside_the_unit_circle = inside_the_unit_circle + 1 pi_values = 4_inside_the_unit_circle/attempts pi_estimates = np.append(pi_estimates, pi_values) attempts += 1 Finally, we just go ahead and write a few more lines of codes to make graphing easier
axis_length = len(pi_estimates) #Tells us how big our pi values array is x_axis = np.linspace(0,axis_length,axis_length) #Creating a x-axis by using numpy's linspace. The first parameter tells us where to start, the second parameter tells us where to end, and the third parameter tells us how many points to make. By default, linspace uses 50 for the third parameter, meaning that if we have a range from zero to ten, it will split those ten entries by 50, turning into a 1/5 step size. Since I want a step size of 1, I use the same parameter that I used in the second, for the third parameter. plt.plot(x_axis, pi_estimates) #Plotting the values plt.axhline(y=np.pi,color='r') #Creating a y-axis line, using pi as the value, and colored red All combined, without annotations, our code looks like this
import random import numpy as np import matplotlib.pyplot as plt trials = 1000 inside_the_unit_circle = 0.0 attempts = 1 pi_estimates = np.array([]) while(attempts <= trials): random_throws = random.random(), random.random() if (random_throws[0]**2 + random_throws[1]*_2 < 1): inside_the_unit_circle = inside_the_unit_circle + 1 pi_values = 4_inside_the_unit_circle/attempts pi_estimates = np.append(pi_estimates, pi_values) attempts += 1
axis_length = len(pi_estimates) x_axis = np.linspace(0,axis_length,axis_length) plt.plot(x_axis, pi_estimates) plt.axhline(y=np.pi,color='r') An important note! When we run this code, we will be using random numbers. Once we re-run the same code, we will have a different output since Python will assign new random numbers to our random_throws value. If you want to receive the same output each time you run this code, all you need to do is write
random.seed(0) You can place this line along with the variables along the beginning and each time you re-run the code, you should get the same output as before.
Overview Once we run this code, I end up getting the following graph
It’s quite interesting to see our estimation start with low accuracy but as we increase our attempts, we start to get a convergence on the the value of pi, as shown by our red line. Statistical methods like this, and other more complex versions, are nice tools to understand and experiment within the world of physics. I highly recommend taking a look into some Statistical Mechanic concepts to see the beauty behind the application of statistics and probability in physics, and maybe take some time play with these concepts in Python!