This preview shows page 1. Sign up to view the full content.
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 (1p)quantile of the chisquare distribution with k degrees of freedom • • •
• NOTE: =NORMINV(p, 0, 1) returns the pquantile 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 Openbook and opennotes (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 AcceptanceRejection 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, F1(⋅) welldeﬁned 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 F1(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 yaxis with the random numbers U1, U2,.... and trace the values of F1(U1), F1(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 F1(U1) F1(U2)
1 2 3 4 4 5 5 6 6 7 7 8 8 9 9 10 10 Key Observation: The values {F1(U1), F1(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 F1(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(1U) } / λ. 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
F1(U) x5
11 Example: Geometric Distribution • • • Suppose X has a geometric distribution with parameter p • Pr{ X = k } = (1p)k1p 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], F1(u) is equal to an integer k such that F(k1) < 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 F1(⋅) may be difﬁcult to compute. Example: Suppose X has a Beta(3,4) distribution. • • • Density Function: f(x) = 60x3(1x)2 for 0 ≤ x ≤ 1 Cumulative Distribution: F(x) = 15x4  24x5 + 10x6 for 0 ≤ x ≤ 1 Computing F1(⋅) requires ﬁnding the root of a 6thorder polynomial. 13 • • AcceptanceRejection 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 • • • AcceptanceRejection Method
STEP 1: Generate independent U1 ~ Uniform[0,1] and U2 ~ Uniform[0,1] STEP 2: Set Z1 = a + (ba) 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 xcoordinate. 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 AcceptanceRejection Method • • • • Let Y denote the output of the AcceptanceRejection Method • • Y corresponds to the xcoordinate 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×(ba)) Lemma 2: Pr { (Z1, Z2) is accepted } = 1 / (M×(ba))
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 Efﬁciency of AcceptanceRejection Method • • • • Let N denote the number of pairs (Z1, Z2) that are needed to generate ONE sample under the acceptancerejection method. • • • Let p = Pr { (Z1, Z2) is accepted } = 1 / (M×(ba)) (by Lemma 2) N has a geometric distribution with parameter p Pr{ N = k } = (1p)k1p for k = 1, 2, .... E[N] = ∑k=1∞ k × Pr{ N = k} = ∑k=1∞ k (1p)k1p = 1 / p = M(ba) So, the acceptancerejection method is “efﬁcient” when M(ba) ≈ 1 Intuition: This means that M ≈ 1/(ba) 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 acceptancerejection 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} = 1p 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).
 '10
 PAAT

Click to edit the document details