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
 '10
 PAAT
 Normal Distribution, Probability theory, probability density function, Cumulative distribution function, Geometric distribution

Click to edit the document details