Random Number and Random Variate Generation

1. Random Number and Random Variate Generation

2. Random Number and Random Variate Generation • How does the computer generate observations from various distributions specified after input analysis? • There are two main components to the generation of observations from probability distributions. • Random number generation. • Random variate generation.

3. Random Number Generation • Random number generation – The generation of U(0,1) random variates (observations from Uniform (0,1) distribution). • This serves as the foundation for the generation of observations from other distributions, which is called random variate generation. • Random Number Generator is the term used to describe the procedure and parameters used to generate the U(0,1) observations.

4. Random Number Generation • Since the “stream” of random numbers generated is reproducible, random number generation procedures are also referred to as pseudo random number generators. • The stream or sequence of numbers produced by a generator should pass statistical tests for randomness. • An outside observer should not be able to tell the difference (statistically) between a stream of pseudo random numbers and an actual random number stream.

5. Pseudorandomness • A pseudorandom process appears random, but isn’t • Pseudorandom sequences exhibit statistical randomness • but generated by a deterministic process • Pseudorandom sequences are easier to produce than a genuine random sequences • Pseudorandom sequences can reproduce exactly the same numbers • useful for testing and fixing software.

6. Pseudorandom Generator Seeds • Random number generators typically compute the next number in the sequence from the previous number • The first number in a sequence is called the seed • to get a new sequence, supply a new seed(current machine time is useful) • to repeat a sequence, repeat the seed

7. Random Number Generation • Desirable Attributes: • Uniformity • Independence • Efficiency • Replicability • Long Cycle Length

8. Random Number Generation (cont.) Each random number Rt is an independent sample drawn from a continuous uniform distribution between 0 and 1 ì1 , 0 £ x £ 1 pdf: f(x) = í î0 , otherwise

9. 1 PDF: f(x) x 0 Random Number Generation (cont.)

10. Random Number Generation • One early method – the midsquare method(von Neumann and Metropolis 1940) • Start with a four digit positive integer Z0. • Square Z0 to get an integer with up to eight digits (append zeros if less than eight). • Take the middle four digits as the next four digit integer Z1. • Place a decimal point to the left of Z1 to form the first “U(0,1)” observation. • Repeat

11. MidSquare • Example: • X0 = 7182 (seed) • = 51581124 • ==> R1= 0.5811 • = (5811) 2 = 33767721 • ==> R2= 0.7677 • etc. Random Number Generation - MidSquare

12. Note: Cannot choose a seed that guarantees that the sequence will not degenerate and will have a long period. Also, zeros, once they appear, are carried in subsequent numbers. • Ex1: X0 = 5197 (seed) = 27008809 • ==> R1 = 0.0088 = 00007744 • ==> R2 = 0.0077 • Ex2: X0 = 4500 (seed) = 20250000 • ==> R1 = 0.2500 = 06250000 • ==> R2 = 0.2500 Random Number Generation- MidSquare

13. Random Number Generation • The prior method does not work well. • Degenerates to zero. • What are good methods? • Linear Congruential Generators (LCGs). • Composite generators. • Tausworthe generators.

14. Random Number Generation - LCG • Linear Congruential Generators (LCGs). • A LCG generates a sequence of integers Z1, Z2, Z3,… using the following recursive formula, • mod m is short for modulo m or the remainder when divided by m.

15. Random Number Generation - LCG • Since the mod m operation is used, all Zi’s will be between 0 and m-1. • To get the “U(0,1)” random observations each Zigenerated is divided by m. • So are the Ui’s really U(0,1) random observations ?

16. In-class Exercise • Let m=63, a=22, c=4 and Z0 =19. • Generate the first five “U(0,1)” observations.

17. Random Number Generation

18. Random Number Generation - LCG • What will happen after the 63rd number is generated? • m, a, and c are the parameters of the random number generator. • There can be an infinite number of different implementations of a LCG. • The values used for m, a, and c determine whether the generator is good or bad.

19. Random Number Generation - LCG • The example LCG demonstrates cycling in the prior table. • Since m=63, it can generate at most 63 numbers before it repeats the same sequence. • This small random number generator has full period since it generates all possible (m=63) numbers before cycling. • A long period (full if possible) is desirable since more observations can be generated before cycling. • No “gaps”.

20. Random Number Generation • The example generator has full period but bad statistical properties (next slide). • A good random number generator will have values for m, a, and c such that full or close to full period is obtained, as well as good statistical properties. • Crystal Ball • m = 231 – 1 • a = 62089911 • c = 0 • Period = 231 – 2

