# More Buffon’s Needle

In my previous blog entry, I discussed Buffon’s needle (a method of measuring $\pi$ using random numbers). I tried to recreate what would occur in physical reality without using trig functions in the measurement of $\pi$. The tricky part was trying to make a uniform probability distribution along the angle $\theta$, which is what you expect with physical needles being dropped. In this entry, I will do something more exact, but less physical.

In reality, the probability density in $\theta$ is uniform due to there being no preferential direction for a needle to point. To simulate this with uniform random numbers without using trigonometric functions (to avoid just measuring python’s $\pi$) is tricky. The way I did this was assign uniform random numbers (uniform probability density) to a projection of a needle, either perpendicular or parallel to the lines (with a probability of 0.5 assigned to each case). This results in a probability distribution in $\theta$ that isn’t quite uniform, but it’s close enough to approximate $\pi$ after a large number of needles are dropped:
$2 \frac{N_T}{N_C} = \pi$
where $N_T$ and $N_C$ are the total number of needles dropped and number of needles crossing a line, respectively.

If instead of assigning a uniform probability to either projection fairly, what if we assign it to just one? This was shown in the previous post. As the number of needles dropped increases, the ratio $2 \frac{N_T}{N_C}$ approaches some number around $2.5$$2.6$, shown in the figure below ($\pi$ shown in red). Why is this the case?

The reason why has to do with a uniform probability density along a projection ($q$) does not translate into a uniform probability density along $\theta$, see figure below. This means needles have a preferred direction. It turns out this can still be used to extract a measurement of $\pi$!

The uniform probability density is assigned to the projection $q$:
$\rho_q = 1/L$, for $0 \leq q \leq L$, and $\rho_q=0$ elsewhere.
To get this in terms of $\theta$, some substitutions are made in the probability density integral:
$\int_0^L \rho_q dq = 1$
The variables $q$ and $\theta$ are related by:
$\text{cos}(\theta) = q/L$
so by taking the derivative we get the substitution:
$dq = -L\text{sin}(\theta) d \theta$

So now the probability density integral can be written as:
$\int_0^{\pi/2} \rho_q L \text{sin}(\theta) d \theta$
which is also:
$= \int_0^{\pi/2} \rho_{\theta} d \theta$
So it gives the form of the probability density for $\theta$ in terms of the one for $q$:
$\rho_{\theta} = \rho_q L\text{sin}(\theta)$

Plug in $\rho_q = 1/L$ to get:
$\rho_{\theta} = \text{sin}(\theta)$
and it is obviously NOT uniform, it’s just the sine function!
Plotting a histogram of the number of needles per angle $\theta$ will give the shape of the sine function:

Following the same steps as in the previous entry, we can find out why when a uniform probability is assigned to one projection, $2 N_T/N_C$ approaches $2.5$ instead of $\pi$.
For simplicity $W=L$, and again for a needle to be crossing a line it has to be below the upper bound:
$x \leq \frac{L}{2}\text{sin}(\theta)$
and the probability of a needle crossing a line is:
$P = \int_0^{\frac{L}{2}\text{sin}(\theta)} \int_0^{\pi/2} \rho_x \rho_{\theta} d \theta d x$
Plugging in the probability densities: $\rho_x = 2/L$ and $\rho_\theta = \text{sin}(\theta)$:
$P = \int_0^{\pi/2} \text{sin}^2(\theta) d \theta = \pi/4$

$P = \frac{N_C}{N_T}$ when the number of needles dropped is large, so by combining these equations:
$2 \frac{N_T}{N_C} = \frac{8}{\pi} = 2.52$
This is the number that is being approached!
But wait, $\pi$ is there, so we can extract it by:
$4 \frac{N_C}{N_T} = \pi$

By plotting this ratio as a function of number of needles dropped, we see it approaching $\pi$. See the figure below.

Even though this isn’t how needles behave in reality, it’s a simpler way to measure $\pi$ by using random numbers in a simulation. The code used to make these plots is shown below.

