Velocity Reviews > Monte Carlo Method and pi

# Monte Carlo Method and pi

Karl Pech
Guest
Posts: n/a

 07-08-2004
Hi,

I got a task there I have to compute pi using the
Method above.

So I wrote the following program:
---
import random
import math

n = long(raw_input("Please enter the number of iterations: "))
sy = 0 # sum of the function-values

for i in range(0, n):
x = random.random() # coordinates
sy += math.sqrt(1-math.sqrt(x)) # calculate y and add to sy

print 4*sy/n # compute pi
---

Unfortunately, even for n = 2000000 the result is ~2,13... .
It converges very slow to pi and I don't know why.

Regards
Karl

Nick Smallbone
Guest
Posts: n/a

 07-08-2004

"Karl Pech" <(E-Mail Removed)> wrote in message
news:ccke8p\$c0h\$06\$(E-Mail Removed)-online.com...
> Hi,
>
>
> I got a task there I have to compute pi using the
> Method above.
>
> So I wrote the following program:
> ---
> import random
> import math
>
> n = long(raw_input("Please enter the number of iterations: "))
> sy = 0 # sum of the function-values
>
> for i in range(0, n):
> x = random.random() # coordinates
> sy += math.sqrt(1-math.sqrt(x)) # calculate y and add to sy

This line is wrong. I think it should be sy += math.sqrt(1 - x ** 2). Doing
that gives me ~3.1 for 1000 iterations.

Nick

Jeff Epler
Guest
Posts: n/a

 07-08-2004
Others spotted the problem with your implementation of the approximatoin
method.

You can greatly speed things up with numarray. It gets "3.141..." as
the approximation using n=2000000 in about 4 seconds on this sub-GHz
laptop. Of course, a method with faster convergence would get to 3.141
in a lot less than 4 seconds, even if written in pure Python.

Jeff

import numarray
import numarray.random_array

def approx_pi(n):
a = numarray.random_array.uniform(0, 1, n)
return 4 * numarray.sum(numarray.sqrt(1 - a**2)) / n

if __name__ == '__main__':
n = int(raw_input("Please enter the number of iterations: "))
print approx_pi(n)

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.2.4 (GNU/Linux)

iD8DBQFA7cbRJd01MZaTXX0RAoxMAJ4wCeN/z6C7LY7uXg21Y/TYYkpApACZAfBa
8w3R5vITWyFxLNBKG4dS10E=
=Ch/k
-----END PGP SIGNATURE-----

Paul Rubin
Guest
Posts: n/a

 07-08-2004
"Karl Pech" <(E-Mail Removed)> writes:
> I got a task there I have to compute pi using the
> Method above.

This sounds like a homework problem...
> for i in range(0, n):
> x = random.random() # coordinates
> sy += math.sqrt(1-math.sqrt(x)) # calculate y and add to sy
>
> print 4*sy/n # compute pi
> ---
>
> Unfortunately, even for n = 2000000 the result is ~2,13... .
> It converges very slow to pi and I don't know why.

Better make sure the formulas are correct, by asking yourself why the
simulation is supposed to work at all.

Karl Pech
Guest
Posts: n/a

 07-08-2004
Thank you very much! It's working now!

Dan Christensen
Guest
Posts: n/a

 07-11-2004
something like:

from random import random
from math import sqrt

def v1(n = 500000):
rand = random
sqr = sqrt
sm = 0
for i in xrange(n):
sm += sqr(1-rand()**2)
return 4*sm/n

This takes about 0.68s on a 2.66GHz P4 (using Python 2.4a) and
about 0.29s with psyco (using Python 2.3.4).

One of the great things about python is that it's possible (and fun)
to move loops into the C level, e.g. with this implementation:

def v2(n = 500000):
a = numarray.random_array.uniform(0, 1, n)
return 4 * numarray.sum(numarray.sqrt(1 - a**2)) / n

