Python Programming for Biology: Bioinformatics and Beyond

Figure 25.2.  Using the Monte Carlo method to estimate π from the area of a circle.

By uniformly selecting random points in a square interval we can estimate the area of a

circle from the proportion of points that fall within a given distance of the centre. This in

turn can be used to give an estimate for the mathematical constant π: for a circle of radius

1, and hence area π, the area of the bounding square is 4. Thus the ratio of circle area

points to total points approaches π/4 as the sample size increases. This is a simple

example of the Monte Carlo approach, which uses random samples of test points to solve

a problem.

Firstly,  we  import  the  random  module  from  NumPy  and  make  a  few  definitions.  The

uniform  function  is  defined  upfront  so  that  we  are  not  repeatedly  using  the  dot  notation

inside  loops  (which  is  slower).  This  will  generate  pseudorandom  numbers  with  an  even

distribution in our desired range. The number of random points sampled for the example is

defined as 100,000, but it is interesting to play with this to see the effect on the accuracy

of estimating π. The number of points inside the circle numInside is naturally defined as

zero at the start.

from numpy import random

uniform = random.uniform

numSamples = 100000

numInside = 0

Then  comes  a  loop  with  the  required  number  of  cycles  where  we  define  the  random

points.  Here  we  use  the  range()  function,  which  works  in  both  Python  2  and  Python  3,

although for the latter, rather than generating a whole list of values, it creates an iterable

object that gives the sequence of numbers on demand (you can use xrange() in Python 2

for this behaviour). Note that the first two arguments to the function uniform() define the

limit  to  the  random  values  (i.e.  between  plus  and  minus  one)  and  the  last  argument

indicates how many values to yield; here we need two, for x and y.

---

for i in range(numSamples): # could use xrange in Python 2

x, y = uniform(-1.0, 1.0, 2)

With  the  coordinates  of  the  random  point  defined  we  now  test  whether  the  point  is

inside the circle; if the sum of coordinates squared is less than one (squared) the point is

sufficiently  close  to  the  centre.  If  the  coordinates  represent  a  point  for  which  this  holds

then the count for internal points is increased by one.

if (x * x) + (y * y) < 1.0:

numInside += 1

At  the  end  the  estimate  for  π  is  simply  four  times  the  proportion  of  internal  points,

relative to the total number of points. This follows from the area of the region being 4.0,

i.e. a circle of radius 1.0 is bounded by a square with sides of length 2.0.

pi = 4.0 * numInside / float(numSamples)

print(pi)

This  is  not  an  especially  efficient  means  of  estimating  π

4

,  but  the  same  idea  can  be

applied to any shape or function that can be defined, however many dimensions it might

have.