By the last drop:
$2 \frac{N_T}{N_C} = 2.54$
$4 \frac{N_C}{N_T} = 3.14$

import os
import sys
import math
import numpy
import pylab
import random

CL = 0
NN = 0

PIa = []
PIb = []
NNa = []
PIc = []
ANGLES = []

for ii in range(100000):

PPy = random.uniform(0.0,1.0) #A-end y position (lines separated by unit distance)

RRx = random.uniform(0.0,1.0) #x projection
RRy = numpy.sqrt(1.0 - RRx**2) #y projection

QQy = PPy + RRy #B-end y position
QQyf = math.floor(QQy)
PPyf = math.floor(PPy)

NN = NN + 1
if (QQyf-PPyf > 0.0): #if crossing line
CL = CL + 1

NNa.append(NN)
PIc.append(math.pi)
RAT = (1.0*CL)/(1.0*NN)
if (CL==0):
PIa.append(0)
else:
PIa.append(2.0/RAT)

PIb.append(4*RAT)

ANGLES.append(abs(math.atan(RRy/(RRx+0.0001)))/math.pi)

print PIa[-1]
print PIb[-1]

pylab.figure(1)
pylab.xlabel("NT")
pylab.ylabel("2 NT/NC")
pylab.plot(NNa,PIc,'r')
pylab.plot(NNa,PIa,'k')
pylab.savefig('DPI1.png')

pylab.figure(2)
pylab.hist(ANGLES, bins=10)
pylab.xlabel("Acute angle in factors of pi")
pylab.ylabel("Needles dropped")
pylab.savefig('DPI2.png')

pylab.figure(3)
pylab.xlabel("NT")
pylab.ylabel("4 NC/NT")
pylab.plot(NNa,PIc,'r')
pylab.plot(NNa,PIb,'k')
pylab.savefig('DPI3.png')

pylab.show()


# Happy Pi Day!

Since it is $\pi$ day, I thought I would show a way to measure $\pi$. The way I have in mind is known as Buffon’s needle. If you have a flat surface with parallel lines of uniform separation and drop needs of uniform length on it, the ratio between the total number needles dropped and the number of needles crossing a line will be a factor of $\pi$.

The figure below shows the needles of length $L$ in blue, and the lines in black with separation $W$. For each needle there are two important values: $x$ and $\theta$. The $x$ value is the distance between the center of the needle and the nearest line. The $\theta$ value is the acute angle between the needle and the nearest line. The span of available values of $x$ is then $0 \leq x \leq W/2$, and the span of $\theta$ is $0 \leq \theta \leq \pi/2$.

Considering all possible values of $x$ and $\theta$, the probability is 1 by definition.
$P = \int_{\Delta \theta} \int_{\Delta x} \rho dx d\theta = 1$
Where in general, if $x$ and $\theta$ are independent: $\rho = \rho(x)\rho(\theta)$.

A uniform probability densities are assumed for both $x$ and $\theta$ since there is no preferred direction or position for needles. So the probability densities $\rho$ are defined by:
$\rho_{\theta} = 2/\pi$, for $0 \leq \theta \leq \pi/2$, and $\rho_{\theta} = 0$ elsewhere.
$\rho_{x} = 2/W$, for $0 \leq x \leq W/2$, and $\rho_{x} = 0$ elsewhere.
A uniform probability density is shown in the figure below for $\theta$, note that the area is 1.

For simplicity, we will only consider the case where the needles are shorter or the same length as the separation distance between the lines: $L \leq W$.
So for needles crossing a line, the projected half-length of the needle perpendicular to the lines is greater than the distance between the center of the needle to the nearest line:
$x \leq \frac{L}{2}\text{sin}(\theta)$
This is the upper bound that is then plugged into the probability integral, to get the probability of crossing a line:
$P = \int_{0}^{\pi/2} \int_{0}^{\frac{L}{2}\text{sin}(\theta)} \rho_{x} \rho_{\theta} dx d\theta = \frac{2 L}{\pi W}$
For the simple case where $L = W$:
$P = \frac{2}{\pi}$
After enough drops the ratio between the number of needles crossing a line $N_C$ and the total number of needles $N_T$ will resemble the probability:
$P = \frac{N_C}{N_T}$
Combining these two equations, we get ratio that can measure $\pi$:
$2 \frac{N_T}{N_C} = \pi$

