In my previous blog entry, I discussed Buffon’s needle (a method of measuring using random numbers). I tried to recreate what would occur in physical reality without using trig functions in the measurement of . The tricky part was trying to make a uniform probability distribution along the angle , 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 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 ) 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 that isn’t quite uniform, but it’s close enough to approximate after a large number of needles are dropped:

where and 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 approaches some number around –, shown in the figure below ( shown in red). Why is this the case?

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

The uniform probability density is assigned to the projection :

, for , and elsewhere.

To get this in terms of , some substitutions are made in the probability density integral:

The variables and are related by:

so by taking the derivative we get the substitution:

So now the probability density integral can be written as:

which is also:

So it gives the form of the probability density for in terms of the one for :

Plug in to get:

and it is obviously NOT uniform, it’s just the sine function!

Plotting a histogram of the number of needles per angle 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, approaches instead of .

For simplicity , and again for a needle to be crossing a line it has to be below the upper bound:

and the probability of a needle crossing a line is:

Plugging in the probability densities: and :

when the number of needles dropped is large, so by combining these equations:

This is the number that is being approached!

But wait, is there, so we can extract it by:

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

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

By the last drop:

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