Random
— 5 min read

Estimating Pi with Monte Carlo method.

Using the Monte Carlo method to estimate the value of pi.

Monte Carlo is statistical method based on a series of random numbers. Monte Carlos are used in a wide-range of physical systems, finance and other research field.

In this post, we are going to use the Monte Carlo method to estimate the value of pi.

First, consider a circle inscribed in a square (as in the figure).

circle-square

f assume that the radius of the circle is , then the Area of the circle and the Area of the squar .

If we throw a dart blindly inside of the square, what is the probability (P) the dart will actually land inside the circle?

P = Area of the circle / Area of the square = Pi / 4

So the chances of hitting the circle are Pi / 4.

In other words,

Let’s try to do this in python:

1 from __future__ import division
2 import numpy
3 
4 NBR_POINTS = 1000000
5 RADIUS = numpy.sqrt(NBR_POINTS)
6 
7 print (len(filter(lambda x:numpy.sqrt(numpy.random.randint(0, RADIUS)**2+numpy.random.randint(0, RADIUS)**2)<RADIUS,xrange(NBR_POINTS)))/NBR_POINTS)*4

If you have problem with the filter, read this!

To go a little bit further with this we can imagine a map-reduce system to make the estimation faster.

 1 from __future__ import division  
 2 import collections
 3 import itertools
 4 import multiprocessing
 5 import numpy
 6 
 7 class MapReduce(object):
 8     """
 9     The map reduce object, should be initialized with:
10         map_fn
11         reduce_fn
12         nbr_workers
13     """
14     def __init__(self, map_fn, reduce_fn, num_workers=None):
15         """
16         initiaize the mapreduce object
17             map_fn : Function to map inputs to intermediate data, takes as
18             input one arg and returns a tuple (key, value)
19             reduce_fn : Function to reduce intermediate data to final result
20             takes as arg keys as produced from the map, and the values associated with it
21         """
22         self.map_fn = map_fn
23         self.reduce_fn = reduce_fn
24         self.pool = multiprocessing.Pool(num_workers)
25 
26     def partition(self, mapped_values):
27         """
28         returns the mapped_values organised by their keys. (keys, associated values)
29         """
30         organised_data = collections.defaultdict(list)
31         for key, value in mapped_values:
32             organised_data[key].append(value)
33         return organised_data.items()
34 
35     def __call__(self, inputs=None, chunk_size=1):
36         """
37         process the data through the map reduce functions.
38         inputs : iterable
39         chank_size : amount of data to hand to each worker
40         """
41         mapped_data = self.pool.map(self.map_fn, inputs, chunksize = chunk_size)
42         partioned_data = self.partition(itertools.chain(*mapped_data))
43         reduced_data = self.pool.map(self.reduce_fn, partioned_data)
44         return reduced_data
45 
46 
47 NBR_POINTS = 1000000
48 RADIUS = numpy.sqrt(NBR_POINTS)
49 NBR_WORKERS = 4
50 NBR_PER_WORKER = int(NBR_POINTS/NBR_WORKERS)
51 
52 
53 def probability_calculation(item):
54     print multiprocessing.current_process().name, 'calculating', item
55     output = []
56     output.append( ('pi', len(filter(lambda x:numpy.sqrt(numpy.random.randint(0, RADIUS)**2+numpy.random.randint(0, RADIUS)**2)<RADIUS,xrange(NBR_PER_WORKER)))) )
57     return output
58 
59 
60 def estimate_pi(items):
61     keys, occurrences = items
62     return (sum(occurrences)/NBR_POINTS)*4
63 
64 mapper = MapReduce(probability_calculation, estimate_pi)
65 pi = mapper([i for i in xrange(NBR_WORKERS)])
66 print pi