**[Al--Bri--xtra--co--nz]**Web: http://valley.150m.com/ Date: 10/August/2002 ] Random-Numbers: function [ "Read a bunch of random numbers from random.org." N [integer!] "The number of random numbers to read." Minimum [integer!] "The minimum value of each number." Maximum [integer!] "The maximum value of each number." ] [Digit Numbers Number] [ Digit: charset "0123456789" if all [ 1 <= N N <= 10000 -1000000000 <= Minimum Minimum <= 1000000000 -1000000000 <= Maximum Maximum <= 1000000000 ] [ Numbers: make block! N foreach Number parse read rejoin [ http://www.random.org/cgi-bin/randnum #"?" "num=" N #"&" "min=" Minimum #"&" "max=" Maximum #"&" "col=1" ] "" [ append Numbers to-integer Number ] ] Numbers ] probe Random-Numbers 5 1 100 Andrew Martin Random Rebolutionary... ICQ: 26227169 http://valley.150m.com/

# Random numbers for Monte Carlo

### [1/9] from: lmecir::mbox::vol::cz at: 10-Aug-2002 10:57

Hi all,
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?
TIA
Ladislav

### [2/9] from: al:bri:xtra at: 10-Aug-2002 22:42

Sunanda wrote:

> No one has yet contributed Rebol code to random.org for pulling down

blocks of numbers. Anyone up for it?
Try this:
Rebol [
Name: 'Random-Numbers
Title: "Random Numbers"
File: %"Random Numbers.r"
Author: "Andrew Martin"
eMail:

### [3/9] from: lmecir:mbox:vol:cz at: 10-Aug-2002 17:37

Thanks, Sunanda and Andrew,
you have been extremely helpful.
-L

### [4/9] 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 ]

### [5/9] from: lmecir:mbox:vol:cz at: 10-Aug-2002 19:24

Thanks, Joel,
your advice is wise. Some subproblems have got an analytical solution, so I
can test the accuracy of my simulation, but the final solution isn't
obtainable (at least for me) analytically. The speed problem is present,
because a million number crunching in the presence of nonlinearities is
really a hard job.
For now I decided to use the fastest generator, after some debugging I may
use a better method.
Ladislav

### [6/9] from: tomc:darkwing:uoregon at: 10-Aug-2002 10:46

Hi Ladislav,
When playing with Monte Carlo for a class a few years back I was content
to periodically reseed Rebols RNG with a very random number.
(when I just now checked, the HotBits site was up but not serving numbers)
comment{
seed the RNG with a radioactive decay induced non-pesudo-random number if
connected to the net(and HotBits server up) or with the current time
otherwise
}
random/seed either
all [ not error? try[
rnd-pg: read
http://www.fourmilab.ch/cgi-bin/uncgi/Hotbits/?fmt=hex&nbytes=2]
parse rnd-pg [thru <pre> "^/" copy seed to "^/"]
]
[make integer! make binary! seed]
[now]
On Sat, 10 Aug 2002, Ladislav Mecir wrote:

### [7/9] from: al:bri:xtra at: 11-Aug-2002 11:44

Sunanda wrote:

> Maybe needs a touch of error trapping when minimum and maximum are

reversed:
Fixed!
Rebol [
Name: 'Random-Numbers
Title: "Random Numbers"
File: %"Random Numbers.r"
Author: "Andrew Martin"
eMail:

**[Al--Bri--xtra--co--nz]**Web: http://valley.150m.com/ Date: 11/August/2002 Usage: [ probe Random-Numbers 5 1 100 probe random-numbers 3 10 5 probe random-numbers 3 -5 -10 ] ] Random-Numbers: function [ "Read a bunch of random numbers from random.org." N [integer!] "The number of random numbers to read." Minimum [integer!] "The minimum value of each number." Maximum [integer!] "The maximum value of each number." ] [Numbers] [ if all [ 1 <= N N <= 10000 -1000000000 <= Minimum Minimum <= 1000000000 -1000000000 <= Maximum Maximum <= 1000000000 ] [ Numbers: load rejoin [ http://www.random.org/cgi-bin/randnum #"?" "num=" N #"&" "min=" min Minimum Maximum #"&" "max=" max Minimum Maximum #"&" "col=1" ] ] Numbers ] And I've sent it off to http://www.random.org as well. Andrew Martin ICQ: 26227169 http://valley.150m.com/### [8/9] 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

### [9/9] from: al:bri:xtra at: 13-Aug-2002 16:00

Sunanda wrote:

> Other than that, you should submit it to them,

:)
Hi Andrew,
Thanks for the client. I've put it on the random.org client web page.
All the best,
--
Mads Haahr <

**[Mads--Haahr--cs--tcd--ie]**> | Lecturer in Computer Science Department of Computer Science | Phone: +353 1 608 1540 Trinity College Dublin | Web: http://www.cs.tcd.ie/ Andrew Martin ICQ: 26227169 http://valley.150m.com/