Created
February 11, 2015 01:00
-
-
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.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
""" | |
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