Monte Carlo estimate for pi with numpy

In this post we will use a Monte Carlo method to approximate pi. The idea behind the method that we are going to see is the following:

Draw the unit square and the unit circle. Consider only the part of the circle inside the square and pick uniformly a large number of points at random over the square. Now, the unit circle has pi/4 the area of the square. So, it should be apparent that of the total number of points that hit within the square, the number of points that hit the circle quadrant is proportional to the area of that part. This gives a way to approximate pi/4 as the ratio between the number of points inside circle and the total number of points and multiplying it by 4 we have pi.

Let's see the python script that implements the method discussed above using the numpy's indexing facilities:

from pylab import plot,show,axis
from numpy import random,sqrt,pi

# scattering n points over the unit square
n = 1000000
p = random.rand(n,2)

# counting the points inside the unit circle
idx = sqrt(p[:,0]**2+p[:,1]**2) < 1

plot(p[idx,0],p[idx,1],'b.') # point inside
plot(p[idx==False,0],p[idx==False,1],'r.') # point outside
axis([-0.1,1.1,-0.1,1.1]) 
show()

# estimation of pi
print '%0.16f' % (sum(idx).astype('double')/n*4),'result'
print '%0.16f' % pi,'real pi'

The program will print the pi approximation on the standard out:

3.1457199999999998 result
3.1415926535897931 real pi

and will show a graph with the generated points:

Note that the lines of code used to estimate pi are just 3!

Approximating pi

This script uses the following formula

to approximate the value of pi with a fixed number of correct digits.

import math
def pi(digits):
 k = 0
 pi = 0.0
 e = 1.0
 tol = pow(10,-digits) # is minimum error that we want to accept
 while e > tol:
  pi_old = pi
  pi += (2*pow(-1,k)*pow(3,.5-k))/(2*k+1)
  e = abs(pi-pi_old) # current error
  k += 1
  print '%0.16f %20.16f' % (pi,e)
 print '\nerror',e,'\niterations ',k
 print '%0.16f' % pi,'result'
 print '%0.16f' % math.pi,'real pi'

pi(8) # approximating pi with 8 correcet digits

During the execution you can see the current approximated value on the left column and the difference with the approximation at the previous step on the right column. At the end, the difference between the approximated value and the value provided by the math python module will be printed.

$ python pi.py 
3.4641016151377544   3.4641016151377544
3.0792014356780038   0.3849001794597506
3.1561814715699539   0.0769800358919501
3.1378528915956800   0.0183285799742738
3.1426047456630846   0.0047518540674045
3.1413087854628832   0.0012959602002014
3.1416743126988376   0.0003655272359544
3.1415687159417840   0.0001055967570536
3.1415997738115058   0.0000310578697218
3.1415905109380802   0.0000092628734256
3.1415933045030817   0.0000027935650015
3.1415924542876463   0.0000008502154354
3.1415927150203800   0.0000002607327336
3.1415926345473140   0.0000000804730660
3.1415926595217138   0.0000000249743999
3.1415926517339976   0.0000000077877162

error 7.78771624965e-09 
iterations  16
3.1415926517339976 result
3.1415926535897931 real pi