Mailing List Archive: 49091 messages

[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
```