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.

Answer 1

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.

Answer 2

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.