Skip to content

Instantly share code, notes, and snippets.

@napsternxg
Created February 11, 2015 01:00
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save napsternxg/9d8fae71247e93765cbc to your computer and use it in GitHub Desktop.
Save napsternxg/9d8fae71247e93765cbc to your computer and use it in GitHub Desktop.
Monte Carlo Method for approximating the value of pi using ratio of area of circle inscribed in unit length square.
"""
Using Monte Carlo to find the value of pi.
Intiution: Ratio of area of circle to area of uniq length squre it is inscribed it can be used to approximate pi.
Area of circle = pi/4
Area of square = 1
Ratio = pi/4
Let us generate random points x,y inside the square. The probabilty of the point being inside circle is equal to the above ratio.
Circle centered on origin.
Given: -0.5 <= x,y <= 0.5
A point is inside circle if x**2 + y**2 <= 1/4
"""
import matplotlib.pyplot as plt
from matplotlib import gridspec
import numpy as np
plt.clf()
MAXN = 10000 # Max number of iterations
n = 1 # Current value of iteration
n_arr = [] # Saving values of iterations
pi_arr = [] # Approx values pi for given iteration
n_circle = 0 # Number of points inside circle
fig = plt.figure(figsize=(10,15))
gs = gridspec.GridSpec(2, 1, height_ratios=[2, 1])
ax = [plt.subplot(gs[0]),plt.subplot(gs[1])]
# Create a circle just for representation
circle = plt.Circle((0, 0), radius=0.5, fc='y')
ax[0].add_patch(circle)
# Run the simulation
while n <= MAXN:
p_color = 'b' # If point is outside circle mark it with blue color
x = np.random.uniform(-0.5,0.5)
y = np.random.uniform(-0.5,0.5)
if (x**2 + y**2) <= 0.25:
n_circle += 1
p_color = 'r' # If point is outside circle mark it with blue color
ax[0].plot(x,y,p_color+'+')
n_arr.append(n)
pi_arr.append(n_circle*4.0/n) # Value of pi = 4*points in circle/points overall
n+= 1
print n, n_circle, pi_arr[-1], n_arr[-1]
ax[1].plot(n_arr,pi_arr,"b-",label="monte carlo")
ax[1].plot(n_arr, np.pi*np.ones(len(n_arr)),"r-", label="actual")
ax[1].legend()
ax[1].set_xlabel("Number of iterations")
ax[1].set_ylabel("Value of pi")
plt.title("Monte Carlo approximation of pi")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment