Transcript Lecture15

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/
24 May 2009
Dr. Anis Koubâa
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, 2007
Chapter 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
Random Inter-Arrival Time: Exp(A)
Random Service Time: Exp(m)
D1
A2
A3
The Concept





Goal: create random variates
for random variables X of a
certain distribution defined with
its CDF 0  F  x   u  1
Approach: Inverse cumulative
distribution function
Cumulative Distribution Function
CDF F(x)=u
u2
u1
0  F x   u  1
Generate a random variable
0 u 1
Determine x such that
F  x   u  F 1 u   x
X1Exp
X2Exp
X2Poisson
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:
a large integer m, let the set m = {1, 2, … m-1}
 Draw in an integer x  m randomly
 Compute: u = x/m
 For
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+1 value 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 x 0   m
The function g(.) is defined as:
g  x   a  x  mod m
The mod function produces the remainder after dividing the first
argument by the second
 u  
 More precisely: u modv  u  v 
v   where  x  is the largest


integer n such that n≤ x

Observations
x i 1  a  x i  mod m



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 prime and initial 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
x i 1  a  x i  mod m

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:
x i 1  a  x i  mod m

Once a value is repeated, all the sequence is then
repeated itself
x 0 , x 1 ,..., x i ,..., x i  p where x i  x i  p
 p is the period: number of elements before the first
repeat
 clearly p ≤ m-1
 Sequence:

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 x 0   m and the sequence x 0 , x 1, x 2 ,... is produced by the
Lehmer generator x i 1  a  x i  mod m where m is prime,
then there is a positive integer p with p  m  1 such that:


x 0 , x 1,..., x p 1 are all different

and x i  x i  p , i  0,1, 2,...
In addition,  m 1 mod p  0
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
m into 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 where i is 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
1;
p =
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
x i  a i  x 0 mod m i  0,1, 2,...
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
a i mod m   m i  0,1, 2,..., m  1
is also a full period multiplier relative to m if 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 a  x  mod m
 Obvious computation is to compute  a  x  first, then do mod
operation



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 a  x  mod m without
overflow
 General Idea: Perform mod operation first,
before multiplication.
 Suppose that m  a  q (not prime)

a  x  mod m  a  x  mod a  q   a   x mod q 
 We have:  x mod q   q  1
 Thus, a   x mod q   a  q  1  a  q  m
No overflow

Implementation Issues: Mathematical Solutions

For the case, m is prime, so let m  a  q  r
q


= quotient; r = remainder
Let

x 
  x   a   x mod q   r   

q 

 x   x    a  x 
    q   m 

It can be shown that (Page 59, Lemis/Park Textbook)
a  x  mod m    x   m    x 
and
0 if   x    m
 x   
1 if -   x    m
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 0  F  x   u  1
Approach: Inverse cumulative
distribution function
Cumulative Distribution Function
CDF F(x)=u
u2
u1
0  F x   u  1
Generate a random variable
0 u 1
Determine x such that
F  x   u  F 1 u   x
X1Exp
X2Exp
X2Poisson
X1Poisson
Generating Discrete Random Variates
PDF
F(x)
1.0
0.8
u
u = F(x)
0.5 = F(2)
x
f(x)
1
0.2
2
0.3
3
0.1
4
0.2
5
0.2
Random variate generation:
1. Select u, uniformly distributed (0,1)
2. Compute F*(u); result is random variate with
distribution f()
F(x)
1.0
0.8
u
0.6
x = F*(u)
2 = F(0.5)
0.6
0.4
0.4
f(x=2)
0.2
0.2
1
2
3
4
5 x
x
Cumulative Distribution Function of X:
F(x) = P(X≤x)
1
2
3
4
5 x
x
Inverse Distribution Function (idf) of X:
F*(u) = min {x: u < F(x)}
Discrete Random Variate Generation
 Bernoulli(p):
 Return 1 with probability p,
 Return 0 with probability 1-p
u  random   ;

if u  1  p  return 0;
else return 1;


 Geometric(p): f  x   p k  1  p 
 Number of Bernoulli trials
until first ‘0’)
u  random   ;

return log 1.0  u  log  p  ;
Uniform (a,b): equally likely to select an integer in interval [a,b]
u  random   ;

return a  u  b  a  1  ;
Exponential Random Variates

Exponential distribution with mean m
u  random   ;

return -m  log 1  u  ;