Buffon's needle simulation in python -
import numpy np import matplotlib.pylab plt class buffon_needle_problem: def __init__(self,x,y,n,m): self.x = x #width of needle self.y = y #witdh of space self.r = []#coordinated of centre of needle self.z = []#measure of alingment of needle self.n = n#no of throws self.m = m#no of simulations self.pi_approx = [] def samples(self): # throwing needles in range(self.n): self.r.append(np.random.uniform(0,self.y)) self.z.append(np.random.uniform(0,self.x/2.0)) return [self.r,self.z] def simulation(self): self.samples() # m simulation j in range(self.m): # n throw hits = 0 #setting succes 0 in range(self.n): #condition hit if self.r[i]+self.z[i]>=self.y or self.r[i]-self.z[i] <= 0.0: hits += 1 else: continue hits = 2*(self.x/self.y)*float(self.n/hits) self.pi_approx.append(hits) return self.pi_approx y = buffon_needle_problem(1,2,40000,5) print (y.simulation()) for unfamiliar buffon's problem, here http://mathworld.wolfram.com/buffonsneedleproblem.html
or
implementing same idea (and output) http://pythonfiddle.com/historically-accurate-buffons-needle/
my expected output should value of pi code give me around 4. can point out logical error?
the sampling of needle's alignment should uniform cosine. see following link method: http://pdg.lbl.gov/2012/reviews/rpp2012-rev-monte-carlo-techniques.pdf
also, there few logical problems program. here working version.
#!/bin/python import numpy np def sample_cosine(): rr=2. while rr > 1.: u1=np.random.uniform(0,1.) u2=np.random.uniform(0,1.) v1=2*u1-1. rr=v1*v1+u2*u2 cc=(v1*v1-u2*u2)/rr return cc class buffon_needle_problem: def __init__(self,x,y,n,m): self.x = float(x) #width of needle self.y = float(y) #witdh of space self.r = [] #coordinated of centre of needle self.z = [] #measure of alignment of needle self.n = n #no of throws self.m = m #no of simulations self.p = self.x/self.y self.pi_approx = [] def samples(self): # throwing needles in range(self.n): self.r.append(np.random.uniform(0,self.y)) c=sample_cosine() self.z.append(c*self.x/2.) return [self.r,self.z] def simulation(self): # m simulation j in range(self.m): self.r=[] self.z=[] self.samples() # n throw hits = 0 #setting success 0 in range(self.n): #condition hit if self.r[i]+self.z[i]>=self.y or self.r[i]-self.z[i]<0.: hits += 1 else: continue est =self.p*float(self.n)/float(hits) self.pi_approx.append(est) return self.pi_approx y = buffon_needle_problem(1,2,80000,5) print (y.simulation()) 
Comments
Post a Comment