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).
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)*4If 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