Doing this measurement in real life by dropping needles would make a mess, so I wrote a python code to it instead. Below is a graph showing this ratio as a function of needle drops (shown in black), $\pi$ is shown in red. You can see the black line approaches $\pi$ as the number of needles increases.

The probability density as a function of $\theta$ is shown below. It is not perfectly flat as in the ideal case. This is because in my code the uniform probability density is along the Cartesian axes ($X,Y$) rather than along $\theta$. (The parallel lines are separated along $Y$). I did this to avoid measuring python’s $\pi$ in its trigonometric functions.

A uniform random number is assigned to the projected needle length along $X$ or $Y$ with a range of 0 to 1. The length of a needle is forced to be 1 by using the Pythagorean theorem: $X = \sqrt{1-Y^2}$ or $Y = \sqrt{1-X^2}$. There is a 0.5 probability that either $X$ or $Y$ is assigned to the uniform random number. This approximates a uniform probability density in $\theta$. It’s not perfect though.

If however, the uniform random number is assigned to only $X$, the ratio doesn’t make it to $\pi$, it stops around 2.6, and the probability density is far from uniform. See the figures below:

The code for these plots is shown below:

import os
import sys
import math
import numpy
import pylab
import random

CL = 0
NN = 0

PIa = []
NNa = []
PIc = []
ANGLES = []

for ii in range(10000):

PPx = random.uniform(1.0,10.0) #1-10 A-end x position
PPy = random.uniform(0.0,10.0) #1-10 A-end y position (lines)

RRx = random.uniform(-1.0,1.0) #B-end x displacement
RRy = numpy.sqrt(1.0 - RRx**2) #B-end y displacement

if (random.uniform(0.0,1.0)&amp;amp;amp;amp;amp;amp;gt;0.5): #balance probability density to make it more uniform in terms of acute angle
SIGN = 1
if (random.uniform(0.0,1.0)&amp;amp;amp;amp;amp;amp;gt;0.5):
SIGN = -1
RRy = random.uniform(0.0,1.0) #B-end y displacement
RRx = numpy.sqrt(1.0 - RRy**2)*SIGN #B-end x displacement

QQx = PPx + RRx #B-end x position
QQy = PPy + RRy #B-end y position
QQyf = math.floor(QQy)
PPyf = math.floor(PPy)

NN = NN + 1
if (QQyf-PPyf &amp;amp;amp;amp;amp;amp;gt; 0.0): #if crossing line
CL = CL + 1

NNa.append(NN)
PIc.append(math.pi)
RAT = (1.0*CL)/(1.0*NN)
if (CL==0):
PIa.append(0)
else:
PIa.append(2.0/RAT)

ANGLES.append(abs(math.atan(RRy/(RRx+0.0001)))/math.pi)

print PIa[-1]

pylab.figure(1)
pylab.xlabel(&amp;amp;amp;amp;amp;amp;quot;NT&amp;amp;amp;amp;amp;amp;quot;)
pylab.ylabel(&amp;amp;amp;amp;amp;amp;quot;2 NT/NC&amp;amp;amp;amp;amp;amp;quot;)
pylab.plot(NNa,PIc,'r')
pylab.plot(NNa,PIa,'k')
pylab.savefig('BPI1q.png')

pylab.figure(2)
pylab.hist(ANGLES, bins=10)
pylab.xlabel(&amp;amp;amp;amp;amp;amp;quot;Acute angle in factors of pi&amp;amp;amp;amp;amp;amp;quot;)
pylab.ylabel(&amp;amp;amp;amp;amp;amp;quot;Needles dropped&amp;amp;amp;amp;amp;amp;quot;)
pylab.savefig('BPI2q.png')

pylab.show()