280 likes | 589 Views
Al-Imam Mohammad Ibn Saud Islamic University College Computer and Information Sciences. CS433 Modeling and Simulation Lecture 15 Random Number Generator. http://www.aniskoubaa.net/ccis/spring09-cs433/. Dr. Anis Koubâa. 24 May 2009. Reading. Required
E N D
Al-Imam Mohammad Ibn Saud Islamic University College Computer and Information Sciences CS433Modeling and Simulation Lecture 15Random Number Generator http://www.aniskoubaa.net/ccis/spring09-cs433/ Dr. Anis Koubâa 24 May 2009
Reading • Required • Chapter 4:Simulation and Modeling with Arena • Chapter 2:Discrete Event Simulation - A First Course, • Optional • Harry Perros, Computer Simulation Technique - The Definitive Introduction, 2007Chapter 2 and Chapter 3
Goals of Today • Understand the fundamental concepts of Random Number Generators (RNGs) • Lean how to generate random variate for standard uniform distributions in the interval [0,1] • Learn how to generate random variate for any probability distribution
Outline • Why Random Number Generators? • Desired Properties of RNGs • Lehmer’s Algorithm • Period and full-period RNGs • Modulus and multiplier selection (Lehmer) • Implementation - overflow
Why Random Number Generators? • Random Numbers (RNs) are needed for doing simulations • Generate random arrivals (Poisson), random service times (Exponential) • Random numbers must be Independent (unpredictable) NOW A1 A3 A2 Random Inter-Arrival Time: Exp(A) Random Service Time: Exp(m) D1
The Concept • Goal: create random variates for random variables X of a certain distribution defined with its CDF • Approach: Inverse cumulative distribution function u2 CDF F(x)=u u1 • Cumulative Distribution Function • Generate a random variable • Determine x such that X2Exp X2Poisson X1Exp X1Poisson
Problem Statement of RNGs Problem: We want to create a function: u = rand(); that • produces a floating point number u, where 0<u<1 AND • any value strictly greater than 0 and less than 1 is equally likely to occur (Uniform Distribution in ]0,1[ ) • 0.0 and 1.0 are excluded from possible values
Problem Statement of RNGs • This problem can be simply resolved by the following technique: • For a large integer m, let the set m = {1, 2, … m-1} • Draw in an integer x mrandomly • Compute: u = x/m • m should be very large • Our problem reduces to determining how to randomly select an integer in m
Lehmer’s Algorithm • The objective of Lehmer’s Algorithm is to generate a sequence of random integers in m: x0, x1, x2, … xi, xi+1, … • Main idea: Generate the next xi+1value based on the last value of random integer xi xi+1 = g(xi) for some function g(.)
Lehmer’s Algorithm • In Lehmer’s Algorithm, g(.)is defined using two fixed parameters • Modulus m:a large, fixed prime integer • Multiplier a: a fixed integer in m • Then, choose an initial seed • The function g(.)is defined as: • The mod function produces the remainder after dividing the first argument by the second • More precisely: where is the largest integer n such that n≤ x
Observations • The mod function ensures a value lower than m is always produced • If the generator produces the value 0, then all subsequent numbers in the sequence will be zero (This is not desired) • Theorem if (m is primeandinitial seed is non-zero)then the generator will never produce the value 0 • In this case, the RNG produces values in m = {1, 2, … m-1}
Observations • The above equation simulates drawing balls from an urn without replacement, where each value in m represents a ball • The requirement of randomness is violated because successive draws are not independent • Practical Fact: The random values can be approximately considered as independent if the number of generated random variates (ball draws) is << m The Quality of the random number generator is dependent on good choices for a and m
The Period of a Sequence • Consider sequence produced by: • Once a value is repeated, all the sequence is then repeated itself • Sequence: where • p is the period: number of elements before the first repeat • clearlyp ≤ m-1 • It can be shown, that if we pick any initial seed x0, we are guaranteed this initial seed will reappear
Full Period Sequences • [LP] Discrete-Event Simulation: A First Course by L. M. Leemis and S. K Park, Prentice Hall, 2006, page. 42 Theorem 2.1.2 • If and the sequence is produced by the Lehmer generator where m is prime, then there is a positive integer p with such that: • are all different • and • In addition,
Full Period Sequences • Ideally, the generator cycles through all values in m to maximize the number of possible values that are generated, and guarantee any number can be produced • The sequence containing all possible numbers is called a full-period sequence(p = m-1) • Non-full period sequences effectively partition minto disjoint sets, each set has a particular period.
Modulus and Multiplier Selection Criteria • Selection Criteria 1:m to be “as large as possible” • m = 2i - 1 whereiis the machine precision (is the largest possible positive integer on a “two’s complement” machine) • Recall m must be prime • It happens that 231-1 is prime (for a 32 bit machine) • Unfortunately, 215-1 and 263-1 are not prime ;-( • Selection Criteria 2:p gives full-period sequence (p = m-1) • For a given prime number m, select multiplier a that provide a full period • Algorithm to test if a is a full-period multiplier (m must be prime):
Modulus and Multiplier Selection Criteria • Criterias • m to be “as large as possible” • p gives full-period sequence (p = m-1) Algorithm for finding if p is full-period multiplier p = 1; x = a; // assume, initial seed is x0=1, thus x1=a while (x != 1) { // cycle through numbers until repeat p++; x = (a * x) % m; // careful: overflow possible } if (p == m-1) // a is a full period multiplier else // a is not a full period multiplier
Other Useful Properties Theorem 2.1.1[LP]: If the sequence x0, x1, x2, … is produced by a Lehmer generator with multiplier a and modulus m, then Note this is not a good way to compute xi! Theorem 2.1.4[LP, p. 45]: If a is any full-period multiplier relative to the prime modulus m, then each of the integers is also a full period multiplier relative to mif and only if the integer i has no prime factors in common with the prime factors of m-1 (i.e., i and m-1 are relatively prime, or co-prime)
Other Useful Properties Generate all full-period multipliers // Given prime modulus m and any full period multiplier a, // generate all full period multipliers relative to m i = 1; x = a; // assume, initial seed is 1 while (x != 1) { // cycle through numbers until repeat if (gcd(i,m-1)==1) // x=aimod m is full period multiplier i++; x = (a * x) % m; // careful: overflow possible }
Implementation Issues • Assume we have a 32-bit machine, m=231-1 • Problem • Must compute • Obvious computation is to compute first, then do modoperation • The multiplication might overflow, especially if m-1 is large! • First Solution: Floating point solution • Could do arithmetic in double precision floating point if multiplier is small enough • Double has 53-bits precision in IEEE floating point • May have trouble porting to other machines • Integer arithmetic faster than floating point
Implementation Issues: Mathematical Solutions • Problem: Compute without overflow • General Idea: Perform modoperation first, before multiplication. • Suppose that (not prime) • We have: • Thus, No overflow
Implementation Issues: Mathematical Solutions • For the case, m is prime, so let • q = quotient; r = remainder • Let • It can be shown that (Page 59, Lemis/Park Textbook) and
Next Random number variants Discrete random variables Continuous random variables
The Concept • Goal: create random variates for random variables X of a certain distribution defined with its CDF • Approach: Inverse cumulative distribution function u2 CDF F(x)=u u1 • Cumulative Distribution Function • Generate a random variable • Determine x such that X2Exp X2Poisson X1Exp X1Poisson
F(x) F(x) 1.0 1.0 0.8 0.8 x = F*(u) 2 = F(0.5) u = F(x) 0.5 = F(2) 0.6 0.6 0.4 0.4 u 0.2 0.2 u x x 1 1 2 2 3 3 4 4 5 5 x x Inverse Distribution Function (idf) of X: F*(u) = min {x: u < F(x)} Generating Discrete Random Variates Random variate generation: 1. Select u, uniformly distributed (0,1) 2. Compute F*(u); result is random variate with distribution f() PDF Cumulative Distribution Function of X: F(x) = P(X≤x) f(x=2)
Discrete Random Variate Generation • Geometric(p): • Number of Bernoulli trials until first ‘0’) • Uniform (a,b): equally likely to select an integer in interval [a,b] • Bernoulli(p): • Return 1 with probability p, • Return 0 with probability 1-p
Exponential Random Variates • Exponential distribution with mean