I got the book Digital dice the other day

It is about using monte carlo methods to solve probability problems that are hard to solve analytically. The first example given is how to estimate pi based on the % of darts that fall inside a cirle that is inside a square with the side length of twice the circles radius.

I decided to try out some languages to see how they got on with this problem. The aim is to test

1. Accuracy for mante carlo simulations

2. Ease of writing

3. Speed

in order of importance

If a language cannot accurately estimate π using this method then it indicates its random number generator (RNG) might not be very good. To test the RNG propery you need to do chi square tests and such as described in Seminumerical Algorithms by Knuth but i think its funny you can use π to test a RNG.

Perl

$N=10000000;

$p=0;

$k=1;

while ($k<$N){

$x=rand();

$y=rand();

if (($x*$x)+($y*$y)<1){

$p=$p+1;

}

$k=$k+1;

}

$ans=((4*$p)/$N);

print "ans".$ans;

Ruby

N=10000000

p=0

k=1

while k < N

x=rand

y=rand

if (x*x)+(y*y)<1

p=p+1

end

k=k+1

end

ans=(4*p).quo(N)

puts "pi is #{ans}"

Scheme

I cannot figure out how to get random numbers with scheme. There is code to do it here but I am not going to write my own rand function.

C

int main(void)

{

long N=10000000;

double x=0;

double y=0;

double p=0;

int k=1;

while (k < N){

x=( (double)rand() / ((double)(RAND_MAX)+(double)(1)) );

y=( (double)rand() / ((double)(RAND_MAX)+(double)(1)) );

if ((x*x)+(y*y)<1){

p=p+1;

}

k=k+1;

}

double mypi=((4*p)/N);

printf ("my pi %f \n",mypi);

return 0;

}

(less then sign) is < (thanks Helium).

I should use srand here to make it reproducible. Or more to be really fair run the test a few times and collect the average error. This is how the Digital dice runs the test with matlab. Anyway I wont tell you what language worked out the best. If you have any suggestions or code please comment.

## 3 comments:

You can use normal HTML entities:

< is <

Don't use rand(3C). It's not random enough for anything more than the proverbial "does the monster attack?" type decision.

Slightly better is drand48, but that is still a linear congruential generator, with well known problems. See for example http://random.mat.sbg.ac.at/~charly/server/server.html .

You want to use the Mersenne twister algorithm. It was designed with Monte Carlo uses in mind. Note that it is *not* a good algorithm for crypto.

http://en.wikipedia.org/wiki/Mersenne_twister

Post a Comment