On 9/12/24 21:10, Krishna Myneni wrote:
...
I ran a simple Monte-Carlo integration of the area in a unit circle
using the 64-bit LCG PRNG (from random.4th in kForth dist) and
Marsaglia's 64-bit KISS PRNG. Hard to say which is better for this problem.
--
Krishna
=== Comparison of 64-bit LCG and KISS PRNGs ===
PRNG: RANDOM (random.4th)
Ntrials Area log(error)
-------------------------------
10^2 3.28 -0.86
10^3 3.18 -1.42
10^4 3.1724 -1.51
10^5 3.14048 -2.95
10^6 3.14350 -2.72
10^7 3.14225 -3.18
10^8 3.14152 -4.14
PRNG: RAN-KISS (kiss.4th)
Ntrials Area log(error)
-------------------------------
10^2 3.24 -1.01
10^3 3.10 -1.38
10^4 3.1252 -1.79
10^5 3.14084 -3.12
10^6 3.14073 -3.06
10^7 3.14184 -3.61
10^8 3.14167 -4.11
A bit more involved test of the same 64-bit PRNGs starts to show the possibility of defects in Marsaglia's 64-bit kiss prng (RAN-KISS)
compared to a simple 64-bit LCG prng (RANDOM). Instead of drawing from a uniform distribution, as assumed for these generators, the tests below
draw from a Maxwell-Boltzmann distribution, which describes the
probability distribution for speeds of ideal gas atoms at a given
temperature (non-interacting atoms except for hard sphere collisions in
three dimensions). The different "moments" of the speed may be computed
from the samples and compared with exact theoretical results: <v>,
<v^2>, <v^3>, etc., where v is a random sample of the speed for the M-B distribution at given temperature and mass of atom. The parameter "a",
given by
a = m/(2*kB*T)
where m is the mass, T is the temperature, and kB is the Boltzmann
constant, determines the shape of the M-B distribution. Unsigned 64-bit
random numbers, generated by the PRNG under test, are converted to
floating point numbers in the interval [0, 1]. Ideally, the PRNG will distribute these samples uniformly over the interval. The samples are
converted to random speed samples with an M-B distribution by an inverse mapping of the cumulative probability distribution for M-B to the speed.
The moments for <v>, <v^2>, and <v^3> are computed from these samples.
In this simulation, the temperature is set at 300 K, and the mass of the
He atom is used. Theory predictions for this m and T are
<v>_th = sqrt( 4/(pi*a) ) = 1258.727 m/s
<v^2>_th = 3/(2*a) = 1869538 (m/s)^2
<v^3>_th = sqrt( 16/(pi*a^3) ) = 3140144 x 10^3 (m/s)^3
Below are the results for the two PRNGs. For N trials of 10^5 through
10^7, the simple 64-bit LCG prng consistently gives less relative error
than Marsaglia's 64-bit kiss prng with respect to the theory predictions.
Source code for these test results is given below.
--
Krishna
=== begin tests ===
Test of prng RANDOM 64-bit (random.4th):
' random test-prng
Moments of speed
N <v> (m/s) <v^2> (m/s)^2 <v^3> (m/s)^3
10^2 1220.2444 1761665.1 2879481740.
10^3 1258.2315 1844889.5 3025976380.
10^4 1250.0696 1837249.9 3055346856.
10^5 1259.9899 1870367.4 3143506292.
10^6 1259.3923 1868383.4 3136761299.
10^7 1259.6343 1869366.9 3140010276.
Relative Errors
N |e1| |e2| |e3|
10^2 3.13e-02 5.77e-02 8.30e-02
10^3 1.19e-03 1.32e-02 3.64e-02
10^4 7.67e-03 1.73e-02 2.70e-02
10^5 2.08e-04 4.44e-04 1.07e-03
10^6 2.66e-04 6.18e-04 1.08e-03
10^7 7.39e-05 9.15e-05 4.27e-05
Test of prng RAN-KISS 64-bit (kiss.4th):
' ran-kiss test-prng
Moments of speed
N <v> (m/s) <v^2> (m/s)^2 <v^3> (m/s)^3
10^2 1212.8900 1702106.0 2666594142.
10^3 1274.7097 1900519.4 3190481667.
10^4 1259.6892 1866276.3 3123781859.
10^5 1260.4578 1872500.5 3147907366.
10^6 1260.2589 1871249.6 3145073404.
10^7 1259.8296 1869882.9 3140954422.
Relative Errors
N |e1| |e2| |e3|
10^2 3.72e-02 8.96e-02 1.51e-01
10^3 1.19e-02 1.66e-02 1.60e-02
10^4 3.03e-05 1.74e-03 5.21e-03
10^5 5.80e-04 1.58e-03 2.47e-03
10^6 4.22e-04 9.16e-04 1.57e-03
10^7 8.11e-05 1.85e-04 2.58e-04
=== end tests ===
=== begin test-prng-mb.4th ===
\ test-prng-mb.4th
\
\ Test a pseudo-random number generator by computing the
\ moments of powers of speed using random samples from a
\ Maxwell-Boltzmann distribution of ideal gas speeds in 3D.
\
\ Copyright 2024 Krishna Myneni
\
\ Permission is granted to use this code for any purpose,
\ provided the original source is referenced.
\
\ Usage Example:
\
\ ' random TEST-PRNG
\
\ where RANDOM is a user-provided PRNG which returns
\ a cell-sized pseudo-random number.
\
\ Notes
\
\ 1. The xt for any PRNG may be passed to the word
\ TEST-PRNG provided it returns a cell-sized
\ bit pattern. All bits of the cell are assumed
\ to be part of the pseudo-random number.
\
\ 2. The default M-B distribution is computed for
\ helium atoms at 300K. If you change the default
\ temperature or atomic mass, the maximum velocity
\ used in the root finder call in CDF^-1 may need
\ to be changed.
\
include ans-words
include strings
include phyconsts
include fsl/fsl-util
include fsl/erf
include fsl/regfalsi
include random
include kiss
[UNDEFINED] FPICK [IF]
cr cr .( This program requires a separate FP stack! ) cr
quit
[THEN]
4.0e0 pi f* fconstant 4pi
-1 0 d>f fconstant F_MAX_U
300.0e0 fconstant DEF_TEMPERATURE \ K
1.0e4 fconstant DEF_VMAX \ m/s
1.0e-9 fconstant DEF_TOL \ root-finder tolerance
4.002602e0 fconstant m_He \ mass in amu for He atom
: fcube ( F: r -- r^3 ) fdup fsquare f* ;
: f3dup 2 fpick 2 fpick 2 fpick ;
[UNDEFINED] FTUCK [IF]
: ftuck ( F: a b -- b a b ) fswap fover ;
[THEN]
\ Relative error between a value and its reference value
: rel-error ( F: a aref -- e ) ftuck f- fswap f/ ;
\ The random number generator to be tested must return a
\ cell-sized bit pattern (unsigned single) -- all bits of
\ the cell are assumed to be randomized.
defer prngU
defer prngU-init
' random is prngU \ default PRNG
:noname ( -- ) 1000 0 do prngU drop loop ; is prngU-init
\ Generate a floating point random number between 0 and 1
\ using the selected unsigned single cell prng.
: randf ( F: -- r ) prngU 0 d>f F_MAX_U f/ ;
fvariable m_atom_kg
: set-atomic-mass ( F: m_amu -- ) kg/amu f* m_atom_kg f! ;
m_He set-atomic-mass \ default value of mass of atom
fvariable a
: set-T ( F: Tk -- )
kB f* 2.0e0 f* m_atom_kg f@ fswap f/ a f! ;
\ Return probability density, p(v), for the M-B distribution
\ with parameter "a"
: pdf ( F: v -- p )
fsquare fdup a f@ f* fnegate fexp f*
4pi f* a f@ pi f/ 1.5e0 f** f* ;
\ Return cumulative probability, C(v1) i.e. P(0 <= v < v1),
\ the integral of p(v) from 0 to v1 for the M-B distribution.
: cdf ( F: v1 -- P )
fdup fsquare a f@ f* fdup fsqrt erf
fswap fnegate fexp frot f* 2e f* a f@ pi f/ fsqrt f* f- ;
\ Invert cdf to obtain velocity from value of cdf
fvariable pv1
: cdf^-1 ( F: P -- v )
pv1 f!
[: cdf pv1 f@ f- ;] 0.0e0 DEF_VMAX DEF_TOL )root ;
\ N random trials of velocities following the M-B distribution
\ are sampled and the moments v, v^2, and v^3 are computed
\ from the samples and compared to theory.
fvariable vsum
fvariable v2sum
fvariable v3sum
: mb-prng ( Ntrials -- ) ( F: -- <v> <v^2> <v^3> )
DEF_TEMPERATURE set-T
0.0e0 fdup fdup vsum f! v2sum f! v3sum f!
dup
0 ?DO
randf \ random draw of cdf value
cdf^-1 \ invert cdf to get random v in M-B distribution
fdup vsum f+!
fdup fsquare v2sum f+!
fcube v3sum f+!
LOOP
s>f
vsum f@ fover f/ fswap
v2sum f@ fover f/ fswap
v3sum f@ fswap f/ ;
\ Theoretical moments of v for the M-B distribution at given
\ temperature and mass, specified by parameter "a".
: <v>_th ( F: -- <v> ) 4.0e0 pi a f@ f* f/ fsqrt ;
: <v^2>_th ( F: -- <v^2> ) 3.0e0 2.0e0 a f@ f* f/ ;
: <v^3>_th ( F: -- <v^3> ) 16.0e0 pi a f@ fcube f* f/ fsqrt ;
\ Compute relative errors of randomly sampled averages with
\ respect to the theoretical moments.
: moment-errors ( F: <v> <v^2> <v^3> -- e1 e2 e3 )
frot <v>_th rel-error
frot <v^2>_th rel-error
frot <v^3>_th rel-error ;
9 constant MAX-POW
MAX-POW 1+ 3 float matrix v_mom{{
MAX-POW 1+ 3 float matrix m_rel{{
: print-tables ( -- )
cr ." Moments of speed"
cr ." N <v> (m/s) <v^2> (m/s)^2 <v^3> (m/s)^3" cr
MAX-POW 1- 2 DO
v_mom{{ I 2 }} f@ v_mom{{ I 1 }} f@ v_mom{{ I 0 }} f@
." 10^" I 1 .r 2 spaces
12 4 f.rd 2 spaces
12 1 f.rd 6 spaces
12 0 f.rd cr
LOOP
cr ." Relative Errors"
cr ." N |e1| |e2| |e3|" cr
precision \ save precision
3 set-precision
MAX-POW 1- 2 DO
m_rel{{ I 2 }} f@ m_rel{{ I 1 }} f@ m_rel{{ I 0 }} f@
." 10^" I 1 .r 2 spaces
fabs fs. 2 spaces
fabs fs. 2 spaces
fabs fs. cr
LOOP
set-precision \ restore precision
cr
;
: test-prng ( xt -- )
is prngU
prngU-init
MAX-POW 1- 2 DO
I s>f falog fround>s
mb-prng
f3dup
v_mom{{ I 2 }} f! v_mom{{ I 1 }} f! v_mom{{ I 0 }} f!
moment-errors
m_rel{{ I 2 }} f! m_rel{{ I 1 }} f! m_rel{{ I 0 }} f!
LOOP
print-tables
;
cr cr .( Usage: ' <prng> test-prng )
cr cr .( e.g. ' random test-prng ) cr cr
=== end test-prng-mb.4th ===
--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)