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!
Source: http://glowingpython.blogspot.com/2012/01/monte-carlo-estimate-for-pi-with-numpy.html
(Note: Opinions expressed in this article and its replies are the opinions of their respective authors and not those of DZone, Inc.)
Comments
Gene Leynes replied on Wed, 2012/01/25 - 4:45pm
This is a really nice example of how to use Python with numpy. Thanks for posting.
For fun, I made a comparison of the code in R. It only took a minute, and I thought it would be interesting to see the similarities and differences between the two languages. It's interesting to me because they are both very simple to code when compared to something like Java or a C language. However, I think that Python is generally going to run faster in many instances. So, it'a another vote for me that it's worthwhile to use Python and its libraries like numpy and scipy.
## MC Simulation of PI using R n = 1000000 p = matrix(runif(n*2), ncol=2) # random numbers / point locations ind = sqrt(p[,1]^2 + p[,2]^2) < 1 # is the point inside the circle? col = ifelse(ind, 'blue', 'red') # blue if inside, otherwise red plot(p, col=col, pch=16, cex=.1) # not as pretty as numpy! sum(rep(1,n)[ind]/n*4) # Simulation result # rep makes a vector of 1s with length n # only the 1s where ind == true are totaled pi # Real pi