21. Random Number Generation - LCG • Theorem (Hull and Dobell 1962) • The LCG Zi = (aZi-1 + c) mod m has full period if and only if the following three conditions hold. • The only positive integer that exactly divides both m and c is 1. • If q is a prime number that divides m, then q divides a-1. • If 4 divides m, then 4 divides a-1. • The parameters of the LCG dictate the period length of the LCG as well as other properties of the numbers generated.

22. Period vs. Autocorrelation • A generator that has the maximum possible period is called a full-period generator. • Lower autocorrelations between successive numbers are preferable. • Both generators have the same full period, but the first one has a correlation of 0.25 between xn-1 and xn, whereas the second one has a negligible correlation of less than 2-18.

23. Random Number Generation • Types of LCGs • When c = 0, the LCG is called a multiplicative generator. • When c ≠ 0, the LCG is called a mixed generator. • Most LCGs implemented are multiplicative • Can’t have full period. • How is m selected. • A large period is desired → m=231 (based on a 32 bit word size). • With m=231 it has been proven that the period can be at most 229 (25% of the values are cycled and gaps may be present).

24. Multiplicative LCG • Multiplicative LCG: c=0 • Two types:

25. Multiplicative LCG Example: Using the multiplicative congruential method, find the period of the generator for a = 13, m = 26, and X0 = 1, 2, 3, and 4. The solution is given in next slide. When the seed is 1 and 3, the sequence has period 16. However, a period of length eight is achieved when the seed is 2 and a period of length four occurs when the seed is 4.

26. Period Determination Using Various seeds i Xi XiXiXi 0 1 2 3 4 1 13 26 39 52 2 41 18 59 36 3 21 42 63 20 4 17 34 51 4 5 29 58 23 6 57 50 43 7 37 10 47 8 33 2 35 9 45 7 10 9 27 11 53 31 12 49 19 13 61 55 14 25 11 15 5 15 16 1 3 Multiplicative LCG

27. Multiplicative LCG with m=2k • Maximum possible period 2k-2 • Period achieved if multiplier a is of the form 8i± 3, and the initial seed is an odd integer • One-fourth the maximum possible may not be too small • Low order bits of random numbers obtained using multiplicative LCG's with m=2k have a cyclic pattern.

28. Examining Bits of a Multiplicative LCG

29. Examining Bits of a Mixed LCG

30. Multiplicative LCG with m≠2k • When the modulus m is a prime number and a >1, the maximum period is m-1, no matter whether c=0 or not. • The maximum period m-1 is obtained if and only if ‘a’ is a primitive element modulo m. • If m is prime then ‘a’ is a primitive element modulo m (or primitive root of m) if and only ifan mod m ≠ 1 for n=1, 2, 3, …,m-2. • Recommended. (Prime moduli are best in terms of sequence randomness.)

31. Multiplicative LCG with m≠2k Example: • Starting with a seed of x0=1: 1, 3, 9, 27, 19, 26, 16, 17, 20, 29, 25, 13, 8, 24, 10, 30, 28, 22, 4, 12, 5, 15, 14, 11, 2, 6, 18, 23, 7, 21, 1, … The period is 30 ⇒ 3 is a primitive root of 31 • With a multiplier of a = 5: 1, 5, 25, 1,… The period is only 3 ⇒ 5 is not a primitive root of 31 • Primitive roots of 31= 3, …????????.

32. Random Number Generation • LCGs are a special case of the form Zi = g(Zi-1, Zi-2, ...) (mod m), Ui = Zi/m, for some function g • Examples: • g(Zi-1) = aZi-1 + c LCG • g(Zi-1, Zi-2, ..., Zi-q) = a1Zi-1 + a2Zi-2 + ... + aqZi-q multiple recursive generator • g(Zi-1) = a'Z2i-1 + aZi-1 + c quadratic CG • g(Zi-1, Zi-2) = Zi-1 + Zi-2 Fibonacci (bad)

33. Random Number Generation Composite Generators Combine two (or more) individual generators in some way. • Differencing LCGs • Z1i and Z2i from LCGs with different moduli • Let Zi = (Z1i – Z2i ) (mod m); Ui = Zi / m • Very good statistical properties • Very portable (micros, different languages) • Wichmann/Hill • Use three LCGs to get U1i, U2i, and U3i sequences • Let Ui = fractional part of U1i + U2i + U3i • Long period, good statistics, portability

