为什么布丰的针法估算圆周率得不到预期值？

Question

这是我的代码

buffon <- function(a,l)
    {
      n<-0
      N<-0
      repeat {N<-N+1
      print(N)
      p<-c(n/N)
      # Sample the location of the needle's centre.
      x<-runif(1,min = 0,max =a/2)
      print(x)
      # Sample angle of needle with respect to lines.
      theta<-runif(1, 0, pi/2)
      print(theta)
      # Does the needle cross a line?
      k<-l/2*sin(theta)
      ifelse(x<=k,n<-n+1,n<-n)
      p<-c(n/N)
      print(p)
      pie<-(2*l)/(p*a)
      print(pie)
      if(N>5000) {break}
      }
    }

我正在尝试使用布冯针的想法来估计 pi 的值，但是，当我尝试时：buffon(2,3)，最终估计是 3.8，这远大于3.1。如果我的代码中有任何错误或者我不能使用 pi 来估计 pi，有人可以向我解释吗？

加法：

我意识到我的代码有很多行是多余的，所以我今天早上修改了一下：

  buffon01 <- function(n,a,l)
  {
    # Sample the location of the needle's centre.
    x<-runif(n,min = 0,max =a/2)
    # Sample angle of needle with respect to lines.
    theta<-runif(n, 0, pi/2)
    # Does the needle cross a line?
    k<-l/2*sin(theta)
    # l is the length of the needle
    # a is the distance between to parallel line
    v<-length(x[x<=k])
    p<-c(v/n)
    pie<-(2*l)/(p*a)
    list("pi"=pie,"error"=abs(pie-pi))
  }

通过设置大于 l 的方式，我可以获得非常接近 3.14 的结果...但是我得到的结果非常不稳定...就像我再次进行相同的实验一样 3.1* 可能是除了 4 之外的任何数字...我是否忽略了设置中的其他一些问题？

Answer 1

据我所知。你的问题是你使用的是 "long" 针（其中 l > a），然后公式 l/2*cos(theta) 不起作用。在那种情况下，您需要使用稍微 more elaborate formula.

所以我稍微清理了你的代码，并修改了确保 l < a:

buffon <- function(a,l) {
  stopifnot(l < a)  
  n<-0
  N<-0
  repeat {
    N <- N+1

    # Sample the location of the needle's centre.
    x <- runif(1, min=0, max=a/2)

    # Sample angle of needle with respect to lines.
    theta <- runif(1, min=0, max=pi/2)

    # Does the needle cross a line?
    n <- ifelse(x <= l/2*cos(theta), n + 1, n)

    # Estimate probability of crossing line
    p <- n/N

    # Compute pi
    pie <- (2*l)/(p*a)

    if (N>50000) { # Note the increased iterations
      break
    }
  }
  return(pie)
}

ans <- buffon(a=3,l=2)
print(ans)
#[1] 3.159621

看来您还把上述公式中的 cos 换成了 sin，以检查针是否碰到一条线。然而——在我的脑海中——这并不重要。最后，在 "repeat" 循环内打印使得函数在每次遇到它时都打印出来（从而填满控制台）。相反，我已将函数修改为 return pie 对象，因此您可以将其存储在变量中。

希望这对您有所帮助。

编辑：

您的新简明函数说明了问题。比较：

buffon01(1e6, a=3, l=2)
#$`pi`
#[1] 3.143023
#
#$error
#[1] 0.00143062

buffon01(1e6, a=2, l=3)
#`pi`
#[1] 3.850819

#$error
#[1] 0.7092266

在这里，您看到使用 l>a 失败，因为使用了错误的公式。关于收敛到 pi，Buffon 的针收敛速度很慢，需要多次投掷才能获得合适的估计值。但这是统计问题，here。

为什么布丰的针法估算圆周率得不到预期值？

Why Buffon's needle method in estimating pi doesn't produce the expected value?

pi

r