Uniform Variate Generation Refs: Chapter 7 in Law, Pierre Lecuyer - PowerPoint PPT Presentation

Uniform Variate Generation Refs: Chapter 7 in Law, Pierre Lecuyer Tutorial, Winter Simulation Conference 2015 Peter J. Haas CS 590M: Simulation Spring Semester 2020 1 / 21

Pseudo-Random Numbers Overview Simple Congruential Generators Combined Generators Other Generators Testing Uniform Random Number Generators 2 / 21

Popular Mechanics 3 / 21

Pseudo-Random Numbers A deterministic sequence that “looks random” I Deterministic recurrence relation: sequence of integer seeds I Each seed converted into a uniform “random number” I A good generator has desirable theoretical properties, passes statistical tests Repeatability is a good thing I Facilitates debugging and verification I Allows “common random numbers” (e ffi ciency improvement) Versus “natural” sources of randomness I Ex: Silicon graphics / Cloudflare “lava lamp” generator I Ex: HotBits site uses radioactive decay I Ex: random.org uses atmospheric noise (radio static) I Non-reproducible and slow (OK for lotteries, games, generating encryption keys) 4 / 21

Pseudo-Random Numbers Overview Simple Congruential Generators Combined Generators Other Generators Testing Uniform Random Number Generators 5 / 21

Linear Congruential Generators (LCGs) Fundamental recurrence: x n +1 = ( ax n + c ) mod m I a = multiplier, c = increment, m = modulus I k mod m is remainder after dividing k by m (e.g., 14 mod 5 = 4 and 2 mod 10 = 2) I x n ’s take values in { 0 , 1 , 2 , . . . , m � 1 } I Return U n = x n / m I In C: rand() returns seed between 1 and RAND MAX I Historically the earliest (e ff ective) rng 6 / 21

Period of an LCG . . . . . . An LCG is periodic I Period  m I Want full period (every number in [0 .. m � 1] appears once) I Maximizes number of available random variates in a simulation I Otherwise, gaps may cause statistical anomalies 7 / 21

Multiplicative Congruential Generators (MCGs) , HE 1,43=3,44--1 - 3Xn mad 4 Anti Xo 't , N - , =3 - - . , Fundamental recurrence: x n +1 = ax n mod m I Special case of LCG with increment c equal to 0 I Full period: all values in { 1 , 2 , . . . , m � 1 } visited in cycle Theorem: An MCG has full period if m is prime and a is a primitive element modulo m I I.e., r ( a ) , min { k > 0 : a k mod m = 1 } = m � 1 is smallest integers 't . X not i' I Ex: x n +1 = 3 x n mod 4 - m - i Prime is evenly divisible by m ai . , I Ex: x n +1 = 3 x n mod 5 toil , Hi 3,1--4,45-2,44=1,45-3 k 1 2 3 4 - - - , 3 k mod 5 3 4 2 1 8 / 21

Classical MCGs Modulus choices I m = 2 b for convenience on binary computer I mod 2 b is simple: retain b lowest-order bits I Ex: IBM RANDU generator with m = 2 31 and a = 2 16 + 3 I For b > 3, period is at most m / 4 I m = 2 31 � 1 is, fortuitously, a (Mersenne) prime number I Because 2 31 � 1 is “almost” 2 31 , can compute mod quickly [Bratley et al., pp. 212–213] Lewis-Learmonth Generator (“Minimal Standard Generator”) x n +1 = 7 5 x n mod (2 31 � 1) = 6807 x n mod 2 , 147 , 483 , 647 Used for many years, but fails modern statistical tests, cycle is too short 9 / 21

Streams and Substreams in an MCG Jumping ahead I Goal: quickly compute x k for k large a k mod m I x k = ( a k x 0 ) mod m = �� � � x 0 mod m I Precompute numbers α k = a k mod m for multiple values of k I Allows partitioning of cycle into streams and substreams I Better than, e.g., setting y n = x 2 n and z n = x 2 n +1 I Caution: For an MCG, non-overlap is not su ffi cient (see demo) stream 1 4 m a e r t s s u b s t r e a m 1 s u b s t r e unigen.ResetNextSubstream() a m 2 unigen.ResetStartStream() unigen.ResetStartSubstream() stream 2 s t r e a m 3 10 / 21

Pitfalls of MCGs (and Other Generators) Short cycles I MCG numbers fall on a lattice I Only want to use O ( p period) numbers Low-order bits Claim: If x n +1 = ax n mod 2 b and r n +1 = x n +1 mod 2 k , where 3 < k < b , then ( r n : n � 1) has period at most 2 k � 2 n x n r n I r n ’s are k low-order bits 1 16049371 11 2 208641823 15 I Ex: x n +1 = 13 x n mod 2 31 with k = 4 3 564860051 3 [period = (2 31 / 4) � 1 ⇡ 537 ⇥ 10 6 ] 4 900729719 7 5 972068107 11 I So avoid algorithm that sets 6 1899467151 15 X = b nU c and V = nU � X 7 1070752835 3 where U D 8 1034884967 7 ⇠ Uniform(0 , 1) t.net 11 / 21

Other Pitfalls (Demo) Starting seeds for Lewis-Learmonth generator I X : use starting seed s = 1 I Y : use starting seed s 0 = 2 I s and s 0 are over 1.3 billion steps apart in cycle I Plot ( X , Y ) pairs I Plot histogram of X + Y Box-Muller and MCG 1. Generate U , V i.i.d. U [01 , ] 2. Set X = p� 2 log u cos(2 π V ) 3. Set Y = p� 2 log u sin(2 π V ) 4. Return X and Y independent N (0 , 1) 12 / 21

