Transcript Document

SIMULATION MODELING AND ANALYSIS
WITH ARENA
T. Altiok and B. Melamed
Chapter 4
Random Number and
Variate Generation
Altiok / Melamed Simulation Modeling and Analysis with Arena
Chapter 4
1
Random Number Generators
• Random Number Generators (RNG) are used
to introduce randomness into Monte Carlo
simulation
• An RNG is an algorithm that generates a stream
of pseudo-random numbers, depending on some
parameters and an initial seed
• Each stream “looks” random, but can be exactly
reproduced from the initial seed
• Basic RNGs typically produce uniformly
distributed sequences, and other distributions are
often obtained by their transformation
Altiok / Melamed Simulation Modeling and Analysis with Arena
Chapter 4
2
Linear Congruential Generators
• Generation algorithm has the form
x i = [ a x i - 1 + c ]( mod m ), i = 1, 2, ...
where
• the initial value x 0 , called the seed (or initial seed),
must be chosen with care (provided as part of the RNG)
• the n ( mod m ) operation returns the remainder of n / m
• the resulting pseudo-random integer sequence
{x n }, n = 0, 1, ¼ is in the range { 0, 1, ¼ , m - 1}
• to convert this sequence to a standard RNG sequence
in the range [0, 1], one uses the transformation u i = x i / m
Altiok / Melamed Simulation Modeling and Analysis with Arena
Chapter 4
3
Example: Linear Congruential RNG
Step i
0
1
2
3
4
5
6
7
8
9
10
xi
ui
3
6
5
0
7
2
1
4
3
6
5
3/8=0.375
0.750
0.625
0.000
0.875
0.250
0.125
0.500
0.375
0.750
0.625
Random numbers generated by a linear congruential RNG with
a = 5, c = 7, m = 8 and x 0 = 3
Altiok / Melamed Simulation Modeling and Analysis with Arena
Chapter 4
4
The Inverse Transform Method
(a) Every cdf, F (x ), is non-decreasing, so the inverse cdf, F -1 (u ) ,
X
X
is always well defined
(b) Applying the distribution function, F , to the underlying variate, X ,
X
results in a variate FX (X ) : Unif (0, 1) , uniform between 0 and 1
(c) Conversely, applying an inverse distribution function to a Unif(0,1)
variate, U , results in a variate X = F - 1 (u ) with cdf F (x )
X
X
(d) The Inverse Transform method uses (a) – (c) in a two-step algorithm:
1. Use your favorite RNG to generate a realization u from a variate,
U : Unif (0, 1)
-1
2. Compute x = FX (u ) as a realization of X
Altiok / Melamed Simulation Modeling and Analysis with Arena
Chapter 4
5
The Inverse Transform Method (Cont.)
FX (x )
1
u
0
x  FX1 (u )
Altiok / Melamed Simulation Modeling and Analysis with Arena
Chapter 4
6
Generation of Uniform Variates
• Suppose we wish to generate a realization x from the uniform distribution
Unif(2,10), for a realization u = 0.65 of the underlying RNG
• Recall that the uniform cdf can be rewritten as
x- a
u=
b- a
where u is given and x is unknown
• Solving the above for x to obtain the inverse cdf readily yields the formula
x = FX- 1 (u ) = (b - a ) u + a
and substituting a = 2, b = 10 and u = 0. 65 into the above
results in the requisite value x = (10 - 2) 0.65 + 2 = 7.2
Altiok / Melamed Simulation Modeling and Analysis with Arena
Chapter 4
7
Generation of Exponential Variates
• Suppose we wish to generate a realization x from the exponential distribution
Expo(0.5) of rate l = 0. 5 (mean 2), for a realization u = 0. 45 of the
underlying RNG.
• Recall that the cdf of the exponential distribution can be written as
u = 1- e - l x , where u is given and x is unknown
• Solving the above for x readily yields the formula
x = F - 1 (u ) = X
1
ln (1 - u )
l
and substituting l = 0. 5 and u = 0. 45 above results in the requisite value
x = - 2 ln (1 - 0.45) = 1.1957
• Since U ~ Unif (0, 1) implies 1 - U = W ~ Unif (0, 1) , it follows that the
equation above may be simplified into the equivalent formula
1
x = F - 1 (u ) = ln (u )
X
l
Altiok / Melamed Simulation Modeling and Analysis with Arena
Chapter 4
8
Generation of Discrete Variates
• Suppose we code the state space by integers, say, S = {1,2,…},
and we wish to generate a realization x from a discrete pdf over S
as in the table below
x
p X (x )
1
0.3
2
0.5
3
0.2
• Recall that the discrete cdf can be rewritten as
if x < 1
ìï 0,
ïï
u = í k
ïï å p , if k £ x < k + 1
ïî i = 1 i
Altiok / Melamed Simulation Modeling and Analysis with Arena
Chapter 4
9
Generation of Discrete Variates (Cont.)
• Utilizing the general formula for an inverse cdf, we deduce that
the inverse cdf can be written as
x = F - 1 (u ) = å k 1 (u )
F
kÎ S
X
where F = [
k
k- 1
k
å p , å p
n =1 n n =1 n
k
)
or equivalently, x = F - 1 (u ) = k such that
X
k- 1
k
n =1
n =1
å pn £ u < å pn )
• In our case, we have
ìï 1, for 0 £ u < 0.3
ïï
x = F - 1 (u ) = ïí 2, for 0. 3 £ u < 0. 8
X
ïï
ïï 3, for 0.8 £ u < 1
ïî
Altiok / Melamed Simulation Modeling and Analysis with Arena
Chapter 4
10
Generation of Step Variates
• Suppose we wish to generate a realization x from a step distribution with J = 4
steps, whose pdf is specified in the table below
i
l
r
i
p
i
i
1
0
2
0.1
2
2
4
0.3
3
4
6
0.4
4
6
8
0.2
• Recall that the step cdf can be rewritten as
ìï
0,
ïï
if
x < l
é
ù
1
p
ïï J
êj - 1
ú
j
u = í å 1
(x ) êå p + ( x - l )
ú, if l j £ x < r j
[
l
,
r
)
i
i
ïï j = 1 j j
r - l ú
êi = 1
j
j ú
ê
ë
û if
ïï
x > r
1
,
J
ïïî
Altiok / Melamed Simulation Modeling and Analysis with Arena
Chapter 4
11
Generation of Step Variates (Cont.)
• Solving the cdf for x yields the inverse step cdf
é
j - 1 ör - l
J
æ
ê
j
j
x = F - 1 (u ) = å 1 (u ) êl + ççu - å l ÷
X
j =1 Fj
ø p
i=1 i ÷
êj è
j
êë
j- 1
j
n =1
n =1
ù
ú
ú,
ú
ú
û
where F = [ å p n , å p n )
j
• It follows that the requisite variate realization x, corresponding to u = 0.5 is
x = F
- 1
X
2
æ
ö r3 - l3
( 0. 5) = l + çç0. 5 - å p ÷
=
÷ p
3
è
i=1 i ø
3
4 + ( 0.5 - [0.1 + 0.3])
6- 4
= 4.5
0.4
Altiok / Melamed Simulation Modeling and Analysis with Arena
Chapter 4
12
Generation of Step Variates (Cont.)
1.00
0.80
0.45
0.30
0
1
2
3
4
The Inverse Transform method for generating a variate with a step pdf
Altiok / Melamed Simulation Modeling and Analysis with Arena
Chapter 4
13
Generation of General Processes
• Generating general processes (non-iid) that follow a prescribed probability law
is hard. More specifically, the problem is two-fold:
1. To specify a probability law for a process, one must specify the temporal
dependence of the process in the form of joint distributions of any
dimension. Except for special cases (when dependence does not extend
arbitrarily far into the past), this approach can be impractical.
2. When temporal dependence cannot be estimated in terms of joint
distributions, one may elect to use a statistical proxy of temporal
dependence, often in the form of autocorrelations. Even when the
process to be modeled is stationary, the problem of devising a model that
fits the observed (empirical) autocorrelations and marginal distribution,
simultaneously, it is non-trivial.
• To illustrate the complexity inherent in the first problem, we outline an
algorithm for generating a general stochastic sequence, with a real valued state
space, S.
Altiok / Melamed Simulation Modeling and Analysis with Arena
Chapter 4
14
General Generation Algorithm
Input: a probability law L for a general random sequence {X n } ¥
n= 0
F (x ) = Pr { X £ x } , x Î S ,
X
0
0
1
0
0
F
X
0
(x 1 x 0 ) = Pr { X 1 £ x 1 X 0 £ x 0 } , x 0, x 1 Î S ,
F
X X
0
, given by
n +1
X
0
,K ,X
n
(x n + 1 x 0, ..., x n ) =
Pr{ X
n +1
£ x
n +1
X £ x , ..., X n £ x n } ,
0
0
x , ..., x
0
n +1
Î S.
Output: A sample path (realization) x 0, ..., x n , K from the probability law L.
Algorithm:
1. Apply the Inverse Transform method to the initial seed, u , to get x 0 = FX- 1 (u 0 )
0
0
2. Suppose x 0, ..., x n have already been generated from seeds u 0, ..., u n .
Then use the Inverse Transform method to generate the next variate
x n +1 = F - 1
(u n + 1 x 0, ..., x n ) , where u is the next seed
X
n +1
X ,K ,X n
n +1
0
3. Continue analogously as necessary…
Altiok / Melamed Simulation Modeling and Analysis with Arena
Chapter 4
15