(In honor of Pi Day)
Buffon’s needle is a mathematical problem posed by George-Louis Leclerc, Comte de Buffon, in the 18th century. It’s considered to be the first problem of geometric probability to be solved, and it can be used as a basis for a Monte Carlo method for calculating π. The resulting method is not one of the best ways to calculate π, but it is simple and illustrates the use of Monte Carlo methods.
The problem can be stated as such: (from the Wikipedia entry)
Suppose we have a floor made of parallel strips of wood, each the same width, and we drop a needle onto the floor. What is the probability that the needle will lie across a line between two strips?
We will think about the experiment in a “probability as frequency of occurrence” way and consider n repetitions of throwing a needle on the floor. As the figure shows, we let L be the length of the needles and d the width of the floor strips. Consider a single needle thrown in the ground; we will call X the distance between the center of the needle and the nearest strip line, and the smallest angle between the needle and the floor lines.
Because is the smallest angle between the needle and the floor line, it will always be acute, so , and because X is the distance to the nearest floor line, it’s at most half the distance between lines (which is the strip width d), so . The needle will cross a line if the center of the needle is closer to the nearest line than to the end of the needle in the direction perpendicular to the strip lines. X is the distance between the center of the needle and the nearest floor line, and the distance (in the direction perpendicular to the lines) between X and the end of the needle is given by (see figure). So the needle will cross a line if
Following the idea of the probability of an event as a measure of the event set in relation to the sample space, the probability of a needle crossing a line can be seen as the area under the curve for from 0 to π/2, divided by the whole space of values for y and . We want X < y and, as said before, X can be at most d/2, and is at most π/2, so the whole space of values is a rectangle with sides d/2 and π/2 (see figure below).
So the probability of a single cross is
which can be calculated to be
This last equation is valid only if d is not smaller than L, but this restriction makes sense for the experiment. Now comes the interesting part: make d = 2L and the probability of a cross becomes
The idea is to use a Monte Carlo method to estimate P(C), and this value is the reciprocal or our estimate for π. So, in Monte Carlo fashion, we will simulate the experiment many times, using any values for d and L such that d = 2L, and estimate the expected value of crosses in n repetitions of the experiment. The simulation draws samples from two uniform distributions: X from 0 to d/2 and from 0 to π/2, and then count the number of repetitions such that . Then, dividing the number of crosses in the simulation by the total number of repetitions, we get an estimate of the probability of a cross, which is the reciprocal of our estimate for π.
We have one last issue here: to sample , the usual way would be to get a sample from a Uniform[0,1] random number generator and then multiply by π/2, which requires us to know the value of π already. If we performed the actual experiment, with a real needle and a floor, we would not need the value of π. We will use a small cheat here and assume we have a function for computing the inverse of the sine function and create a variable halfpi that is equal to . Here’s the code in Julia for estimating π using this idea:
The results are not very good though. Even after millions of iterations, depending on the state of the RNG, I don’t get much more than one or two decimal places correct. In contrast to MATLAB, in Julia vectorized code doesn’t run faster, so here is an iterative version:
Finally, a version in OCaml. Of course, it would be easy to translate the iterative version to OCaml using imperative code, but just to do it differently here’s a pure functional version:
Results aren’t that great too. After almost 1 billion iterations I got π = 3.1415627. There are better ways to calculate π. There are even more efficient ways to use the idea of the Buffon’s needle experiment (see the end of the Wikipedia entry). A statistical analysis of the generated samples would also be beneficial, but this post is already running long. I any case, the idea is neat and serves as an illustration of an interesting application of Monte Carlo methods.