Buffon在python中的针模拟

Question

import numpy as np
import matplotlib.pylab as plt

class Buffon_needle_problem:

    def __init__(self,x,y,n,m):
        self.x = x #width of the needle
        self.y = y #witdh of the space
        self.r = []#coordinated of the centre of the needle
        self.z = []#measure of the alingment of the needle
        self.n = n#no of throws
        self.m = m#no of simulations
        self.pi_approx = []

    def samples(self):
        # throwing the needles
        for i 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
        for j in range(self.m):
            # n throw
            hits = 0 #setting the succes to 0
            for i in range(self.n):
                #condition for a 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())

对于那些不熟悉布冯问题的人，这里是http://mathworld.wolfram.com/BuffonsNeedleProblem.html

或者

实现相同的想法（和输出） http://pythonfiddle.com/historically-accurate-buffons-needle/

我的预期输出应该是pi的值，但是我的代码给出了大约4。有人能指出逻辑错误吗？在

Answer 1

针对中的取样应为均匀余弦。方法请参见以下链接：http://pdg.lbl.gov/2012/reviews/rpp2012-rev-monte-carlo-techniques.pdf

另外，程序也有一些逻辑问题。这是一个有效的版本。在

#!/bin/python
import numpy as 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 the needle
        self.y = float(y)  #witdh of the space
        self.r = [] #coordinated of the centre of the needle
        self.z = [] #measure of the alignment of the 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 the needles
        for i 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
        for j in range(self.m):
            self.r=[]
            self.z=[]
            self.samples()
            # n throw
            hits = 0 #setting the success to 0
            for i in range(self.n):
                #condition for a 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())

Answer 2

我想说的是，问题是你用一个简单的线性函数来定义针的对中，而实际上，针距中心的有效长度是由正弦函数定义的。在

你想通过使用一个函数来计算针的有效长度（与线成90度角）。在

比如：

self.z.append(np.cos(np.random.uniform(-np.pi/2, np.pi/2))*self.x)

在90°和-90°之间的随机余弦角度。在

作为参考，cos(+/-90) = 0和{}，因此在90°时，针的有效长度为零，而在0°时，其全长。在

我在这台机器上既没有安装mathplotlib也没有安装numpy，所以我不知道这是否可以修复它，但它绝对是必要的。在

Answer 3

看起来你犯了一个简单的舍入错误。下面的代码可以工作，虽然结果不是很接近pi。。。在

import numpy as np
import matplotlib.pylab as plt

class Buffon_needle_problem:
def __init__(self,x,y,n,m):
    self.x = x #width of the needle
    self.y = y #witdh of the space
    self.r = []#coordinated of the centre of the needle
    self.z = []#measure of the alingment of the needle
    self.n = n#no of throws
    self.m = m#no of simulations
    self.pi_approx = []

def samples(self):
    # throwing the needles
    for i 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 simulations
    for j in range(self.m):
        self.r=[]
        self.z=[]
        for i in range(self.n):
            self.r.append(np.random.uniform(0,self.y))
            self.z.append(np.random.uniform(0,self.x/2.0))
        # n throws
        hits = 0 # setting the succes to 0
        for i in range(self.n):
            # condition for a 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.0*(float(self.x)/self.y)*float(self.n)/float(hits)
        self.pi_approx.append(hits)

    return self.pi_approx

y = Buffon_needle_problem(1,2,40000,5)
print (y.simulation())

还请注意，您对所有模拟都使用相同的样本！在