which clocks in at 0.16s with Python 2.4a. (psyco doesn't help.)

The problem with the numarray approach is that it consumes large
amounts of memory for no good reason. So, now that generator
expressions have arrived, the next thing to try is:

def v6(n = 500000):
rand = random
sqr = sqrt
sm = sum(sqr(1-rand()**2) for i in xrange(n))
return 4*sm/n

Unfortunately, this takes 0.64s, not much better than the most direct
method. (I don't have psyco installed for 2.4, so I haven't tried
that, but I wouldn't be surprised if it was worse than the direct
method with psyco.)

Maybe the reason is that the loop is still in python? I tried some
awkward ways of eliminating the for loop, e.g.

from itertools import *
from operator import pow, sub, add

def v4(n = 500000):
rand = random
pw = pow
sqr = sqrt
sm = sum(imap(lambda dummy: sqr(1-rand()**2), repeat(None, n)))
return 4*sm/n

def v3(n = 500000):
rand = random
pw = pow
sqr = sqrt
sm = sum(imap(sqr,
imap(sub, repeat(1),
imap(pw,
imap(lambda dummy: rand(), repeat(None, n)),
repeat(2)))))
return 4*sm/n

but they are both slow (0.88s and 0.97s respectively).

Is there a faster way to get Python to compute a sequence of repeated
operations like this and sum them up?

Is there a more direct way to make an iterator like this?

itrand = imap(lambda dummy: rand(), repeat(None, n))

Wouldn't it be great if one could write iterator code as concisely as
the numarray code, e.g.

sm = sum(sqrt(1-itrand**2))

i.e. if operations like +, - and ** were overloaded to handle
iterators producing numbers ("numiterators"?) Could a method like
this be fast?

Dan

Fernando Perez
Guest
Posts: n/a

 07-11-2004
Dan Christensen wrote:

> Is there a faster way to get Python to compute a sequence of repeated
> operations like this and sum them up?

Since the killer here is the dynamic type checks on the tight loops, until we
get optional static typing in python you'll need C for this. Fortunately,
weave makes this a breeze. Here's your v1() along with a trivially C-fied
version called v2:

################################################## #######################
import random
import math
import weave

def v1(n = 100000):
rand = random.random
sqrt = math.sqrt
sm = 0.0
for i in xrange(n):
sm += sqrt(1.0-rand()**2)
return 4.0*sm/n

def v2(n = 100000):
"""Implement v1 above using weave for the C call"""

support = "#include <stdlib.h>"

code = """
double sm;
float rnd;
srand(1); // seed random number generator
sm = 0.0;
for(int i=0;i<n;++i) {
rnd = rand()/(RAND_MAX+1.0);
sm += sqrt(1.0-rnd*rnd);
}
return_val = 4.0*sm/n;"""
return weave.inline(code,('n'),support_code=support)
################################################## #######################

Some timing info:

In [74]: run montecarlo_pi.py

In [75]: genutils.timings 1,v1
-------> genutils.timings(1,v1)
Out[75]: (1.4000000000000057, 1.4000000000000057)

In [76]: genutils.timings 10,v2
-------> genutils.timings(10,v2)
Out[76]: (0.13999999999998636, 0.013999999999998635)

In [77]: __[1]/_[1]
Out[77]: 100.00000000001016

A factor of 100 improvement is not bad for a trivial amount of work weave is
truly a cool piece of machinery. Just so you trust the values are OK:

In [80]: v1()
Out[80]: 3.1409358710711448

In [81]: v2()
Out[81]: 3.141602852608508

In [82]: abs(_-__)
Out[82]: 0.00066698153736322041

In [83]: 1.0/sqrt(100000)
Out[83]: 0.003162277660168379

The difference is within MonteCarlo tolerances (the usual 1/sqrt(N)).

Note that the above isn't 100% fair, since I'm calling the stdlib rand, a
C-coded simple generator, while v1 calls python's random.py, a Python-coded
Wichmann-Hill one. But still, for tight numerical loops this is pretty typical
of weave's results. Combined with the fact that via blitz, weave gives you
full access to Numeric arrays, it's a bit of a dream come true.

Cheers,

f

Paul Rubin
Guest
Posts: n/a

 07-11-2004
Fernando Perez <(E-Mail Removed)> writes:
> Note that the above isn't 100% fair, since I'm calling the stdlib
> rand, a C-coded simple generator, while v1 calls python's random.py,
> a Python-coded Wichmann-Hill one.

random.random in python 2.2 and later calls the Mersenne Twister generator,
a very fast generator implemented in C.

Fernando Perez
Guest
Posts: n/a

 07-11-2004
Paul Rubin wrote:

> Fernando Perez <(E-Mail Removed)> writes:
>> Note that the above isn't 100% fair, since I'm calling the stdlib
>> rand, a C-coded simple generator, while v1 calls python's random.py,
>> a Python-coded Wichmann-Hill one.

>
> random.random in python 2.2 and later calls the Mersenne Twister generator,
> a very fast generator implemented in C.

Sure?

[python]> ipython
Python 2.2.3 (#1, Oct 15 2003, 23:33:35)

IPython 0.6.1.cvs -- An enhanced Interactive Python.
? -> Introduction to IPython's features.
@magic -> Information about IPython's 'magic' @ functions.
help -> Python's own help system.
object? -> Details about 'object'. ?object also works, ?? prints more.

In [1]: import random

In [2]: random.random??
Type: instance method
Base Class: <type 'instance method'>
String Form: <bound method Random.random of <random.Random instance at
0xa166834>>
Namespace: Interactive
File: /usr/lib/python2.2/random.py
Definition: random.random(self)
Source:
def random(self):
"""Get the next random number in the range [0.0, 1.0)."""

# Wichman-Hill random number generator.
#
# Wichmann, B. A. & Hill, I. D. (1982)
# Algorithm AS 183:
# An efficient and portable pseudo-random number generator
# Applied Statistics 31 (1982) 188-190
#
# Correction to Algorithm AS 183
# Applied Statistics 33 (1984) 123
#
# McLeod, A. I. (1985)
# A remark on Algorithm AS 183
# Applied Statistics 34 (1985),198-200

# BEGIN CRITICAL SECTION
x, y, z = self._seed
x = (171 * x) % 30269
y = (172 * y) % 30307
z = (170 * z) % 30323
self._seed = x, y, z
# END CRITICAL SECTION

# Note: on a platform using IEEE-754 double arithmetic, this can
# never return 0.0 (asserted by Tim; proof too long for a comment).
return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0

In [3]:

This looks strangely like python to me

You are correct under 2.3:

[python]> python2.3 `which ipython`
Python 2.3.3 (#1, Mar 6 2004, 22:39:33)

IPython 0.6.1.cvs -- An enhanced Interactive Python.
? -> Introduction to IPython's features.
@magic -> Information about IPython's 'magic' @ functions.
help -> Python's own help system.
object? -> Details about 'object'. ?object also works, ?? prints more.

In [1]: import random

In [2]: random.random??
Type: builtin_function_or_method
Base Class: <type 'builtin_function_or_method'>
String Form: <built-in method random of Random object at 0x907614c>
Namespace: Interactive
Docstring:
random() -> x in the interval [0, 1).

In [3]:

The fact that ipython fails to show source under 2.3 is an indication that the
method is loaded from a C extension.

But the tests I posted where done using 2.2.

Cheers,

f

ps. My main point was to illustrate weave's usage, so this doesn't really
matter.

Paul Rubin
Guest
Posts: n/a

 07-11-2004
Fernando Perez <(E-Mail Removed)> writes:
> > random.random in python 2.2 and later calls the Mersenne Twister generator,
> > a very fast generator implemented in C.

>
> Sure?
> Python 2.2.3 (#1, Oct 15 2003, 23:33:35)
> Source:
> def random(self):
> """Get the next random number in the range [0.0, 1.0)."""
>
> # Wichman-Hill random number generator.

Hmm, I guess MT only became the default in 2.3. Thanks.