[REBOL] Re: Random numbers for Monte Carlo
From: joel:neely:fedex at: 10-Aug-2002 10:26
Hi, Ladislav,
Ladislav Mecir wrote:
> do you know how much suitable are Rebol random numbers for a
> Monte Carlo simulation? (I need more than one million reasonably
> random numbers). Would you suggest to use /SECURE for this?
>
It depends on how your simulation uses the pseudo-random values.
(It would also be nice to know what algorithm RANDOM uses, but I
haven't had the time to do the detective work, and don't recall
seeing that described. I'd be grateful for any corrections to
my present ignorance on that score.)
It also depends on your definition of "reasonable". At the end
of this note is a little off-the-cuff script that will allow you
to see how much divergence from pure uniformity RANDOM produces.
On a more practical level, the best advice I've ever heard on
using Monte-Carlo methods goes something like this:
1) If your problem has an analytic solution, DON'T use
Monte-Carlo simulation, as it will only be worse.
2) If you don't have an analytic solution, find a simpler
version of your problem (one that uses the same kind
of sampling) which DOES have an analytic solution and
test your Monte-Carlo approach on that problem to see
how well it does. Only use the M-C simulation on the
full problem if you can satisfy yourself that the
solution for the reduced problem is good enough, and
that the two programs will REALLY be doing comparable
calculations.
In that spirit, 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!
-jn-
rnn: func [
n [integer!] samples [integer!]
/local stats prev curr pure fmt fmtlen dev hdev ldev t
][
samples: n * n * pure: to-integer samples / n / n
fmtlen: 1 + length? join " " pure
fmt: func [x [integer!] /local r] [
head insert/dup r: join copy "" x " " fmtlen - length? r
]
stats: copy []
loop n [
append/only stats copy/deep array/initial n [[0]]
]
curr: random n
loop samples [
prev: curr
curr: random n
stats/:prev/:curr/1: 1 + stats/:prev/:curr/1
]
print ["Samples:" samples ", expecting:" pure]
hdev: ldev: 0
foreach row stats [
t: 0
foreach val row [
prin [" " fmt dev: val/1 - pure]
t: t + dev
hdev: max hdev dev
ldev: min ldev dev
]
print [": " fmt t]
]
print ["Deviations from" ldev "to" hdev]
]
It checks for consecutive correlation in the simplest form
possible, but the output could be used for chi-squared analysis
if you wish.
>> rnn 8 800000
Samples: 800000 , expecting: 12500
8 25 10 -1 -8 6 2 -17: 25
-7 13 10 -23 -19 21 6 17: 18
22 1 -11 5 -8 -5 -16 12: 0
7 -4 -1 28 -2 -24 -8 14: 10
13 -3 13 -14 -2 -16 6 -19: -22
-12 12 0 12 13 -1 10 -25: 9
-16 -17 -1 9 4 21 28 -19: 9
10 -9 -21 -6 0 8 -19 -12: -49
Deviations from -25 to 28
--
; Joel Neely joeldotneelyatfedexdotcom
REBOL [] do [ do func [s] [ foreach [a b] s [prin b] ] sort/skip
do function [s] [t] [ t: "" foreach [a b] s [repend t [b a]] t ] {
| e s m!zauafBpcvekexEohthjJakwLrngohOqrlryRnsctdtiub} 2 ]
