title: "Lecture 15: MCM Simulation" author: "Junaid Hasan" date: "August 2, 2021"¶
In [2]:
import math
import random
import matplotlib.pyplot as plt
%matplotlib inline
If this does not work !pip install matplotlib or similar
In [88]:
random.random?
Estimate pi using circle and square¶
We want to find the probability that a point lands inside a circle inscribed in a unit square.
In [89]:
N = 1000
insidecircle = 0
listofpi = []
listofstep = []
estimatepi = 0.0
for i in range(1, N+1):
x = random.random()
y = random.random()
if ((x-0.5)**2) + ((y-0.5)**2) < 0.25 :
insidecircle = insidecircle +1
estimatepi = 4.0*insidecircle/i
if i%100 :
listofpi.append(estimatepi)
listofstep.append(i)
In [90]:
approximatePi = 4.0*insidecircle/N
In [91]:
approximatePi
Out[91]:
3.108
In [92]:
plt.plot(listofstep, listofpi)
plt.show()
Buffon Needle simulation¶
In [93]:
math.pi?
In [5]:
N = 10000
hits = 0
currentestimate = 1
listofestimates = []
listofsteps = []
needle = []
for i in range(1, N+1):
x1 = random.random()
y1 = 0
angle = random.random()*2*math.pi
x2 = x1 + 1.0*math.cos(angle)
y2 = y1 + 1.0*math.sin(angle)
if x2 > 1.0 or x2 < 0 :
needle.append(1)
hits = hits +1
currentestimate = 1.0*hits/i
else:
needle.append(0)
if i% (N/100) == 0 :
#print("Current estimate", currentestimate)
listofestimates.append(currentestimate)
listofsteps.append(i)
In [95]:
plt.plot(listofsteps, listofestimates)
Out[95]:
[<matplotlib.lines.Line2D at 0x7fe467a18eb0>]
In [84]:
listofestimates[-1]
Out[84]:
0.6361990859725791
In [85]:
2/math.pi
Out[85]:
0.6366197723675814
In [6]:
plt.hist(needle)
Out[6]:
(array([3548., 0., 0., 0., 0., 0., 0., 0., 0.,
6452.]),
array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ]),
<BarContainer object of 10 artists>)
In [ ]: