Sunday, February 17, 2013

Estimating Pi Using Random Samples

Programming languages usually come with defined values of pi. However, if you want to estimate it yourself, there is a very simple way to do so.
The area of a circle is given by \[A_c = \pi r^2\]
The area of a square is given by \[A_s = a^2\]
Choose the circle to be the unit circle and the square to be enclosing the circle and draw random numbers that are uniformly distributed over the square. By calculating the ratio of the number n of samples within the circle and all samples N you can get an estimation for pi. a and r are chosen according to above descriptions.
\[pi = \frac{n}{N}\frac{a^2}{r^2} = \frac{n}{N}\frac{2^2}{1^2} = 4\frac{n}{N}\]
In Python numpy this can be done within a few lines of code, in fact you can even condense it to one line.
import numpy as np

def estimate_pi(N):
    rand = 2.0*np.random.rand(2, N) - 1.0
    dist = rand[0,:]*rand[0,:] + rand[1,:]*rand[1,:]
    return np.sum(dist <= 1)*4.0/N

You should choose a large N, as you will get a better and more stable estimate:

Calculating pi using 1000 random samples
pi = 3.092
derivation from numpy.pi =  0.0157858319197

Calculating pi using 100000 random samples
pi = 3.1438
derivation from numpy.pi =  0.000702620184601

Calculating pi using 10000000 random samples
pi = 3.1424128
derivation from numpy.pi =  0.000261060710487

Note: Even if you do not know how to calculate the area of a circle, you still can calculate the area of a circle with an arbitrary radius using this method and then you can figure out that the area is a quadratic function of the radius.

No comments:

Post a Comment