[REBOL] Re: Random numbers for Monte Carlo
From: lmecir:mbox:vol:cz at: 12-Aug-2002 16:15
Hi all,
<<Joel>>
...
here are two functions to compute PI via
Monte-Carlo methods, one using RANDOM and the other using
RANDOM/SECURE for equivalent purposes.
mcpi: func [
pts [integer!]
/local x y inside
][
inside: 0
loop pts [
x: 0.000001 * random 1000000
y: 0.000001 * random 1000000
if (x * x) + (y * y) <= 1.0 [inside: inside + 1]
]
print inside: inside / pts * 4
print pi
print rejoin ["" inside - pi * 1000000 "ppm"]
]
mcspi: func [
pts [integer!]
/local x y inside
][
inside: 0
loop pts [
x: 0.000001 * random/secure 1000000
y: 0.000001 * random/secure 1000000
if (x * x) + (y * y) <= 1.0 [inside: inside + 1]
]
print inside: inside / pts * 4
print pi
print rejoin ["" inside - pi * 1000000 "ppm"]
]
Both show the estimated value, the "true" value according to REBOL,
and the variance in parts-per-million.
Two QAD runs of each show a decided improvement from /SECURE...
>> mcpi 100000
3.14976
3.14159265358979
8167.346410207ppm
>> mcpi 100000
3.15064
3.14159265358979
9047.34641020699ppm
>> mcspi 100000
3.14172
3.14159265358979
127.34641020673ppm
>> mcspi 100000
3.13988
3.14159265358979
-1712.65358979333ppm
However, MCPI is about ten times faster than MCSPI, so you pay
for what you get...
HTH!
<</Joel>>
I computed the theoretical error estimate for the above Monte Carlo method
to see, whether Random/Secure is good enough and found out, that the
theoretical error estimate for Joel's test is 1298 ppm ((square-root 1 - (pi
/ 4) * pi / 4 / 100'000) * 1'000'000), which shows that Random/Secure worked
as well as possible.
Congratulations to RT!
-L
