Skip to main content

Estimating Pi with Python

π is one of the fundamental constants in mathematics. It is an irrational number and it's value has been calculated up to an extremely large number of digits. In most cases a value of 3.14 will suffice as an approximation but there are many ways in which it can be approximated to a higher degree of accuracy. A particularly fun and simple example can be completed with Python in approximately 15 lines of code.

For this little project I am going to use numpy and I begin by importing this Python module as np,

import numpy as np

If we considering a quarter circle with unit radius we can use N uniform randomly generated points between 0 and 1 to estimate it's area. We can generate our circle for plotting later using the following code,

x = np.linspace(0, 1, 100)
y = np.sqrt(1 - x**2)

Our N randomly generated points are produced like so,

r_x = []
r_y = []
for i in range(N):
    r_x.append(np.random.uniform(0, 1))
    r_y.append(np.random.uniform(0, 1))

Two examples of this are shown below for N = 100 and N = 10000 randomly generated points. The orange line shows the parameter of  our quarter circle.


By summing up the number of points inside the quarter circle we can approximate it's area and it is clear that increasing the value of N improves the accuracy of our estimate. We perform this sum by the following,

counts_inside = 0
counts_outside = 0
for i in range(len(r_x)):
    distance = r_x[i]**2 + r_y[i]**2
    if distance <= 1:
        counts_inside += 1
    else:
        counts_outside +=1

Where distance delineates the radius of a circle drawn through each random point centered on the origin. If this distance exceeds 1 then the point is outside our quarter circle.

Now we note that the area of our quarter circle is given by π/4 since our radius is 1. Consequently we can evaluate π using the following formula,

$\frac{A_{quarter~~circle}}{A_{unit~~square}} = \frac{counts\_inside}{counts\_inside + counts\_outside}$

which becomes,

$\pi = 4 \frac{counts\_inside}{counts\_inside + counts\_outside}$

and in Python this translates too,

pi_estimate = 4*(counts_inside/(counts_outside+counts_inside))

As N increases the estimate of the area of both the unit square and the quarter circle become more accurate and consequently so does the estimate of $\pi$. This is shown in the left hand figure below where the blue line tracks the estimate of $\pi$ as N increases. The black line is numpy's stored value of $\pi$ which is commonly used. The right hand figure shows the ratio of the numpy value and the estimate calculated using the quarter circle which tends to 1 with increasing N. For $N = 10^7$ this method approximates $\pi$ to be 3.1414856 which to 4 decimal places matches the accepted value!
The full code including that used to produce the figures can be found at https://github.com/htjb/Blog/tree/master/Estimating_Pi_with_Python.
Labels:

Comments

Post a Comment

Popular posts from this blog

Contour Plotting with Matplotlib

Visualising data well is an important part of any analysis and a good handle on the Python package Matplotlib is essential for any Python data analyst. I hope to provide a few tutorials on some of the more complex concepts in data visualisation and how to produce clear and tidy graphs with Matplotlib. I will assume some basic knowledge of the Matplotlib package but will try and explain the code as clearly as possible. Comments are always welcome! I am going to begin with this piece on contour plotting which is an area I have a particular interest in. Specifically I am interested in plotting parameter spaces for fitted functions with contours defined by an objective function like $\chi^2$. We will get to an example of this shortly but I first want to look at how we make contour plots with a simpler example. Basic Example with Radii For our simple example we will define variables $x$ and $y$ over a given range and plot the corresponding radius, $Z$ from 0 for each data point $(x, y)$ as ...
Post a Comment
Read more

Random Number Generation: Box-Muller Transform

Knowing how to generate random numbers is a key tool in any data scientists tool box. They appear in multiple different optimisation routines and machine learning processes. However, we often use random number generators built into programming languages without thinking about what is happening below the surface. For example in Python if I want to generate a random number uniformally distributed between 0 and 1 all I need to do is import numpy and use the np.random.uniform() function. Similarly if I want gaussian random numbers to for example simulate random noise in an experiment all I need to do is use np.random.normal(). But what is actually happening when I call these functions? and how do I go about generating random numbers from scratch? This is the first of hopefully a number of blog posts on the subject of random numbers and generating random numbers. There are multiple different methods that can be used in order to do this such as the inverse probability transform method and I ...
Post a Comment
Read more

Binary Search with Python

The Binary Search is an important algorithm for any Data Analyst to know. Frequently, we may need to identify if a particular value can be found in a data set and if it can we may need to know it's index. For example we may have some conditional task that only gets acted upon if our data set features a specific value.  Perhaps, another example, we have collected some data but we know having plotted it in some informative way that there is a spurious value in the set that isn't following the expected trend. We can most likely attribute it to a typo, 349 instead of 345 say, when entering our data but either way we consider it an outlier and we need to identify it's position in the list and remove it from the data set. A possible solution is to loop through our data set until we find the entry that equals 349 but our data set could be big and this could be very time consuming. The alternative is to use the Binary Search which allows us to find the index of our typo in a sorted...
Post a Comment
Read more