Wednesday, 30 December 2015

Buffon Needle Problem

Buffon Needle problem is a classical mathematical problem first suggested on 18th century by the French naturalist Georges-Louis Leclerc, Comte de Buffon. Although it was not the original aim of Buffon, the problem turned out to be a method to estimate the number π.

The problem is as follows: we have a set of evenly distanced parallel bars on a board. If we drop a needle on the board, what is the probability that the needle will touch one of the bars?

Assume the distance between the bars is t, and the length of the needle is l. Assume also that t is greater than l. If that is not the case, the problem can still be easily solved but the solution is different.

This problem is supposed to the very first problem of geometrical probability, and can be easily solved using simple integral calculus. At the same time, it turns out that the ubiquitous number π is in the solution of the problem, so, it can be equally used as a Monte-Carlo estimation of this number. Although it is not the most efficient way to estimate the digits of π, we will solve the problem and write a Python program to check it, just for the fun of it.

Assume $x$ is the distance from the middle point of the needle to the closest bar, then assume, between any two bars, $x \in [0, \frac{t}{2}]$. Assume also that $\theta$ is the angle of the needle with the bar. Then, $x, \theta$ will be uniform random variables for each of the random needle throw. The p.d.f of each variable will be $f_X(x) = \frac{2}{t} = , f_{\Theta}(\theta) = \frac{2}{\pi}$. Since both variables will be independent, the joint distribution density probability function will be just the product of the two, i.e. $f_{X,\Theta}(x, \theta) = \frac{4}{\pi t}$.

Then, once we have drawn the $x, \theta$ random variables at each trial, the needle will touch the bar if the condition $x < \frac{l}{2}sin ( \theta )$.

Thus, the probability that the needle touches a bar will be $P = \int_{\theta = 0}^{\frac{\pi}{2}} \int_{x=0}^{\frac{l}{2}sin(\theta)}  f_{X, \Theta}(x, \theta) dx d \theta$. With a little of calculus, we get $  P = \frac{4}{\pi t} \int_{\theta = 0}^{\frac{\pi}{2}} \frac{l}{2}sin(\theta) d \theta = \frac{2 l}{\pi t}$.

So, the probability of the needle touching one bar is surprisingly related to number $\pi$ through the expression $P = \frac{2 l}{\pi t}$.

Assume that we throw indeed $n$ needles in a simulation, then, let $p$ be the number of trials in which the needle will touch one of the bar. An estimate of the probability above, will be $\frac{p}{n}$ , i.e. $\frac{p}{n} \approx \frac{2 l}{\pi t}$.

Now we have a relation between number $\pi$ and a probabilistic simulation, we can write a snippet of Python code to run simulations and see how good the estimates of number $\pi$ are using this method.

The Python code is included at the end of this blog post. The code will simulate $m$ times the Buffon problem, each problem throwing $n$ needles in a three bars board. The image resulting from the Python code gives us an idea of the needles touching the bars, drawn in red, and the ones not touching any bar, drawn in green.

Example of simulation with n = 100, command line argument:
python -n 150 -m 1

Executing the program with the parameter $m = 1000$, we get the mean and standard deviation of this estimator for number $\pi$.

n Mean Std
100 3.21196453418 0.502342987432
1000 3.16177668689 0.152438640479
10000 3.15099821931 0.0462329285661

As for the result, we see that it might not be the most efficient estimator for number $\pi$, anyway, we had a lot of fun thinking and solving the problem, which was the point of all this.

The code, as usual, can be found on my GitHub page.


No comments:

Post a Comment