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

user_8307869 0
Viewed 598

Here's my code

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}
      }
    }

I am trying to estimate the value of pi using the idea of Buffon's needle, however, when I try:buffon(2,3), the final estimation is 3.8, which is far greater than 3.1. Could someone explain to me if there's any mistakes in my code or I cannot use pi to estimate pi?

Addition:

I realized that many lines of my code are redundant so I modified it a bit this morning:

  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))
  }

By setting a way larger than l I am able to get a fairly close result to 3.14... but the result I get is very unstable... as in if I do the same experiment again 3.1* could be any number other than 4... Did I ignore some other problems in my setting?

pi r
1 Answers

As far as I can tell. Your problem is that you are using a "long" needle (where l > a), and then the formula l/2*cos(theta) does not work. In that case, you need to use the slightly more elaborate formula.

So I have cleaned your code a bit, and modified ensure that 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

I appears you have also switched cos with sin in the stated formula for checking if the needle hits a line. However---off the top of my head---that does not matter. Lastly, printing inside the "repeat"-loop makes the function print every time it is encountered (and thus filling the console). Instead, I have modified the function to return the pie object, so you can store it in a variable.

Hope this helps.

Edit:

Your new concise function illustrates the problem. Compare:

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

Here, you see that using l>a fails because the wrong formula is used. In regards to convergence to pi, Buffon's needle is quite slow to converge and a lot of throws are needed to get a decent estimate. But that is statistical question, which is touched upon here.

user_15683060
Related