Velocity Reviews - Computer Hardware Reviews

Velocity Reviews > Newsgroups > Programming > Python > Monte Carlo Method and pi

Reply
Thread Tools

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.

Please help me!

Regards
Karl


 
Reply With Quote
 
 
 
 
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


 
Reply With Quote
 
 
 
 
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-----

 
Reply With Quote
 
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.
 
Reply With Quote
 
Karl Pech
Guest
Posts: n/a
 
      07-08-2004
Thank you very much! It's working now!

 
Reply With Quote
 
Dan Christensen
Guest
Posts: n/a
 
      07-11-2004
I got to thinking about this problem. The most direct implementation is
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
 
Reply With Quote
 
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
 
Reply With Quote
 
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.
 
Reply With Quote
 
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)
Type "copyright", "credits" or "license" for more information.

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
#
# see also:
# 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

# This part is thread-unsafe:
# 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)
Type "copyright", "credits" or "license" for more information.

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.
 
Reply With Quote
 
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.
 
Reply With Quote
 
 
 
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are Off


Similar Threads
Thread Thread Starter Forum Replies Last Post
Shariffa Carlo xoxoxo xoxoxo C++ 0 12-07-2009 08:24 PM
Monte Carlo Simulation raja Python 0 08-12-2008 09:29 AM
DVD Verdict reviews: ORCHESTRA WIVES, GANKUTSUOU: THE COUNT OF MONTE CRISTO, CHAPTER 3, and more! DVD Verdict DVD Video 0 03-30-2006 09:19 AM
A simple program with Monte Carlo Aleramo C Programming 3 12-04-2005 06:03 PM
random number generation and Monte Carlo simulation in C++ mescaline C++ 4 09-10-2003 09:01 PM



Advertisments