34. Tausworthe Generator • Originated in cryptography • Can achieve very long periods • Theoretical appeal: for properly chosen parameters, can prove that over a cycle, • mean 1/2 (as for true U(0,1)) • Variance  1/12 (as for true U(0,1)) • Autocorrelation  0 (as for true IID sequence) • Define a sequence of binary digits B1,B2, . . ., by where cj = 0 or 1.

35. Tausworthe Generator • Looks a bit like a generalization of LCG’s. • Let D = delay operator such that Db(n)=b(n+1) • or • Since in mod 2 arithmetic subtraction is equivalent to addition, the preceding equation is equivalent to

36. Tausworthe Generator • The polynomial on the left-hand side of this equation is called a characteristic polynomial and is traditionally written using x in place of D • The period of a Tausworthe generator depends upon the characteristic polynomial. In particular, the period is the smallest positive integer n for which xn - 1 is divisible by the characteristic polynomial. The maximum possible period with a polynomial of order q is 2q - 1. The polynomials that give this period are called primitive polynomials.

37. Tausworthe Generator • Example:Consider the following polynomial: x7 + x3 + 1 Using the D operator in place of x, we get or or using the XOR operator or Substituting n-7 for n, we get

38. Tausworthe Generator Starting with b0 = b1 = ... = b6 = 1, we get the following bit sequence:

39. Tausworthe Generator • The complete sequence is: 1111111 0000111 0111100 1011001 0010000 0010001 0011000 1011101 0110110 0000110 0110101 0011100 1111011 0100001 0101011 1110100 1010001 1011100 0111111 1000011 1000000. • Period = 127 or 27-1 bits ⇒ The polynomial x7+x3+1 is a primitive polynomial.

40. Tausworthe Generator • A Tausworthe sequence can be easily generated in hardware using Linear-Feedback Shift Registers (LFSRs). • For example, the polynomial x5 + x3 + 1 results in the generator bn = bn-2bn-5. This can be implemented using the LFSR shown in the Figure presented next slide. • The circuit consists of six registers, each holding one bit. On every clock cycle, each register’s content is shifted out, and the new content is determined by the input to the register.

41. Tausworthe Generator • Linear Feedback Shift Register: x5+x3+1 ⇒ bn= bn-2 ⊕ bn-5 This can be easily implemented using shift registers:

42. Tausworthe Generator Generating U(0,1): • Divide the sequence into successive groups of s bits and use the first l bits of each group as a binary fraction: xn = 0.bsnbsn+1bsn+2bsn+3 ...bsn+l-1 • Here, sis a constant greater than or equal to l and is relativelyprime to 2q-1. • s ≥ l ⇒ xn and xj for n ≠ j have no bits in common. • Relative prime-ness guarantees a full period 2q-1 for xn.

43. Tausworthe Generator Example:bn = bn-4 ⊕ bn-7 • The period 27-1=127 • l=8, s=8: x0 = 0.111111102 = 0.9921910 x1 = 0.000111012 = 0.1132810 x2 = 0.111001012 = 0.8945310 x3 = 0.100100102 = 0.2968810 x4 = 0.000001002 = 0.3632810 x5 = 0.010011002 = 0.4218810 ….

44. Tausworthe Generator List of Primitive Trinomials x2 + x + 1 x3+ x + 1 x4 + x + 1 x5+ x2 + 1 x6 + x + 1 x7+ x + 1 x7+ x3 + 1 x9+ x4 + 1 x10 + x3 + 1 x11+ x2 + 1 x15+ x + 1 x15+ x4 + 1 x15 + x7 + 1 x17+ x3 + 1 x17+ x5 + 1 x17+ x6 + 1 x18 + x7 + 1 x20+ x3 + 1 x21+ x2 + 1 x22+ x + 1 x23 + x5 + 1 x23+ x9 + 1 x25+ x3 + 1 x25+ x7 + 1 x28 + x3 + 1 x28+ x9 + 1 x28+ x13 + 1 x29+ x2 + 1 x31 + x3 + 1 x31+ x6 + 1 x31+ x7 + 1 x31+ x13 +1 If xq + xr + 1 is listed, then xq + xq-r +1 is also primitive.

45. Tausworthe Generator Homework: • Generate random numbers using the primitive polynomial x5+x2+1. (use l=4) • Generate the same sequence using LFSR.