Pseudo-Random Numbers Overview Simple Congruential Generators Combined Generators Other Generators Testing Uniform Random Number Generators 13 / 21

Combined Generators Example: rngStream (used in Arena and other packages) x n = (1403580 x n − 2 � 810728 x n − 3 ) mod 4294967087 y n = (527612 y n − 1 � 1370589 y n − 3 ) mod 4294944443 z n = ( x n � y n ) mod 4294967087 u n = z n / 4294967087 I Seed = vector of six 32-bit integers Fun Fact I Cycle length ⇡ 2 191 ⇡ 10 57 (1 octodecillion) Time to mostly use I # streams = 2 64 ⇡ 10 19 (10 quintillion) up a generator with period of 2 191 with 1 I Stream length = 2 127 ⇡ 10 38 (100 undecillion) trillion computers I # substreams = 2 51 = 10 15 (1 trillion) generating one seed per nanosecond: I Substream length = 2 76 ⇡ 10 22 (10 sextillion) > 10 38 years I Well-behaved up to at least 45 dimensions 14 / 21

Pseudo-Random Numbers Overview Simple Congruential Generators Combined Generators Other Generators Testing Uniform Random Number Generators 15 / 21

Mersenne Twister General Form (Bit-wise Generator): X n = ( A 1 X n � 1 + · · · + A k X n � k ) mod 2 Y n = BX n mod 2 = ( y n , 1 , . . . , y n , W ) where W = 32 or 64 j =1 y n , j 2 � j or u n = 0 . y n , 1 y n , 2 · · · y n , W u n = P W I X n = ( x n , 1 , . . . , x n , L ) > with each x n , i = 0 or 1 I Binary matrices A 1 , . . . , A k ( L ⇥ L ) and B ( W ⇥ L ) I Fast: XOR operations and bit shifting Mersenne Twister I Default generator in Python and many other systems I Seed = vector of 623 integers (32 bit) I Period = 2 19937 � 1 ⇡ 10 6002 and well-behaved up to d = 623 I Drawbacks: large, slow initialization (demo), some stat. issues I WELL generator (Lecuyer et al.) scambles bits better 16 / 21

More Generators Cryptographic generators I “Impossible” to guess the next number in the sequence I Ex: RC4 (open source: Arc4Random), Threefish, ChaCha20 I Good for security, slow and statistically shaky for simulation Counter-based generators I Trivial seeds: s n = n for n � 1 (great for substreams!) I u n = f ( n ) where f is a “weak” but fast encryption function I Perform well, equidistribution properties not well understood Permuted congruential generator (PCG) I Use “improving” transformation of fast but shaky LCG I Under evaluation 17 / 21

Pseudo-Random Numbers Overview Simple Congruential Generators Combined Generators Other Generators Testing Uniform Random Number Generators 18 / 21

Testing Uniform Random Number Generators A simple generator: x n +1 = ( x n + 1) mod m x n + 1 = ( x n + 1 ) mod 31 Properties 1.0 0.8 I Full period 0.6 u n + 1 I Values uniformly spread out over [0 , 1] 0.4 I Yet: this is a terrible generator 0.2 0.0 Two ways of showing poor quality 0.0 0.2 0.4 0.6 0.8 1.0 u n I Compare to expected statistical behavior of uniform sequence 1 2 3 U 1 = m � 1 , U 2 = m � 1 , U 3 = m � 1 , . . . (not very random) I Look at possible values in higher dimensions (see plot) 19 / 21

x n + 1 = 3x n mod 31 Structural (Theoretical) Tests 1.0 0.8 Possible values in d dimensions 0.6 u n + 1 0.4 I Group the cycle into d -vectors: 0.2 V i = ( U i , . . . , U i + d � 1 ) 0.0 0.0 0.2 0.4 0.6 0.8 1.0 I Want equidistribution over [0 , 1] d u n x n + 1 = 3x n mod 31 1.0 Structural tests for MCGs 0.8 I Points of MCG lie on lattice: how even? 0.6 u n + 1 0.4 (spectral test) 0.2 I For modulus m , points lie on at most 0.0 ( d ! m ) 1 / d hyperplanes 0.0 0.2 0.4 0.6 0.8 1.0 u n d Upper bound I RANDU (demo) 2 31 1 2 16 2 For other generators: 3 2344 I Not a lattice; 32/64-bit U values can 4 476 appear multiple times in cycle 5 192 6 108 I Want same # of points in each “grid cell” 20 / 21

Statistical (Empirical) Tests How does generator jump from point to point? I Do U i numbers look i.i.d. uniform to a statistician? Many kinds of statistical tests I General tests for goodness of fit (e.g., χ 2 test) ( f j � n k ) 2 χ 2 = P k 1. Divide [0 , 1] into k ( > 100) equal intervals j =1 n / k 2. Generate U 1 , . . . , U n (where n ⇡ 10 k ) 3. Count number f 1 , . . . , f k that fall into each interval 4. Compute likelihood under i.i.d. uniform hypothesis I Serial test: essentially d -dimensional version of χ 2 test I Runs-up test (see homework) Test suites (the PRNG arms race) I Gold standard: TestU01 suite (incl. “SmallCrush”, “Crush”) [Lecuyer et al., simul.iro.umontreal.ca/testu01/tu01.html] 21 / 21