## [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 1000003.14976 3.14159265358979 8167.346410207ppm>> mcpi 1000003.15064 3.14159265358979 9047.34641020699ppm>> mcspi 1000003.14172 3.14159265358979 127.34641020673ppm>> mcspi 1000003.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 800000Samples: 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 ]