Simulating Buffon's Needle in R [closed]

Question

I have the following question I'm trying to write a simulation for:

Let T_1 be the number of crossings in n tosses of the needle, then $$E_1 = T_1d/(nl)$$ is an unbiased estimator of 2/π. Write a program to simulate E_1 using n = 100,000 needle tosses.

I have found, hopefully correctly, the minimum l in terms d as the problem asks (not shown here). My question, however, is how do I choose T_1 to get the correct estimate of 2/π? My code below seems to work, but the answer is incorrect in terms of decimal places dependent on choice of T_1. The question doesn't give any additional comments regarding T_1 either. I feel like I'm missing something.

Est <- function(n, T, d){
  E <- c()
  l <- (d*pi)/4
  for(i in 1:n) {
  E[i] <- (T*d)/(n*l)
  }
  mean <- mean(E)
  c.int <- quantile(E, c(0.025, 0.975))
  return(c(mean, c.int))
}

Additionally, how do I also get the function to return c.int?

**Updated code from OP **

Est <- function(n, d){ 
  E <- c() 
  l <- (d*pi)/4 
  for(i in 1:n) {
    x <- runif(n, 0, pi/2)
    y <- runif(n, 0, d/2)
    ncross <- y <= d/2 * sin(x)
#    E[i] <- (ncross*d)/(n*l) 
    E[i] <- mean(ncross) * d/l
  } 
  mean <- mean(E)
  c.int <- quantile(E, c(0.025, 0.975))
  return(c(mean, c.int))
}

Est(10000, 1)

Answer 1

This is how I approached the problem, based on the Wikipedia article and your code:

## Let el be the length of the needle and te be the distance between two lines.
## https://en.wikipedia.org/wiki/Buffon%27s_needle_problem

Est <- function(n, el, te) {
  theta <- runif(n, 0, pi/2)
  x <- runif(n, 0, te/2)
  # cross is a vector of true/false. Take the mean to find a proportion
  cross <- x <= el/2 * sin(theta)

  return(mean(cross))
}

I then set up some more code to call the function a few times so I could get some repeated samples:

el = 1
te = 5
E = c()

for(i in 1:1000) {
# only if el <= te
  E[i] = 2*el/(Est(1e5, el, te)*te)
}

mean(E)
quantile(E,  c(0.025, 0.975))
hist(E)

And this was my result:

> mean(E)
[1] 3.141163
> quantile(E,  c(0.025, 0.975))
    2.5%    97.5% 
3.090951 3.189290 

Added analysis: So if $T_1$ is the number of crossings in $n$ trials, we have the probability a needle crosses is $ P = \frac{T_1}{n} $ and from the Wikipedia article, $ P = \frac{2l}{d\pi} $, where $d$ is the distance between lines (aka te in the code) and l is the length of the needle (aka el in the code). The code calculates $E = \frac{2l}{d\cdot P} \approx\pi $ or $ E=\frac{2l\cdot n}{d\cdot T_1} \approx \pi \Rightarrow \frac{2}{\pi} \approx \frac {T_1\cdot d} {n \cdot l}$

Tldr: change the calculation of E to

        # only if el <= te
  E[i] = (Est(1e5, el, te)*te)/el

which is equivalent to $\frac {T_1\cdot d} {n \cdot l}$

Answer 2

One major problem with your code is that involves the use of $\pi$, and so it cannot be regarded as a genuine implementation of Buffon's needle to estimate $\pi$. A genuine implementation of the needle algorithm would generate the position and direction of the needles, and determine whether a needle crosses a line, without ever using $\pi$ (or any trigonometric functions) in the algorithm. Generation of a series of indicators from needles would allow you to estimate $\pi$, with an appropriate confidence interval, using the central limit theorem.

---

Buffon's needle algorithm without using $\boldsymbol{\pi}$: Buffon's needle experiment can be implemented using a rejection-sampling method that does not require the use of $\pi$. This can be done over any appropriately sized lined space. For simplicity, we will consider the simplest case of a space that is a unit square $\mathcal{S} = [0,1]^2$, where the left and right boundaries are the "lines" in the experiment. Let $\mathbf{M} \sim \text{U}(\mathcal{S})$ denote the midpoint of the thrown needle, which is uniformly distributed on the unit square. (Note that this means that the needle may lie partly off the unit square, and may cross the left or right boundary lines.)

