As with the Buffon needle problem, the relevant stochastic behaviour here depends only on the proportionate size of the needle and the grid relative to one another. Therefore, without loss of generality, we can stipulate the side-length $a=1$ and let $0 < \ell \leqslant 1$ be the relative length of the needle compared to the length of one side of the grid. Since we are simulating the outcome using a computer algorithm, we can consider just one grid square, and drop the needle so that its midpoint is equally likely to land at any point in this square.
As in the case of Buffon's needle problem, since the goal is to obtain a Monte Carlo estimate of $\pi$, it is appropriate to develop an algorithm that does not use $\pi$ as an input, nor any trigonometric functions or other derivative objects. The Buffon-Laplace needle experiment can be implemented by simulating the midpoint of the needle uniformly on the grid-square, and simulating the direction of the needle uniformly over the unit circle. Since we need to do this without using $\pi$, we will use a rejection-sampling method.
There are a number of possible Monte Carlo estimators for $\pi$ that can be used in the Buffon-Laplace experiment. In our algorithm, we will use the estimator recommended in Perlman and Wichura (1975). This estimator determines the proportion $\bar{H}_n$ of needles that cross one or more of the gridlines, and estimates $\pi$ using the equation:
$$\frac{1}{\hat{\pi}_n} = \frac{(4 \ell - \ell^2)}{\bar{H}_n}.$$
Perlman and Wichura show that this estimator is the UMVUE of $1/\pi$, and is asymptotically efficient. The estimator is optimal when we use a needle that is the same length as a side of the grid (i.e., when $\ell=1$ in our notation). In this case, they show that this estimator provides $12.08$ times as much statistical information as the corresponding estimator in the Buffon needle experiment with only one set of lines. Thus, amazingly enough, the addition of the second set of lines in the grid is equivalent to having more than twelve times as much data!
Buffon-Laplace needle algorithm without using $\boldsymbol{\pi}$: We define the unit square $\mathcal{S} = [0,1]^2$ to represent a single square of the grid. Each of the boundaries of this square are gridlines 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 a gridline.
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.$^\dagger$ Since the 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 a gridline if and only if one or both of the endpoints are outside the unit square $\mathcal{S}$. Thus, the indicator for crossing one or both of the boundary lines can be written as:
$$H \equiv \mathbb{I} \Big( (\mathbf{E}_1 \notin \mathcal{S}) \vee (\mathbf{E}_2 \notin \mathcal{S}) \Big).$$
This gives an algorithm to generate a single indicator value for the needle crossing one or more boundaries of the grid-square. 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}((4 \ell - \ell^2) / \pi)$, so we can estimate $\pi$ by generating a large number of indicator values $H_1,...,H_n$ and taking:
$$\hat{\pi}_n = \frac{4 \ell - \ell^2}{\bar{H}_n},$$
where $\bar{H}_n$ is the sample mean of the indicators. For the optimal case where $\ell = 1$, Perlman and Wichura show that this estimator has asymptotic variance $\mathbb{V}(\hat{\pi}_n) \simeq 0.466 / n$.
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 data frame of values, with one row for each needle. For each needle we have columns giving the coordinates of the endpoints, and an indicator value for whether the needle crosses one or more gridlines.
BUFFON_LAPLACE_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 vectors for needle coordinates and cross indicator
E11 <- rep(0, n); E12 <- rep(0, n);
E21 <- rep(0, n); E22 <- rep(0, n);
CROSS <- rep(0, n);
#Generate needle
for (i in 1:n) {
#Generate midpoint and direction of needle
M <- runif(2);
D <- c(1,1);
while (norm(D, type = '2') > 1) { D <- runif(2); }
DD <- D/norm(D, type = '2');
#Determine endpoints and coordinates of needle
E1 <- M + (l/2)*DD; E2 <- M - (l/2)*DD;
E11[i] <- E1[1]; E12[i] <- E1[2];
E21[i] <- E2[1]; E22[i] <- E2[2];
#Determine whether needle crosses lines
CROSS[i] <- (min(E1) < 0)|(max(E1) > 1)|(min(E2) < 0)|(max(E2) > 1); }
#Generate data frame
data.frame(E11 = E11, E12 = E12, E21 = E21, E22 = E22, Cross = CROSS); }
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 <- 1;
n <- 10^7;
#Generate indicators
set.seed(1);
HHH <- BUFFON_LAPLACE_NEEDLES(l, n);
#Estimate pi
PI_EST <- (4*l - l^2)/mean(HHH$Cross);
PI_EST;
[1] 3.14153
PI_EST - pi;
[1] -6.231363e-05
Taking $n = 10^7$ gives an asymptotic standard error $\hat{\text{se}}_n \simeq 0.000215$, which is sufficient to get quite a good estimate of $\pi$. In this particular case, the error in the estimate is $6.231363 \times 10^5$.
$^\dagger$ It is worth noting here that one can obtain a Monte Carlo estimator for $\pi$ using the fact that the rejection-sampling algorithm has acceptance probability $\pi/4$. This is obviously much simpler than estimation by the needle algorithm. This merely reflects the fact that needle algorithms are not an efficient way to simulate an estimate of $\pi$.
