Tag Archives: pi

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.52.6, shown in the figure below (\pi shown in red). Why is this the case?
DPI1

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!
BuffonS2
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:
DPI2

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.
DPI3

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()
Advertisements

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.

Buffon3
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.

buffon5
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.

BPI

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.

BPI2

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:

BPI1q

BPI2q

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)>0.5): #balance probability density to make it more uniform in terms of acute angle
        SIGN = 1
        if (random.uniform(0.0,1.0)>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 > 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("NT")
pylab.ylabel("2 NT/NC")
pylab.plot(NNa,PIc,'r')
pylab.plot(NNa,PIa,'k')
pylab.savefig('BPI1q.png')

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

pylab.show()