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


Imaginary Powers

Imaginary numbers ($i = \sqrt{-1}$) are used to relate Cartesian coordinates to polar coordinates as shown in the figure below. The argument $\varphi$ is the angle from the real axis along the unit circle, and can be expressed as two terms: a real and an imaginary, giving the Cartesian (x,y) position.
$e^{i\varphi} = \text{cos}(\varphi) + i \text{sin}(\varphi) = x + i y$
(If you want to expand beyond the unit circle just multiply by the desired radius: $e^{i\varphi} \rightarrow re^{i\varphi}$)

This is known as the Euler relation, and it gives us some neat results, like:
$e^{i\pi} = -1$
$e^{i \pi /2} = i$
But what happens when you exponentiate $i$ with itself? $i^i$.

From the Euler relation we know that $i=e^{i \pi /2}$, so it’s really a matter of substitution:
$i^i = e^{i \pi / 2 i} = e^{-\pi / 2} = 0.20788...$
So $i^i$ is a real number!

But wait a second, this isn’t the whole answer. We can add or subtract $2\pi$ from any angle and wind up in the same place on the circle. So $i = e^{i \pi /2 + i 2 \pi n}$, where $n$ is any integer. So we end up with a set of values:
$i^i = e^{- \pi /2 - 2 \pi n}$
which are all real!

The curve above goes through the set of values for $i^i$, where $n$ is an integer.
$i^i = \{...,111.32,0.20788,0.00038820,...\}$
But 0.20788… will be looked at as the primary value for now since it is the result when $n=0$, and having one value keeps things simple.

So what happens when we repeat this exponentiation? $i^{i^i}$
First, this can be written in a simpler way:
$^3i = i^{i^i}$
This operation is called “tetration” or “power tower”.
You start from the top and work your way down. So $^3i = i^{i^i} = i^{0.20788}$
Now use the Euler relation again (ignoring $2 \pi n$ terms in the exponential):
$i^{0.20788} = e^{i \pi/2 \times 0.20788} = e^{i 0.32653}$
So the result is complex:
$^3i = 0.94715+0.32076i$

We can keep going using results from the last iteration:
$^4i = i^{0.94715+0.32076 i} = 0.05009+0.60211 i$

Eventually you find that with $^mi$ as $m$ increases, the results start to converge on one number. This implies there’s a value for $^\infty i$. After plugging and chugging we can eventually conclude that:
$^\infty i = 0.43828 + 0.36059 i$
Convergence can be seen in both the real and imaginary components, around $m \approx 70$ if we require the difference between iterations to be less than 0.1%. (See figures below).

Each result has a different radius associated with it, since the imaginary component of the previous result puts a real term in the exponential of the next. This real term can be expressed as a factor $r$:
$e^{i\varphi} = a + bi \rightarrow i^{a+bi} = e^{i \pi/2 (a + bi)} = e^{-b\pi/2}e^{i a \pi/2} = r e^{i a \pi/2}$
The radius can be found by multiplying the result by its complex conjugate and then find the square root (i.e. the Pythagorean theorem).

We can then plot radius of $^mi$ as a function of $m$ (see plot below).

The python code for making these plots is:

import os
import sys
import math
import numpy
import pylab

REAL = []
IMAG = []
ITERA = []

pp = 1j
for n in range(100):
ppL = pp #save previous value for cut off
pp = (1j)**pp #interate power tower
RA = 0.001 #ratio for cut off
if (abs(pp.real - ppL.real)/(ppL.real + 0.00001) < RA) and (abs(pp.imag - ppL.imag)/(ppL.imag + 0.00001) < RA): #establish a cut off
print n
print pp
break
else:

REAL.append(pp.real) #save real component
IMAG.append(pp.imag) #save imaginary component
ITERA.append(n+2) #save iteration "m"

ITERA = [1]+ITERA
REAL = [0]+REAL
IMAG = [1]+IMAG

pylab.figure(1)
pylab.xlabel("m")

# pylab.show()

pylab.figure(2)
pylab.xlabel("m")
pylab.ylabel("imaginary component")
pylab.plot(ITERA,IMAG,'k')

pylab.savefig('Imfig.png')

pylab.figure(3)
pylab.xlabel("m")
pylab.ylabel("real component")
pylab.plot(ITERA,REAL,'k')

pylab.savefig('Refig.png')