This preview has intentionally blurred sections. Sign up to view the full version.
View Full DocumentThis preview has intentionally blurred sections. Sign up to view the full version.
View Full DocumentThis preview has intentionally blurred sections. Sign up to view the full version.
View Full DocumentThis preview has intentionally blurred sections. Sign up to view the full version.
View Full DocumentThis preview has intentionally blurred sections. Sign up to view the full version.
View Full DocumentThis preview has intentionally blurred sections. Sign up to view the full version.
View Full DocumentThis preview has intentionally blurred sections. Sign up to view the full version.
View Full Document
Unformatted text preview: Generating Samples of Uniformly Distributed Random Variables Performing Monte Carlo or discreteevent simulations requires generating events having cer—
tain random characteristics. For example, we may have to throw darts at a 3 X 3 square,
where the locations of these darts are uniformly distributed over this region. We may have
to generate the arrival times of the customers at a check—in counter, where the time intervals between successive arrivals have a certain probability distribution. We start by discussing how to generate samples from the uniform distribution over [0, 1].
It turns out that one can “transform” a uniformly distributed random variable into any other
random variable with an arbitrary probability distribution. (Remember one of the previous
lectures where we transformed a uniformly distributed random variable over [0,1] into the
prize random variable, which was uniformly distributed over [1000, 3000].) Therefore, being
able to generate samples from the uniform distribution over [0,1] provides the basis for all randomness in a simulation. For the remainder, we will refer to samples from the uniform distribution over [0,1] as
“random numbers.” 1. Physical Methods The earliest approaches for generating random numbers were manual methods, such as throwﬁ
ing dice, tossing a coin, drawing numbers from an urn or picking cards from a deck. Later
on, objects that “appear” to behave randomly were utilized for random number generation,
such as picking digits from the Ithaca phone book or from the expansion of the number 7r, or
reading the milliseconds from your computer clock. In the middle of the 20th century, people
built devices whose sole purpose was to generate random numbers. For example, one such
device was based on obtaining signals from (apparently randomly) pulsating vacuum tubes. If one utilizes physical methods, then the random numbers can either be generated as
needed (slow!) or stored in memory before the simulation study (expensive!) Also, besides
being slow and memory—intensive, an important drawback of using physical methods is that
bias may exist in the device. Hem example, one can imagine a roulette wheel not being
able to produce numbers that are truly uniformly distributed over {1, . . . ,36} due to some
inaccuracy during construction. Another drawback of these methods is that if the random
numbers are generated as needed (without storing the sequence), one cannot later replicate
the random input sequence. This precludes one from using a variance reduction technique
called “common random numbers” that is often useful in debugging (more about this later in the course). Advantages of using physical methods are, assuming that the uniform model is correct for
the device, we obtain a sequence of “true” random numbers. Of course, you should realize that there are those who believe that the laws of nature are fundamentally deterministic.
(Remember the famous Einstein quote “God does not play dice with the universe”) If one
takes this point of view, then no device is capable of producing true random numbers. 22 2. Mathematical Algorithms Virtually all computer simulations today use mathematical algorithms to generate the required
random numbers. To understand what we mean by a “mathematical algorithm,” consider the following approach for generating random numbers: 1. Start with a k—digit integer.
2. Square the integer. If necessary, supply 0’s to the left of the result to make it a Qk—digit integer.
3. Extract the middle 19 digits and return to Step 2. For example, if one starts with the integer 8234, then one obtains the following sequence
of “random” integers: 8234 x 8234 = 67(7987)56
7987 X 7987 = 63(7921)69
7921 x 7921 : 62(7422)41
7422 x 7422 = 55(0860)84 .... To obtain a random number, one then divides the elements of the sequence by 10000. This
method was proposed by John von Neumann, and is called the mid—square method. A funda—
mental objection to the mid—square method is that the random numbers it generates are not
random at all! If one knows the initial integer, then the whole sequence is deterministically
deﬁned. This criticism applies to all the mathematical algorithms we consider for random
number generation. Speciﬁcally, all the mathematical algorithms we consider are determinis
tic recursive methods that generate a sequence of numbers that appear to be sampled from
the uniform distribution over [0, I]. For this reason, some people refer to the random numbers generated by mathematical algorithms as pseudorandom numbers. The obvious advantages of using mathematical algorithms are they are potentially very
fast and can be replicated at will. However, one needs to carry out extensive tests in order
to make sure that the sequence of random numbers generated by a mathematical algorithm
appear to be uniformly distributed over [0,1] and do not exhibit any correlation with each
other. Otherwise, the results of a simulation study based on these random numbers will not
be valid. (It has been shown that the mid—square method produces random numbers with
undesirable statistical properties and should never be used in a real simulation study. Also, if
one ever obtains a k—digit number such that after squaring, the middle is digits are zero, then the sequence clearly degenerates and produces 0 indeﬁnitely.) 3. Linear Congruential Generators The most commonly used random number generators in practice are linear congruential gen—
erators (LCG’s). An LOG generates a sequence X0, X1, . . . of integers through the recursion XnH I (aXn + 0) mod m, where a, c, and m are ﬁxed integers, respectively known as the multiplier, increment and 23 modulus. Note that the ﬁrst integer X0 has to be supplied by the user and is usually referred
to as the seed. Reminder — k mod m is the integer obtained by dividing k by m and computing the remain
der. For example, 14 mod 5 z 4, 26 mod 21 : 5, —8 mod 3 : 1. It follows that each Xn is an integer in the set {0, l, . . . ,m — 1}. To obtain a random
number, one sets
U A : X“
M Ewample — Suppose that m = 8, a. = 5 and c = 3. Then, the sequence of numbers produced
by starting from the seed 3 is ESQl Sill/U30! it?! ‘3
(515+?) Madge (It, 25
I. (sxrmmmc b
5 (9 v. (Ex Lumbar! 8c 1
(fxl13)m33: O From the two examples abOVe, one can make the following observations: 3 24 p.
1_ 1. An LCG produces periodic output, with a period no larger than m. 2. If the period of an LCG is equal to m with a certain starting seed, then the period of the
LOG is always m (irrespective of the starting seed). 3. If an LCG does not have full period (that is, its period is less than m), then its period may
be dependent on the choice of the starting seed. In general, full period is desirable due to following reasons: 1. One should never use the whole period of an LCG, otherwise dependencies between the
random numbers will occur. (In fact, one should use only a fraction of the whole period.)
Having full period increases the potential number of random numbers one can reliably use in a simulation. 2. If the generator does not have full period, then there are gaps in the output sequence. The nature of such gaps for an LCG is not well understood, and may create distinctly non—random
or correlated behavior in the output. 80, at a minimum, one should choose a, c and m to have full period. Even if an LCG has full period, the “granularity” of the generated random numbers will
be 1/m. That is, two random numbers U1, U2 generated by an LOG will always satisfy
U1 7 U2] 2 1 This is usually not a big problem, since when one uses a large m (which is always the case in practice), 1 /m is small. The following important theorem characterizes the conditions under which an LCG has full period. Theorem 1 e An LCG has full period if and only if all of the following three conditions
hold: 1. m and c are relatively prime (that is, the only positive integer that divides both m and c
is l). 2. If g is a prime number that divides, m, then (1 divides a. — 1. 3. If 4 divides m, then 4 divides a — 1. Note that m = 8, a = 5 and c = 3 satisfy the three conditions above. Therefore, this LCG
has full period. Some commonly used LCG’s use a modulus of the form m = 25, because computers use
binary arithmetic and computing the remainder after division by an integer of the form 21’ just amounts to retaining the last b bits of the binary integer. Consider an LOG with m = 2b and b 2 2. For this LCG to have full period, the three
conditions in the theorem above must be satisfied. Since m = 2b and b 2 2, 4 divides m.
Therefore, in order to satisfy Condition 3, a must be of the form a=4k+l, 25 for some positive integer it. Note that when m = 2b, the only prime number that divides m
is 2. Since a — 1 = 4k, Condition 2 is automatically satisﬁed. Finally, when m 2 21’, m and
c are relatively prime if and only if c is odd. Thus, when m = 2b, the theorem above reduces to the following one:
Theorem 2 7 An LCG with m = 2" has full period if and only if c is odd and 4 divides
a — 1. LCG’S run faster if the increment c is set to 0. Such an LCG is called a multiplicative
generator In practice, most LCG’S are MG’s. 4. Multiplicative Generators A multiplicative generator is of the form X7,“ = (aXn) mod m. For a multiplicative generator, 0 is a “trap” state. That is, if Xn : 0, then X” : n+1 ~
Xn+2 : . . . = 0 indeﬁnitely. Therefore, for an MG, one says that the generator has full period if the sequence X0, X1,X2, . . . visits all m — 1 integers {1, . . . , m — 1} before repeating itself. Theorem 3 — An MG has full period if m is prime and if the smallest integer 3' such that
a5 — 1 is divisible by m is m — 1. However, as the next theorem states, if an MG uses a modulus of the form m : 21', then the
largest period that can be obtained is 711/4. Theorem 4 — The largest possible period for an MG with m = 2b is m/4. This is achieved
when the starting seed is an odd number and a is of the form a23+8k or az5+8k, for some positive integer is. It turns out (fortuitously!) that while m : 231 is clearly not prime, m : 231 — 1 is. On a
32—bit machine, the largest representable integer is 231, since one bit is reserved for the sign.
The most widely—used multiplicative generators today use this choice of modulus. Recall that
computing remainders of the mod 21’ operator is easy. One can take advantage of the fact that
231 — 1 is close to 231 to speed up the implementation of such generators (see Bratley, Fox and Schrage (1987), A Guide to Simulation, pages 212213). 5. Theoretical Analysis of Linear Congruential Generators Since LCG’s produce random numbers according to a mathematical algorithm, it is not sur—
prising that they possess certain theoretical deﬁciencies. The most disturbing of these deﬁciencies is the fact that the random numbers generated
by LCG’s fall in the planes. In order to see what we mean by this, let U0, U1, U2, . . . be the
sequence of random numbers generated by an LOG (that is, Un = Xn/m, where m is the 26 modulus). We would expect U0, U1, U2, . . . to be uniformly distributed over [0,1] and exhibit
no correlation with each other. Therefore, if we take the points (U0, U1, U2,    , Udulja (UdﬂUd+1}Ud+2i"'1U2d—1)} (U2dyU2d+1;U2d—i2,v;U3d—1);~ in the d—dimensional space, then these points should be uniformly distributed over the d—
dimensional hypercube [0, 1]d. (That is, for example, the points (Uo,U1), (Dabs), (U4, (15),... should lie uniformly over the square [0, 1] X [0, Unfortunately, this is not the case for the random numbers generated by an LCG. If
U0, U1, U2, . . . is the sequence of random numbers generated by an LCG, then the points (U0,U1,U2,:Ud—1), (UdaUd+1aUd+2awUQdil): (U2diU2d+laU2d+2_a:U3dil)a lie in a relatively small number of parallel (d — l)—dimensional hyperplanes. (That is, for
example, the points (U0: U1): (U2: U3), (U4, U5), . . .
lie on a relatively small number of parallel lines.) In simulations of geometric phenomena, this hyperplane problem could create serious dif—
ﬁculties. For discrete—event simulations, it appears that this problem, while a concern, does
not seriously inhibit the use of LCG’S. Speciﬁcally, LOG’s are widely used in such simulations with good results, provided the generator has passed an appropriate battery of statistical
tests. Example 7 The following two ﬁgures plot {(U2,,,U2,,+1) : "n, 2 0, l, . . for two LCG’s, each
with full period. The ﬁrst one has m = 210, a = 37, c = 1, and the second one has m=210, (1:57, 0:1. 0.5 U_2n+1 U_2n+1 27 Another curious deﬁciency of LCG’s is their interaction with the Box—Muller method for
normal variate generation. If N0, N1, N2, . . . are samples from the N (0, 1) generated by using the random numbers U0, U1, U2, . . . coming from an LOG, then the pairs (N0, N1), (N2, N3), (N4, N5), . .. must lie on a spiral in the two—dimensional space. This is clearly undesirable. Example 7 The following is a plot of {(Ngn, N2n+1) : n = 0,1, . . obtained by the BoxMuller
method, where the underlying random numbers are generated by the LCG with a = 9, m = 221, c = 1. N 2n+1 However, it can be forcefully argued that, because LCG’s are relatively tractable from a
mathematical viewpoint, it is possible to understand where their theoretical defects lie. For
more complex generators, their corresponding theoretical pathologies may be much harder to
detect. Also one can use methods such as combining or shufﬂing to effectively tackle these problems. 6. Combining Generators As indicated earlier, a popular choice as random number generator is an MG with modulus
m : 231 — 1. Such generators, with appropriately chosen multipliers, can have a period of as
long as m — 1 = 231 — 2. While this period may seem enormous, it really is not all that big.
Given that 210 = 1024, we are talking about a period in the order of 2 billion. This may not be large enough for some applications for the following reasons: 1. Some models may require many billions of random numbers to execute properly. For ex—
ample, consider a detailed simulation of trafﬁc in the downtown core of a large city. Assuming
that we model individual vehicle movements, we potentially may need to track tens of thou
sands of vehicles. Over the course of a day, it may be that each vehicle is itself subjected to 28 many thousands of random disturbances, each one potentially requiring a random number.
So, in aggregate, such a simulation may require hundreds of billions of random numbers to simulate a single day, let alone many days. 2. Over the full period of an LCG, each integer in the period is visited once and only once.
Clearly, over the full period, the output of an LCG is far too regular, as compared with a truly
random sequence (in which random clustering effects are present). 80, one should never use
the full period of an LCG in a simulation. In fact, similar regularity problems can arise with
LCG’s whenever too large a fraction of the total period is used (say, greater than 1%). For
an LCG with a period in the billions, we are therefore talking about a “usable” portion on
the order of millions. It is easy to write simple simulations that require such a large number of uniform random variates. As a consequence, we sometimes may need generators with much larger periods. One way
of doing this is to take two MG’s Xn+1 = (aan) mod m1
Yn+l = (azyn) mOd m2: and to set
Zn 2 (Xn + Yn) mod m3. If the parameters (11, :12, m1, mg and 7713 are chosen appropriately, the sequence Zg, Z1, 22, . . .
can potentially exhibit an enormous period (basically on the order of mlmg). For example,
set (1.1 = 40014, a2 = 40692, m1 = 2147483563, m2 = 2147483399 and m3 2 m1. For com—
puters with 32—bit word length, fast and portable implementations are available (see L’Ecuyer
(1988), Efﬁcient and portable combined random number generators, Comm. of the ACM, 31, 742—749). 7. Shuﬂiing Generators Appearing to have come from the uniform distribution is not the only property that one should
look for in a sequence of random numbers generated by an LCG. An LOG should also be able
to generate sequences of random numbers that are not correlated with each other and do not
possess excessively long or short “runs” (by a run, we mean a monotonically increasing or decreasing sequence of random numbers). Shuffling can be an effective way of fixing an LOG that do not exhibit the later two types
of behavior. In shufﬂing, one LCG is used to generate the random numbers and another LCG
is used to change their positions in the sequence. Here is one simple shuffling algorithm: 1. Fix a vector length k. k 2 128 was originally suggested, but It 2 2 seems to work as well
as anything else. 2. Initialize the k—dimensional vector V 2 (V1, V2, . . . , l/Eg), such that W2U11V2=U2auka=Uka 29 i_l’ where U1, U2, . . . , U}, are the ﬁrst 13 random numbers generated by the ﬁrst LOG. Set the
iteration counter 11. 2 1. 3. By using the second LCG, generate an index number 3' uniformly distributed over {1, 2, . . . , 4. Return as the random number. 5. Set to be equal to UM”, where UH” is the (k + n)—th random number generated by the
ﬁrst LCG. Increment 7?. by 1. G0 to Step 3. Example — Let k = 4. Assume that the ﬁrst LCG has m = 21“, a. = 37, c = 1, and the
second LCG has m = 211, a = 57, c = 3. Note that both of these LCG’s have full periods. Starting with the seeds 600 and 300, the ﬁrst 10 random numbers provided by these LCG’s
are as follows: 0.586 0.681 0.186 0.866 0.051 0.880 0.557 0.597 0.078 0.892 0.990
0.146 0.351 0.013 0.725 0.332 0.927 0.854 0.708 0.330 0.816 0.509 Step 2. Then, the vector V is initialized to V = (0.681, 0.186, 0.866, 0.051) (assuming that we discard the ﬁrst random number from each LCG, since its value is ﬁxed by
the seed assignment). Step 3. Generate an index number from the second LCG by [0.351 X 4] = 2.
Step 4. Return V2 as the current number. That is, return 0.186. Step 5. Replace V2 by the next random number from the ﬁrst LOG. V becomes V = (0.681, 0.880, 0.866, 0.051) Step 3. Generate an index number from the second LCG by [0.013 X 4] = 1.
Step 4. Return V1 as the current number. That is, return 0.681. Step 5. Replace V1 by the next random number from the ﬁrst LCG. V becomes 1/ x (0.557, 0.880, 0.866, 0051) Note that the returned random numbers are now 0.186, 0.681, 0.866. . .. We have shufﬂed the
sequence of the ﬁrst LCGI 8. Testing Random Number Generators Using an LCG with full period is not enough to ensure that the LOG in question provides
a “desirable” sequence of random numbers. One must, at minimum, address the following
concerns before using a particular LCG for a serious simulation study: 30 1. If the full period of the LCG is used, then the generated random numbers will of course look like they truly come from the uniform distribution, because each number in {0, l / m, 2 / m, . . . , (m—
1) /m} is exactly produced once. But we already discussed that the full period of an LCG
should never be exhausted. Therefore, the question is “do the random numbers generated by the ﬁrst m% of the full period look like they are coming from the uniform distribution?” 2. We should also check whether they are correlated. There are also other concerns, such as the length of the “runs,” but we will not discuss them
here. Testing uniformity — Here, we are interested in the following statistical test: HO : {U1,... ,Un} come from UlO, 1]
H1 : {U1,. . . , UH} do not come from U[0,1]. Remember that the failure to reject H0 does not mean that {Uh . . . , Un} do really come from
U [0, 1]. It only means that the non—conformance to uniformity has not been detected by the
test. This is the nature of all statistical tests. We can never accept the null hypothesis through a statistical test, we can only fail to reject it. Also, for each test one should specify a level of signiﬁcance, say or. 0: is the probability of
rejecting the null hypothesis by mistake. That is,
o: = P{reject H0  Hg is true}. 1.1fyousetatoolarge,__ HAL i; Wm “M” q [ll o'i WU“ ex...
View
Full Document
 '08
 TOPALOGLU

Click to edit the document details