Week04-Lecture (20100914)

Week04-Lecture (20100914) - Simulation Modeling and...

Info iconThis preview shows pages 1–5. Sign up to view the full content.

View Full Document Right Arrow Icon

Info iconThis preview has intentionally blurred sections. Sign up to view the full version.

View Full DocumentRight Arrow Icon

Info iconThis preview has intentionally blurred sections. Sign up to view the full version.

View Full DocumentRight Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: Simulation Modeling and Analysis (ORIE 4580/5580/5581) Week 4: Generating Samples of Uniformly Distributed Random Variables (09/14/10 - 09/16/10) 1 Announcement and Agenda • • • • HW 3 is available on Blackboard and will be due by 11am on 9/23. For ORIE 5581 Students ONLY: There will be a mid-midterm at 7:30pm on Thursday, September 30 at Rhodes 471. Sample mid-midterm will be given on 9/23. Quiz Generating samples of uniformly distributed random variables 2 Motivation • • Performing Monte Carlo or Discrete Event simulation requires generating events having certain random characteristics. • • • • • Capacity Allocation Game Example: We might need to generate samples of demands having binomial distribution. Simplest and Most Fundamental Question: How do we simulate a number that is uniformly distributed over the interval [0,1]? Basis for all randomness in simulations In theory, we can transform a Uniform[0,1] RV into a random variable with any arbitrary distribution. Why should you care? Just use RAND() command in Excel? • RAND() in earlier versions of Excel sometimes produces unacceptably short cycles. Need to be careful! Terminology: For the remainder of this lecture, we will refer to samples from the uniform distribution over [0,1] as “random numbers”. 3 Physical Methods • Earliest approach for generating random numbers were manual • • Earlier Example: throwing dice, tossing a coin, drawing number from an urn, or picking cards from a deck. Later Examples: Use objects that “appear” to behave randomly such as digits in the phone book, digits in the expansion of π, or reading the milliseconds from your computer clock. • • Advantages: If the physical device is indeed random, then we obtain a sequence of “true” random numbers! Drawbacks of Physical Methods: • • • • Slow Expensive (need to store these numbers before simulation study) Bias: construction defect associated with the physical device Cannot replicate the random input sequence 4 Mathematical Algorithms • • Virtually all computer simulations today use mathematical algorithms to generate the required random numbers. Mid-Square Method: An example of a mathematical algorithm for generating random numbers proposed by von Neumann: • • • STEP 1: Start with a k-digit integer STEP 2: Square the integer. If necessary, pad 0’s to the left of the number so that the result is a (2k)-digit integer STEP 3: Extract the middle k digits and return to STEP 2 • • • • • Example: Suppose we start with 8234. Then, we have the following sequence of “random” integers: 8234 × 8234 = 67(7987)56 7987 × 7987 = 63(7921)69 7921 × 7921 = 62(7422)41 7422 × 7422 = 55(0860)84 To obtain a random number, we divide the elements in the sequence by 10000 5 Criticism: Pros and Cons of Math Alg for RNG • Pros: Speed and replication • Cons: The random numbers generated are NOT random at all! • If we know the initial integer, then the whole sequence is deterministically defined! • This criticism applies to ALL the mathematical algorithms that we consider for random number generation. • • • All mathematical algorithms that we will consider are deterministic recursive methods that generate a sequence of numbers that appear to be uniformly distributed over [0,1]. For this reason, some people refer to the random numbers generated by mathematical algorithms as pseudo-random numbers Need to conduct extensive tests in order to make sure that the sequence of numbers generated appear to be uniformly distributed and do not exhibit any correlation with each other. Otherwise, the simulation study based on these numbers will not be valid. 6 Linear Congruential Generators (LCG) • • Most commonly used random number generators (RNG) Inputs to LCG: • • • • • • • Multiplier: a Increment: c Modulus: m Seed: X0 LCG generates a sequence of integers X1, X2,... using the recursion: Xn+1 = (aXn + c) mod m NOTE: k mod m is an integer in the set {0, 1, 2, ..., m-1} obtained by dividing k by m and computing the remainder • Example: 14 mod 5 = 4, 26 mod 21 = 5, -8 mod 3 = 1 7 To obtain a random number, we set Un = Xn / m LCG Examples m = 8, a = 5, c = 3 n 0 1 2 3 4 5 6 7 8 9 Xn 3 2 5 4 7 6 1 0 3 2 Un 3/8 2/8 5/8 4/8 7/8 6/8 1/8 0/8 3/8 2/8 Xn+1 2 = (5×3 + 3) mod 8 5 = (5×2 + 3) mod 8 4 = (5×5 + 3) mod 8 7 = (5×4 + 3) mod 8 6 = (5×7 + 3) mod 8 1 = (5×6 + 3) mod 8 0 = (5×1 + 3) mod 8 3 = (5×0 + 3) mod 8 2 = (5×3 + 3) mod 8 5 = (5×2 + 3) mod 8 n 0 1 2 3 4 5 6 7 8 9 10 11 12 m = 32, a = 11, c = 0 (under different seeds) Xn 1 11 25 19 17 27 9 3 1 11 25 19 17 n 0 1 2 3 4 5 6 7 8 9 10 11 12 Xn 2 22 18 6 2 22 18 6 2 22 18 6 2 n 0 1 2 3 4 5 6 7 8 9 10 11 12 Xn 4 12 4 12 4 12 4 12 4 12 4 12 4 8 • • • • An LCG produces periodic output, with a period no larger than m. If the period of an LCG is equal to m (a “full period) with a particular starting seed, then it has a full period for ALL starting seeds. If an LCG does NOT have a full period (that is, its period is less than m), then its period may depend on the choice of the seed value. We prefer full period because we only use a fraction of the period generated (using the entire period causes dependencies between outputs). Observations • • • • Having a full period increases the potential choices of random numbers that we can use in the simulation If LCG does NOT have full periods, there are gaps in the output sequence. These gaps are NOT well understood, and may create non-random or correlated behavior in the output. GOAL: Choose a, c, and m so that the LCG has a full period! NOTE: The gaps between two random numbers under LCG is always at least 1/m (Why?) • In practice, this is ok because m is chosen to be very large! 9 LCG Theory and Application • Theorem: 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, their only common divisor is 1 2. If q is a prime number that divides m, then q divides a-1 3. If 4 divides m, then 4 divides a-1 • • • NOTE: m=8, a=5, and c=3 satisfy this condition Application: Commonly used LCG has a modulus m = 2b for some integer b because computers use binary arithmetic. • • • • The remainder after dividing an integer by 2b corresponds to the last b bits of the integer a = 4k + 1 for some positive integer k (Condition 3) Condition 2 is automatically satisfied because q = 2 Condition 1 holds when c is odd 10 Suppose m = 2b and b ≥ 2. To have a full period, we must have Multiplicative Generator (MG) • • • • • • A multiplicative generator (MG) is a special case of LCG with c = 0 • • It’s more computationally efficient Most LCGs in practice are MGs MG generates a sequence of integers X1, X2,... using the recursion: Xn+1 = (aXn ) mod m For MG, 0 is a “trap” state because if Xn = 0, then all subsequent terms will also be zero! So, we say that an MG has a “full” period if the sequence X0, X1 ,X2 .... visit all m-1 integers {1, 2, ..., m-1} before repeating itself. Unfortunately, for MG, when m = 2b, we do not have a full period! Theorem: 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 form a = 3+8k or a = 5+8k for some positive integer k. 11 When does an MG have a full period? • • Theorem: An MG has a full period if m is prime, and the smallest integer j such that aj - 1is divisible by m is m-1. Application: Fortunately, m = 231 - 1 = 2,147,483,647 is a prime number! • • The most commonly used MG uses this particular modulus, with a = 630,360,016 We can use the fact that 231 - 1 is close to 231 to speed up the computation. 12 • • • • Deficiencies of LCG Since LCGs produce random numbers based on a mathematical algorithm, it is not surprising that it contains theoretical deficiencies Key Problem: The numbers generated by LCGs “fall in the planes” Let U0, U1, U2, .... denote the random numbers under LCG (Ui = Xi / m) U0, U1, U2, .... should be independent and uniformly distributed on [0,1] • Consider the order pairs (U0, U1), (U2, U3), (U4, U5), .... . These points should be uniformly distributed over the square [0,1] × [0,1] 1 0.9 0.8 0.7 0.6 U_2n+1 Scatter Plot of { (U2n, U2n+1) : n = 0, 1,.... } for an LCG with a full period with m = 210, a = 37, and c = 1. These points fall on a small number of parallel lines in the square! 0.5 0.4 0.3 0.2 0.1 0 0 0.2 0.4 0.6 0.8 1 13 Hyperplane Problem • • • Suppose we consider the triplets { (U2n, U2n+1, U2n+2) : n = 0, 1,.... }. These triplets will fall on a relatively small number of 2-dimensional plane in the unit cube [0,1] × [0,1] × [0,1] In general, the collection { (U2n, U2n+1, ....,U2n+d-1) : n = 0, 1,.... } will fall on a relatively small number of parallel (d-1)-dimensional hyperplanes. Implication: For simulation of geometric phenomenon, this hyperplane problem can cause serious difficulties. • • For discrete-event simulations, however, this is not too serious. LCGs are widely used in discrete-event simulations with good results Caveats: The LCG parameters m, a, and c must be chosen appropriately and pass a battery of statistical tests! 14 • • • Spiral Problem LCG interacts strangely with the Box-Muller method for generating normal random variables (to be discussed later) If N0, N1, N2, .... are samples from N(0,1) generated by using the random numbers U0, U1, U2, .... coming from an LCG, then the pairs (N0, N1), (N2, N3), (N4, N5), .... must lie on a spiral in the 2-dimensional space Example: Scatter plot of { (N2n, N2n+1) : n = 0, 1,.... } obtained by BoxMuller method, where the underlying random numbers are generated by an LCG with a = 9, m = 221, and c = 1. 4 3 2 1 0 -4 -2 -1 0 2 4 -2 -3 Given these problems with LCGs, what should we do? Should we use it? How can improve it? 15 N_2n+1 -4 How to increase LCG Periods? • • With modulus m = 231 - 1 and appropriately chosen multiplier, we can get a period of m - 1 = 231 - 2 = 2,147,483,646 (roughly two billions) This is actually not that much because of the following reasons: • Consider a simulation model of traffic in Manhattan where we track the individual movement of each vehicle. • • • • Over the course of a day, each vehicle is subject to thousands of random disturbances, each requiring a random number. In aggregate, over the course of a single day, we may require hundreds of billions of random numbers! We NEVER use the full period of an LCG. Using more than 1% of the full period makes the numbers “too regular” So, the usable portion of the numbers is only in millions. 16 Combining Generators Xn+1 = (a1 Xn ) mod m1 Yn+1 = (a2 Yn ) mod m2 Zn+1 = (Xn + Yn ) mod m3 • • If the parameters a1, a2, m1, m2, and m3 are chosen appropriately, the sequence Z0, Z1, Z2, .... can have periods of order m1 × m2 Example: a1 = 40014, a2 = 40692, m1 = 2,147,483,563, m2 = 2,147,483,399, and m3 = m1. In this case, m1 × m2 ≈ 1020 17 • Shuffling Generators: Reducing Correlations Shuffling = Changing the positions of the numbers generated by an LCG • • • Reduce correlation among the numbers Reduce the number of excessively long “runs” (a run is a monotonically increasing or decreasing sequence of numbers) Let U1, U2, U3,.... denote the sequence of random numbers in [0,1] generated by the FIRST LCG. • • • • • • A Simple Shuffling Algorithm: Require 2 LCGs STEP 1: Fix a vector of length k STEP 2: Initialize the k-dimensional vector V = (V1,V2,...,Vk) so that V1 = U1,V2 = U2,....,Vk = Uk . Set iteration counter n = 1. STEP 3: Using the SECOND LCG, generate an index j that is uniformly distributed over {1, 2, ..., k} STEP 4: Output Xn = Vj as the random number STEP 5: Update the vector V by replacing Vj by Uk+n, which is the (k+n)th random number under the FIRST LCG. Increment n by 1. Go to STEP 3. 18 • • • • • • • Shuffling Example with k = 4 FIRST LCG: m = 210, a = 37, and c = 1. Seed = 600 SECOND LCG: m = 211, a = 57, c = 3. Seed = 300 The sequence of random numbers generated by each LCG 1st 2nd 0.586 0.146 0.681 0.351 0.186 0.013 0.866 0.725 0.051 0.332 0.880 0.927 0.557 0.854 0.597 0.708 0.078 0.330 0.892 0.816 Step 2: V is initialized to be (0.681, 0.186, 0.866, 0.051) and set n = 1. • We are ignoring the first random number from each LCG (shaded in gray) because they are fixed by the seed assignment. Step 3: Using 2nd LCG to generate an index j = 0.351 × 4 = 2 in {1, 2, 3, 4} Step 4: Output X1 = V2 = 0.186 Step 5: Update V by replacing V2 with the next random number from the first LCG. So, V becomes (0.681, 0.880, 0.866, 0.051). Increment n to 2. • Step 3: Using 2nd LCG to generate j = 0.013 × 4 2 1 =1 • Step 4: Output X = V = 0.681 • Step 5: Update V to (.557, 0.880, 0.866, 0.051). n = 3. We have shuffled the 19 sequence of the first LCG! Testing for Uniformity of Random Number Generators Statistical Test H0: {U1,..., Un} come from a uniform distribution on [0,1] H1: {U1,..., Un} do NOT from a uniform distribution on [0,1] • Important: Failure to reject H0 does NOT mean that {U1, ..., Un} actually come from a uniform distribution on [0,1]. • • • It only means that the non-conformance to uniformity has NOT been detected by the test We can never accept a null hypothesis through a statistical test. We can only fail to reject it! Key Idea: Divide the interval [0,1] into k equal subintervals. If H0 is true (that is, {U1,..., Un} do really come from a uniform distribution), then we would expect n/k numbers falling in each interval! Let fi = # of random numbers falling in the subinterval [(i-1)/k, i/k] k k ￿ (fi − (n/k ))2 ￿￿ k n ￿2 = fi − Test Statistics: T = n/k n i=1 k i=1 20 Testing for Uniformity (cont.) Theorem: Under H0, T has chi-square distribution with k-1 degrees of freedom Density and Cumulative Distribution Functions for Chi-Square RVs • Need to specific the significance level α • • This is probability of rejecting the null hypothesis by mistake α = Pr { Reject H0 | H0 is true } • • If we set α too large, then our test will make a lot of errors If α is too small, then the test will be “weak”, that is, it will almost always fail to reject the null hypothesis. Test: Reject H0 if T ≥ χ2 −1,1−α k (1-α)-quantile of the chi-square distribution with k-1 degrees of freedom. 21 • • • Testing for Correlation If X1, X2,... is a sequence of identically distributed RVs, then the covariance at lag j (at n) is defined as Cov(Xn, Xn+j) =E[ Xn Xn+j ] - E[Xn] E[Xn+j] The correlation at lag j is defined by: ρj = Cov(Xn, Xn+j) / Var(Xn) Application: Suppose the random variables {U1, U2,..., Un} are uniformly distributed on [0,1], then E[Un] = 1/2 and Var(Un) = 1/12 ￿1 1￿ E[Un Un+j ] − 2 · 2 ρj = = 12E[Un Un+j ] − 3 1 12 • Suppose {U1, U2,..., Un} are generated from an LCG, we want to test whether these numbers exhibit correlations among each other. Statistical Test: H0: ρj = 0 and H1: ρj ≠ 0 1 Tj = 12 U1+nj U1+(n+1)j − 3 h + 1 n=0 22 Test Statistics: Estimate the correlation at lag j by considering the sequence U1, U1+j, U1+2j,...., U1+hj, U1+(h+1)j, where h is the largest integer such that 1+(h +1)j ≤ n. ￿ h ￿ ￿ Testing for Correlation (cont.) Theorem: Under H0, the test statistic Tj has mean zero and variance (13h + 7)/(h+1)2. Moreover, for large n, Tj is approximately normally distributed. Therefore, ￿h 12 n=0 U1+nj U1+(n+1)j − 3(h + 1) D Tj − E[Tj ] ￿ √ = ≈ N (0, 1) 13h + 7 Var(Tj ) Test: Reject H0 if Tj − E[Tj ] ￿ > z1−α/2 Var(Tj ) or (1-α/2)-quantile of the standard normal Tj − E[Tj ] ￿ < zα/2 Var(Tj ) Note: You should repeat the same test for the sequences: U2, U2+j, U2+2j, U2+3j.... and for every lag j. 23 U3, U3+j, U3+2j, U3+3j.... U4, U4+j, U4+2j, U4+3j.... ...
View Full Document

This note was uploaded on 10/26/2010 for the course OR&IE 5580 at Cornell.

Page1 / 23

Week04-Lecture (20100914) - Simulation Modeling and...

This preview shows document pages 1 - 5. Sign up to view the full document.

View Full Document Right Arrow Icon
Ask a homework question - tutors are online