Sunday, June 08, 2008

Testing Random Number Generators with π

"Imagine that you have a dart board that is 2 units square. It inscribes a circle of unit radius. The center of the circle coincides with the center of the square. Now imagine that you throw darts at that dart board randomly. Then the ratio of the number of darts that fall within the circle to the total number of darts thrown is the same as the ratio of the area of the circle to the area of the square dart board. The area of a circle with unit radius is just π square unit. The area of the dart board is 4 square units. The ratio of the area of the circle to the area of the square is π / 4."

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:

helium said...
This comment has been removed by the author.
helium said...

You can use normal HTML entities:

< is &lt;

Andrew Dalke said...

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