By Evgenia "Jenny" Nitishinskaya and Delaney Granizo-Mackenzie
Notebook released under the Creative Commons Attribution 4.0 License.
A stochastic way to estimate the value of $\pi$ is to sample points from a square area. Some of the points will fall within the area of a circle as defined by $x^2 + y^2 = 1$, we count what percentage all points fall within this area, which allows us to estimate the area of the circle and therefore $\pi$.
In [34]:
# Import libraries
import math
import numpy as np
import matplotlib.pyplot as plt
In [45]:
in_circle = 0
outside_circle = 0
n = 10 ** 4
# Draw many random points
X = np.random.rand(n)
Y = np.random.rand(n)
for i in range(n):
if X[i]**2 + Y[i]**2 > 1:
outside_circle += 1
else:
in_circle += 1
area_of_quarter_circle = float(in_circle)/(in_circle + outside_circle)
pi_estimate = area_of_circle = area_of_quarter_circle * 4
In [46]:
pi_estimate
Out[46]:
We can visualize the process to see how it works.
In [47]:
# Plot a circle for reference
circle1=plt.Circle((0,0),1,color='r', fill=False, lw=2)
fig = plt.gcf()
fig.gca().add_artist(circle1)
# Set the axis limits so the circle doesn't look skewed
plt.xlim((0, 1.8))
plt.ylim((0, 1.2))
plt.scatter(X, Y)
Out[47]:
Finally, let's see how our estimate gets better as we increase $n$. We'll do this by computing the estimate for $\pi$ at each step and plotting that estimate to see how it converges.
In [39]:
in_circle = 0
outside_circle = 0
n = 10 ** 3
# Draw many random points
X = np.random.rand(n)
Y = np.random.rand(n)
# Make a new array
pi = np.ndarray(n)
for i in range(n):
if X[i]**2 + Y[i]**2 > 1:
outside_circle += 1
else:
in_circle += 1
area_of_quarter_circle = float(in_circle)/(in_circle + outside_circle)
pi_estimate = area_of_circle = area_of_quarter_circle * 4
pi[i] = pi_estimate
plt.plot(range(n), pi)
plt.xlabel('n')
plt.ylabel('pi estimate')
plt.plot(range(n), [math.pi] * n)
Out[39]: