## Introduction to Montecarlo Method¶

In this notebook, we will calculate value of pi by implementation of Monte Carlo Simulation

In [1]:
import random
import matplotlib.pyplot as plt
import numpy as np
import pylab
import seaborn as sns
sns.set()
%matplotlib inline


#### Lets calculate the value of PI directly:¶

In this method we randomly generate the points and calculate the ratio for value of PI.

In [12]:
n_trials = 10000
n_hits = 0
plt.figure(figsize = [10,10])
for iter in range(n_trials):
x = random.uniform(-1.0, 1.0)
y = random.uniform(-1.0, 1.0)
if x**2 + y**2 < 1.0:
plt.scatter(x,y,color = "blue",marker='.')
n_hits += 1
else:
plt.scatter(x,y,color = "red",marker ='.')
print(4.0 * n_hits / float(n_trials))

3.174

• Lets construct the function to run above code multiple time
In [3]:
def direct_pi(N):
n_hits = 0
for i in range(N):
x = random.uniform(-1.0, 1.0)
y = random.uniform(-1.0, 1.0)
if x ** 2 + y ** 2 < 1.0:
n_hits += 1
computed_pi = 4.0*n_hits / float(n_trials)
return computed_pi

In [8]:
n_runs = 100
n_trials = 1000
plt.figure(figsize = [12,3])
for run in range(n_runs):
pi = direct_pi(n_trials)
plt.scatter(run,pi)
plt.ylim([0,4.0])
plt.xlabel("No of Run")
plt.ylabel("Value of Pi")
plt.show()


#### Marcob chain Calculation of PI:¶

In this method we perform a random walk and accept or reject the move.

In [13]:
x, y = 1.0, 1.0
delta = 0.1
n_trials = 5000
n_hits = 0

plt.figure(figsize = [10,10])
for i in range(n_trials):
del_x = random.uniform(-delta, delta)
del_y = random.uniform(-delta, delta)

'''to make sure they are inside square'''
if (abs(x + del_x) < 1.0 and abs(y + del_y) < 1.0):
x = x + del_x
y = y + del_y

'''to make sure they are inside circle'''
if x**2 + y**2 < 1.0:
n_hits += 1
plt.scatter(x,y,color = "blue", marker='.')
else:
plt.scatter(x,y,color = "red", marker='.')

print(4.0 * n_hits / float(n_trials))

3.2256

• Lets construct the function to run above code multiple time.
In [14]:
def markov_pi(N, delta):
x, y = 1.0, 1.0
n_hits = 0
for i in range(N):
del_x, del_y = random.uniform(-delta, delta),\
random.uniform(-delta, delta)

if abs(x + del_x) < 1.0 and abs(y + del_y) < 1.0:
x, y = x + del_x, y + del_y
if x**2 + y**2 < 1.0: n_hits += 1
computed_pi = 4.0 * n_hits / float(n_trials)
return computed_pi

In [17]:
n_runs = 100
n_trials = 10000
delta = 0.1
plt.figure(figsize = [12,3])
for run in range(n_runs):
pi =  markov_pi(n_trials, delta)
plt.scatter(run,pi)
plt.ylim([0,4.0])
plt.xlabel("No of Run")
plt.ylabel("Value of Pi")
plt.show()