CODEWITHSUNDEEP
pythonmontecarlo
Matt Johnson
Matt Johnson

Reputation: 91

Monte Carlo Method in Python

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.

Upvotes: 9

Views: 20281

Answers (3)

Hooked
Hooked

Reputation: 88248

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.

Upvotes: 15

pm007
pm007

Reputation: 383

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

Upvotes: 1

RocketDonkey
RocketDonkey

Reputation: 37279

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.

Upvotes: 2

Related Questions