To determine the direction of the needle, use the following rejection-sampling method. Generate a proposed value $\mathbf{D} \sim \text{U} (\mathcal{S})$ and accept this value if $||\mathbf{D}|| \leqslant 1$. This gives a value that is uniformly distributed in the unit circle, so the value $\mathbf{D}/||\mathbf{D}||$ is uniformly distributed on the boundary of the unit circle. We take the needle to be aligned in the direction of this vector. For simplicity, we consider the case of a "short" needle, with length $\ell \leqslant 1$. Since this needle has half-length $\ell/2$, the two end-points of the needle are:

$$\mathbf{E}_1 \equiv \mathbf{M} + \frac{\ell}{2} \cdot \frac{\mathbf{D}}{||\mathbf{D}||} \quad \quad \quad \mathbf{E}_2 \equiv \mathbf{M} - \frac{\ell}{2} \cdot \frac{\mathbf{D}}{||\mathbf{D}||}.$$

The needle crosses one of the left or right boundary lines if and only if the horizontal coordinate of one of these vectors falls outside the unit interval. (They cannot both fall outside this interval, since the midpoint of the needle is in the unit square.) The indicator for crossing one of the boundary lines is denoted as $H$, and can be written as:

$$H \equiv \mathbb{I}(E_{1,1} < 0) + \mathbb{I}(E_{1,1} > 0) + \mathbb{I}(E_{2,1} < 0) + \mathbb{I}(E_{2,1} > 0).$$

This gives an algorithm to generate a single indicator value for the needle crossing the boundary. Note that this algorithm did not involve the use of $\pi$, since the direction of the needle was determined using a simple rejection-sampling method that only involved the generation of uniform random variables. It can be shown that $H \sim \text{Bern}(2 \ell / \pi)$, so we can estimate $\pi$ by generating a large number of indicator values $H_1,...,H_n$ and taking $\hat{\pi} = 2 \ell / \bar{H}$ where $\bar{H}$ is the sample mean of the indicators. Since you want to estimate $2/\pi$, you would use the point estimate $2/\hat{\pi} = \bar{H} / \ell$.

---

Implementation in R: We will use the above algorithm to generate a vector $H_1,...,H_n$ for $n$ needles. We construct a function where we specify a needle length l and the number of needles n. The function generates a vector of indicator values for these needles, indicating whether they crossed the boundary lines.

BUFFON_NEEDLES <- function(l, n = 1) {

#Check inputs
if (!is.numeric(n))     { stop('Error: Number of needles must be numeric') }
if (length(n) != 1)     { stop('Error: Number of needles should be a scalar') }
if (n != as.integer(n)) { stop('Error: Number of needles must be an integer') }
if (n < 1)              { stop('Error: Number of needles must be at least one') }
if (!is.numeric(l))     { stop('Error: Needle length must be numeric') }
if (length(l) != 1)     { stop('Error: Needle length should be a scalar') }
if (l < 0)              { stop('Error: Needle length must be positive') }
if (l > 1)              { stop('Error: Needle length cannot be greater than one') }

#Set output vector
H <- rep(0, n);

#Generate values for output vector
for (i in 1:n) {

  #Generate midpoint of needle
  M  <- runif(2);

  #Generate direction of needle
  D  <- c(1,1);
  while (norm(D, type = '2') > 1) { D <- runif(2); }
  DD <- D/norm(D, type = '2');

  #Determine endpoints of needle
  E1 <- M + (l/2)*DD;
  E2 <- M - (l/2)*DD;

  #Determine whether needle crosses lines
  H[i] <- (E1[1] < 0)|(E1[1] > 1)|(E2[1] < 0)|(E2[1] > 1); }

#Give output
H; }

We can implement this function for a large value of n and use this to get a point estimate of $\pi$ as follows.

#Set number of needles and needle length
l <- 0.4;
n <- 10^7;

#Generate indicators
set.seed(1);
HHH <- BUFFON_NEEDLES(l, n);

#Estimate pi
2*l/mean(HHH);
[1] 3.142286

#Estimate 2/pi
mean(HHH)/l;
[1] 0.6364792

As you can see, taking $n = 10^7$ is sufficient to get quite a good estimate of $\pi$. Since you propose to use $n = 10^5$ your estimator will have less accuracy, but it should still give you a rough estimate of the true value.