Monte Carlo Method in Python

Question

I've been attempting to use Python to create a script that lets me generate large numbers of points for use in the Monte Carlo method to calculate an estimate to Pi. The script I have so far is this:

import math
import random
random.seed()

n = 10000

for i in range(n):
    x = random.random()
    y = random.random()
    z = (x,y)

    if x**2+y**2 <= 1:
        print z
    else:
        del z

So far, I am able to generate all of the points I need, but what I would like to get is the number of points that are produced when running the script for use in a later calculation. I'm not looking for incredibly precise results, just a good enough estimate. Any suggestions would be greatly appreciated.

Do you wish to count how many random pairs are inside the circle? If that's the case just use a counter... — zenpoy, Nov 19 '12 at 20:28

score 15 · Answer 1 · answered Nov 19 '12 at 20:44

If you're doing any kind of heavy duty numerical calculation, considering learning numpy. Your problem is essentially a one-linear with a numpy setup:

import numpy as np

N   = 10000
pts = np.random.random((N,2))

# Select the points according to your condition
idx = (pts**2).sum(axis=1)  < 1.0
print pts[idx], idx.sum()

Giving:

[[ 0.61255615  0.44319463]
 [ 0.48214768  0.69960483]
 [ 0.04735956  0.18509277]
 ..., 
 [ 0.37543094  0.2858077 ]
 [ 0.43304577  0.45903071]
 [ 0.30838206  0.45977162]], 7854

The last number is count of the number of events that counted, i.e. the count of the points whose radius is less than one.

score 2 · Answer 2 · answered Nov 19 '12 at 20:26

Not sure if this is what you're looking for, but you can run enumerate on range and get the position in your iteration:

In [1]: for index, i in enumerate(xrange(10, 15)):
   ...:     print index + 1, i
   ...:
   ...:
1 10
2 11
3 12
4 13
5 14

In this case, index + 1 would represent the current point being created (index itself would be the total number of points created at the beginning of a given iteration). Also, if you are using Python 2.x, xrange is generally better for these sorts of iterations as it does not load the entire list into memory but rather accesses it on an as-needed basis.

@ahans No problem - I remember seeing it for the first time and thinking 'Ah, figured there was a way to do that.' Have fun using it! — RocketDonkey, Nov 19 '12 at 20:31

score 1 · Answer 3 · answered Nov 19 '12 at 20:34

Just add hits variable before the loop, initialize it to 0 and inside your if statement increment hits by one.
Finally you can calculate PI value using hits and n.

import math
import random
random.seed()

n = 10000
hits = 0  # initialize hits with 0

for i in range(n):
    x = random.random()
    y = random.random()
    z = (x,y)

    if x**2+y**2 <= 1:
        hits += 1
    else:
        del z

# use hits and n to compute PI

Monte Carlo Method in Python

3 Answers3