Mailing List Archive: 49091 messages

# 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: [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/

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

> 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/