Week05-Lecture (20100921) - Simulation Modeling and...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right 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 5: Generating Samples of Arbitrary Random Variables (09/21/10 - 09/23/10) 1 • Announcement and Agenda For HW 3 Q 2(d), the function =CHIINV(p, k) returns the (1-p)-quantile of the chi-square distribution with k degrees of freedom • • • • NOTE: =NORMINV(p, 0, 1) returns the p-quantile of the standard normal To get the 95% quantile, you need to use =CHIINV(0.05,49) Please do not worry about this if you get it wrong. I’ve asked the TAs to disregard this problem in their grading. For ORIE 5581 Students ONLY: There will be a midterm at 7:30pm on Thursday, September 30 at Rhodes 471. • • • • Sample midterm is available on Blackboard (see the “Exam” Folder in Course Documents) Review Session: Tuesday (9/28) at 7:30pm in Rhodes 253 • I’ll go over the sample midterm questions Open-book and open-notes (no laptop, iPad, iPhone, iTouch etc) Please bring a calculator! 2 Motivation • • • • We know how to generate samples of uniformly distributed random variables over [0,1]. How can we use this to generate samples from more complex distributions such as exponential, normal, etc. ? In theory, we can “transform” a uniform [0,1] random variable into any other random variable with an arbitrary probability distribution. 3 Methods for Transformation: • • • Inversion Acceptance-Rejection Convolution 3 Inversion Method: Continuous RVs • Suppose X is a continuous random variable with a cumulative distribution function F(⋅) • • Goal: Generate samples of X Since F(⋅) is continuous and increasing, F-1(⋅) well-defined 1 0.9 0.8 0.7 0.6 0.5 F ( x) = u 0.4 0.3 0.2 0.1 0 -5 -3 -1 1 3 5 7 9 x = F 1 (u ) Theorem: If U is uniformly distributed on [0,1], then F-1(U) has the same distribution as X 4 Proof of the Inversion Method (Continuous RV) Theorem: If U is uniformly distributed on [0,1], then -1(U) has the same distribution as X F Proof: Pr F ￿ −1 What is the intuition here? Why does this really work? 1 ￿ ￿ ￿ −1 ￿ ￿ (U ) ≤ x = Pr F F (U ) ≤ F (x) = Pr {U ≤ F (x)} = F (x) Suppose the density and the cumulative distribution functions of X are given by 0.3 0.25 0.2 0.6 0.8 0.15 0.4 0.1 0.05 0 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 -5 -4 -3 -2 -1 0.2 0 0 1 2 3 4 5 6 7 8 9 10 Key Observation: Since the p.d.f. is the derivative of the c.d.f., the c.d.f. will have the steepest rise when the p.d.f. takes large value! 5 Intuition for the Inversion Method • Suppose we “bombard” the y-axis with the random numbers U1, U2,.... and trace the values of F-1(U1), F-1(U2), ..... 1 1 0.8 0.8 U2 0.6 0.6 0.4 0.4 U1 0.2 0.2 0 0 -5 -5 -4 -4 -3 -3 -2 -1 -1 0 F-1(U1) F-1(U2) 1 2 3 4 4 5 5 6 6 7 7 8 8 9 9 10 10 Key Observation: The values {F-1(U1), F-1(U2), ..... } will “cluster” around the points where the c.d.f. is steep, corresponding to where the p.d.f. takes large values. 6 Example: Exponential Distribution • • Suppose X has an exponential distribution with parameter λ Density and cumulative distribution functions: −λx f (x) = λe For any u ∈ [0,1], how do we compute F-1(u)? ￿ ￿ −1 u = F F −1 (u) u = F F (u) ￿ ￿ ⇔ ⇔ ⇔ ⇔ ⇔ u=1−e u=1− −1 −λF −1 (u) e−λF (u) ∀x ≥ 0 F (x) = 1 − e −λx ∀x ≥ 0 • • 1−u=e −1 −λF −1 (u) = ln(1 − u) 1 −1 F −1 (u) = − ln(1 − u) λ −λF −1 (u) −λF −1 (u) e To generate a sample from an exponential distribution with parameter λ, generate U ~ Uniform[0,1] and return -{ ln(1-U) } / λ. You can also use -{ ln(U) } / λ 7 Example: Weibull Distribution • • Suppose X has an Weibull distribution with shape and scale parameters α and β, respectively, where (both α and β are positive real numbers) Density and cumulative distribution functions: f (x) = F (x) = ￿ ￿ αβ (β x) 0 1−e 0 α−1 −(β x)α e if x ≥ 0, otherwise. −(β x)α if x ≥ 0, otherwise. F −1 1 1 /α (u) = (− ln(1 − u)) , β ∀u ∈ [0, 1) 8 Example: Piecewise Linear CDF • Suppose the p.d.f. of a random variable X is given by: 1/3, 2/3, f (x) = 0, if 0 ≤ x ≤ 1, if 1 < x ≤ 2, otherwise. 2/3 1/3 0 0 1 2 3 1 x/3 (2x − 1)/3, F (x) = 1, F −1 if 0 ≤ x ≤ 1, if 1 < x ≤ 2, if x > 2 2/3 1/3 0 (u) = ￿ 0 1 2 3 3u (1 + 3u)/2, if 0 ≤ u ≤ 1/3, if 1/3 < u ≤ 1 9 Inversion Method: Discrete Random Variables • Suppose that X is a discrete random variable taking values in the set {x , 1 x2,..., xn} where x1 ≤ x2 ≤⋅⋅⋅≤ xn and Pr{ X = xi } = p(xi). 0 p(x ) 1 p(x ) + p(x ) 1 2 . . . F (x) = p(x1 ) + p(x2 ) + · · · + p(xi ) . . . 1 Cumulative Distribution Function of X if x < x1 , if x1 ≤ x < x2 , if x2 ≤ x < x3 , if xi ≤ x < xi+1 , if xn ≤ x p( x1 ) + p ( x2 ) + p( x3 ) + p ( x4 ) + p ( x5 ) p( x1 ) + p( x2 ) + p( x3 ) + p ( x4 ) p( x1 ) + p( x2 ) + p ( x3 ) p( x1 ) + p( x2 ) p( x1 ) x1 x2 x3 x4 x5 This is a “rightcontinuous” step function! 10 Inversion Method for Discrete RVs • • To generate a sample of X, we generate U ~ Uniform[0,1] Return xi such that xi is the smallest value in which F(xi) ≥ U, that is, p(x1 ) + p(x2 ) + · · · + p(xi−1 ) < U ≤ p(x1 ) + p(x2 ) + · · · + p(xi ) Illustration: p( x1 ) + p ( x2 ) + p( x3 ) + p ( x4 ) + p ( x5 ) p( x1 ) + p( x2 ) + p( x3 ) + p ( x4 ) U p( x1 ) + p( x2 ) + p ( x3 ) p( x1 ) + p( x2 ) p( x1 ) x1 x2 x3 x4 F-1(U) x5 11 Example: Geometric Distribution • • • Suppose X has a geometric distribution with parameter p • Pr{ X = k } = (1-p)k-1p for k = 1, 2, ... k ￿ i=1 Cumulative Distribution Function: for k = 1, 2, ... F (k ) = Pr {X ≤ k } = Pr{X = i} = 1 − (1 − p)k For any u ∈ [0,1], F-1(u) is equal to an integer k such that F(k-1) < u ≤ F(k) F (k − 1) < u ≤ F (k ) ⇔ ⇔ ⇔ ⇔ 1 − (1 − p) k−1 (1 − p)k ≤ 1 − u < (1 − p)k−1 k ln(1 − p) ≤ ln(1 − u) < (k − 1) ln(1 − p) ln(1 − u) k≥ > (k − 1) ln(1 − p) < u ≤ 1 − (1 − p) k F −1 ln(1 − u) (u) = ln(1 − p) ￿ ￿ 12 Drawback of the Inversion Method • • Although the inversion method can “in theory” be applied to any probability distribution, it may be computationally expensive or intractable because F-1(⋅) may be difficult to compute. Example: Suppose X has a Beta(3,4) distribution. • • • Density Function: f(x) = 60x3(1-x)2 for 0 ≤ x ≤ 1 Cumulative Distribution: F(x) = 15x4 - 24x5 + 10x6 for 0 ≤ x ≤ 1 Computing F-1(⋅) requires finding the root of a 6th-order polynomial. 13 • • Acceptance-Rejection Method Suppose a random variable X has a probability density function f(⋅) taking positive values only over the interval [a,b], that is, we have f(x) = 0 whenever x ∉ [a,b] • Assume that f(⋅) is bounded above by M, that is, maxx∈[a,b] f(x) ≤ M We can “enclose” the p.d.f. in the rectangle [a,b] × [0,M] M f ( x) a b x 14 • • • Acceptance-Rejection Method STEP 1: Generate independent U1 ~ Uniform[0,1] and U2 ~ Uniform[0,1] STEP 2: Set Z1 = a + (b-a) U1 and Z2 = M U2 • • • (Z1, Z2) is a point that is uniformly distributed over [a,b] × [0,M] If the point is below the p.d.f., we accept it and return the x-coordinate. NOTE: f(Z1) is large, then the probability of accepting Z1 is also high. M f ( x) STEP 3: If Z2 ≤ f(Z1), the return Z1. Otherwise, return to STEP 1. rejected accepted a b x 15 Analysis of Acceptance-Rejection Method • • • • Let Y denote the output of the Acceptance-Rejection Method • • Y corresponds to the x-coordinate Z1 when the point (Z1, Z2) is accepted For any x, Pr{ Y ≤ x} = Pr{ X ≤ x} Goal: We wish to show that Y has the same distribution as X Lemma 1: For any x, Pr{ Z1 ≤ x, Z2 ≤ f(Z1) } = F(x) / (M×(b-a)) Lemma 2: Pr { (Z1, Z2) is accepted } = 1 / (M×(b-a)) Pr {Y ≤ x} = Pr {Z1 ≤ x | (Z1 , Z2 ) is accepted } Proof of the correctness of the Pr {Z1 ≤ x and (Z1 , Z2 ) is accepted } = AcceptancePr {(Z1 , Z2 ) is accepted } Rejection Method: Pr {Z1 ≤ x, Z2 ≤ f (Z1 )} = Pr {(Z1 , Z2 ) is accepted } F (x) = (Lemma 1) M (b − a) Pr {(Z1 , Z2 ) is accepted } = F (x) (Lemma 2) 16 Proof of Lemma 1 M Pr{Z1 ≤ x, Z2 ≤ f(Z1)} corresponds to the probability that the point (Z1, Z2) falls in the shaded region! Pr {Z1 ≤ x, Z2 ≤ f (Z1 )} ￿x f (u)du F (x) a = = M (b − a) M (b − a) b a x Proof of Lemma 2 Pr {(Z1 , Z2 ) is accepted } = Pr {Z2 ≤ f (Z1 )} = Pr {Z1 ≤ b, Z2 ≤ f (Z1 )} F (b) 1 = = M (b − a) M (b − a) 17 Efficiency of Acceptance-Rejection Method • • • • Let N denote the number of pairs (Z1, Z2) that are needed to generate ONE sample under the acceptance-rejection method. • • • Let p = Pr { (Z1, Z2) is accepted } = 1 / (M×(b-a)) (by Lemma 2) N has a geometric distribution with parameter p Pr{ N = k } = (1-p)k-1p for k = 1, 2, .... E[N] = ∑k=1∞ k × Pr{ N = k} = ∑k=1∞ k (1-p)k-1p = 1 / p = M(b-a) So, the acceptance-rejection method is “efficient” when M(b-a) ≈ 1 Intuition: This means that M ≈ 1/(b-a) so that the density function is “close” to a uniform distribution 18 Example: Beta(2,2) Distribution • Suppose X has a Beta(2,2) distribution. Then, the p.d.f. of X is given by 1.5 f (x) = ￿ 6x(1 − x) 0 if 0 ≤ x ≤ 1, otherwise. 1.0 0.5 0 0 0.5 1.0 • • • • Apply the acceptance-rejection method with a = 0, b = 1, M = 1.5 STEP 1: Generate U1 ~ Uniform[0,1] and U2 ~ Uniform[0,1] • Z1 = U1 and Z2 = 1.5 U2 STEP 2: If 1.5 U2 ≤ 6 U1 (1 - U1), return U1. Otherwise, repeat STEP 1. The expected number of pairs (U1, U2) that are needed to generate a single sample is 1.5. 19 Convolution Method • • • Suppose the random variable X can be written as X = Y1 + Y2 + ⋅ ⋅ ⋅ + Yn • Y1,Y2,....,Yn are independent random variables To generate a sample of X under the convolution method, we generate samples of Y1,Y2,....,Yn independently and sum them up. Example: X ~ Binomial(n,p) • • • • X = Y1 + Y2 + ⋅ ⋅ ⋅ + Yn Y1,Y2,....,Yn are independent Bernoulli(p) random variables, that is, Pr{ Yi = 1} = p and Pr{ Yi = 0} = 1-p Generating a sample of each Yi is straightforward: If U ≤ p, return 1, otherwise, return 0. This method works well when n is not too large 20 Example: Gamma Distribution • Let X be a Gamma random variable with scale parameter λ and shape parameter n. The p.d.f. of X is given by f (x) = ￿ λn xn−1 e−λx (n−1)! 0 if x ≥ 0, otherwise. • If the shape parameter n is an integer, then X can be written as a sum of i.i.d. random variables, each having an exponential distribution with parameter λ • • • X = Y1 + Y2 + ⋅ ⋅ ⋅ + Yn Y1,Y2,....,Yn are independent random variables, each with an exponential distribution with parameter λ We can generate samples of Yi using the inversion method. 21 ...
View Full Document

This note was uploaded on 10/26/2010 for the course OR&IE 5580 at Cornell University (Engineering School).

Ask a homework question - tutors are online