- 203 Views
- Uploaded on
- Presentation posted in: General

CS433 Modeling and Simulation Lecture 15 Random Number Generator

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - 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

- 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

- 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

- Why Random Number Generators?
- Desired Properties of RNGs
- Lehmer’s Algorithm
- Period and full-period RNGs
- Modulus and multiplier selection (Lehmer)
- Implementation - overflow

- 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

- 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: 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

- 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

- 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(.)

- 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

- 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}

- 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

- 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

- [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,

- 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.

- 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):

- 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

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)

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

}

- 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

- Problem: Compute without overflow
- General Idea: Perform modoperation first, before multiplication.
- Suppose that (not prime)
- We have:
- Thus, No overflow

- For the case, m is prime, so let
- q = quotient; r = remainder

- Let
- It can be shown that (Page 59, Lemis/Park Textbook)
and

Random number variants

Discrete random variables

Continuous random variables

- 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)}

Random variate generation:

1. Select u, uniformly distributed (0,1)

2. Compute F*(u); result is random variate with

distribution f()

Cumulative Distribution Function of X:

F(x) = P(X≤x)

f(x=2)

- 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 distribution with mean