Drawing from distributions Michel Bierlaire michel.bierlaire@epfl.ch Transport and Mobility Laboratory Drawing from distributions – p. 1/33
Discrete distributions • Let X be a discrete r.v. with pmf: P ( X = x i ) = p i , i = 0 , . . . , where � i p i = 1 . • The support can be finite or infinite. • The following algorithm generates draws from this distribution 1. Let r be a draw from U (0 , 1) . 2. Initialize k = 0 , p = 0 . 3. p = p + p k . 4. If r < p , set X = x k and stop. 5. Otherwise, set k = k + 1 and go to step 3. Inverse transform method Drawing from distributions – p. 2/33
Inverse Transform Method: illustration 0 0 . 24 0 . 66 0 . 77 1 p 1 = 0 . 24 p 2 = 0 . 42 p 3 = 0 . 11 p 4 = 0 . 23 Drawing from distributions – p. 3/33
Discrete distributions Acceptance-rejection technique • Attributed to von Neumann. • Mostly useful with continuous distributions. • We want to draw from X with pmf p i . • We know how to draw from Y with pmf q i . Define a constant c ≥ 1 such that p i ≤ c ∀ i s.t. p i > 0 . q i Algorithm: 1. Draw y from Y 2. Draw r from U (0 , 1) 3. If r < p y cq y , return x = y and stop. Otherwise, start again. Drawing from distributions – p. 4/33
Acceptance-rejection: analysis Probability to be accepted during a given iteration: P ( Y = y, accepted ) P ( accepted | Y = y ) = P ( Y = y ) = q y p y /cq y p y = c Probability to be accepted: � P ( accepted ) = y P ( accepted | Y = y ) P ( Y = y ) p y � = cq y q y y = 1 /c. Probability to draw x at iteration n c ) n − 1 p x (1 − 1 P ( X = x | n ) = c Drawing from distributions – p. 5/33
Acceptance-rejection: analysis Therefore, + ∞ � P ( X = x ) = P ( X = x | n ) n =1 � n − 1 p x + ∞ � 1 − 1 � = c c n =1 cp x = c = p x . Reminder: geometric series + ∞ 1 x n = � 1 − x n =0 Drawing from distributions – p. 6/33
Acceptance-rejection: analysis Remarks: • Average number of iterations: c • The closer c is to 1, the closer the pmf of Y is to the pmf of X . Drawing from distributions – p. 7/33
Continuous distributions Inverse Transform Method Idea: • Let X be a continuous r.v. with CDF F X ( ε ) • Draw r from a uniform U (0 , 1) • Generate F − 1 X ( r ) . Motivation: • F X is monotonically increasing • It implies that ε 1 ≤ ε 2 is equivalent to F X ( ε 1 ) ≤ F X ( ε 2 ) . Drawing from distributions – p. 8/33
Inverse Transform Method 1 F X ( ε ) 0 . 9 0 . 8 0 . 7 0 . 6 0 . 5 0 . 4 0 . 3 0 . 2 0 . 1 0 − 4 − 2 0 2 4 Drawing from distributions – p. 9/33
Inverse Transform Method More formally: • Denote F U ( ε ) = ε the CDF of the r.v. U (0 , 1) • Let G be the distribution of the r.v. F − 1 X ( U ) Pr( F − 1 G ( ε ) = X ( U ) ≤ ε ) Pr( F X ( F − 1 = X ( U )) ≤ F X ( ε )) = Pr( U ≤ F X ( ε )) = F U ( F X ( ε )) = F X ( ε ) Drawing from distributions – p. 10/33
Inverse Transform Method Examples: let r be a draw from U (0 , 1) Name F X ( ε ) Draw 1 − e − ε/b Exponential( b ) − b ln r µ − σ ln( 1 Logistic( µ , σ ) 1 / (1 + exp( − ( ε − µ ) /σ )) r − 1) ( ε/σ ) n σr 1 /n Power( n , σ ) Note: the CDF is not always available (e.g. normal distribution). Drawing from distributions – p. 11/33
Continuous distributions Rejection Method • We want to draw from X with pdf f X . • We know how to draw from Y with pdf f Y . Define a constant c ≥ 1 such that f X ( ε ) f Y ( ε ) ≤ c ∀ ε Algorithm: 1. Draw y from Y 2. Draw r from U (0 , 1) 3. If r < f X ( y ) cf Y ( y ) , return x = y and stop. Otherwise, start again. Drawing from distributions – p. 12/33
Rejection Method: example Draw from a normal distribution • Let ¯ X ∼ N (0 , 1) and X = | ¯ X | 2 π e − ε 2 / 2 , 0 < ε < + ∞ 2 • Probability density function: f X ( ε ) = √ • Consider an exponential r.v. with pdf f Y ( ε ) = e − ε , 0 < ε < + ∞ • Then f X ( ε ) 2 e ε − ε 2 / 2 f Y ( ε ) = √ 2 π • The ratio takes its maximum at ε = 1 , therefore f X ( ε ) f Y ( ε ) ≤ f X (1) � f Y (1) = 2 e/π ≈ 1 . 315 . 2 = e − ( ε − 1)2 √ e e ε − ε 2 / 2 = e ε − ε 2 f X ( ε ) 2 − 1 1 • Rejection method, with cf Y ( ε ) = 2 Drawing from distributions – p. 13/33
Rejection Method: example Draw from a normal 1. Draw r from U (0 , 1) 2. Let y = − ln(1 − r ) (draw from the exponential) 3. Draw s from U (0 , 1) 4. If s < e − ( y − 1)2 return x = y and go to step 5. 2 Otherwise, go to step 1. 5. Draw t from U (0 , 1) . 6. If t ≤ 0 . 5 , return x . Otherwise, return − x . Note: this procedure can be improved. See Ross, Chapter 5. Drawing from distributions – p. 14/33
Draws from the exponential Exponential 9000 8000 7000 6000 5000 4000 3000 2000 1000 0 Drawing from distributions – p. 15/33
Rejected draws Rejected draws 9000 8000 7000 6000 5000 4000 3000 2000 1000 0 Drawing from distributions – p. 16/33
Accepted draws Accepted draws 9000 8000 7000 6000 5000 4000 3000 2000 1000 0 Drawing from distributions – p. 17/33
Rejected and accepted draws Accepted draws 9000 Rejected draws 8000 7000 6000 5000 4000 3000 2000 1000 0 Drawing from distributions – p. 18/33
The polar method Draw from a normal distribution • Let X ∼ N (0 , 1) and Y ∼ N (0 , 1) independent • pdf: 1 1 e − y 2 / 2 = 1 e − x 2 / 2 2 π e − ( x 2 + y 2 ) / 2 . f ( x, y ) = √ √ 2 π 2 π • Let R and θ such that R 2 = X 2 + Y 2 , and tan θ = Y/X . ( X, Y ) R θ Drawing from distributions – p. 19/33
The polar method Change of variables (reminder): • Let A be a multivariate r.v. distributed with pdf f A ( a ) . • Consider the change of variables b = H ( a ) where H is bijective and differentiable • Then B = H ( A ) is distributed with pdf � dH − 1 ( b ) � �� f B ( b ) = f A ( H − 1 ( b )) � � � det � . � � db Here: A = ( X, Y ) , B = ( R 2 , θ ) = ( T, θ ) � � � � 1 2 T − 1 1 1 2 cos θ dH − 1 ( B ) 2 cos θ 2 sin θ T − T H − 1 ( B ) = = 1 2 T − 1 1 1 2 sin θ 2 sin θ 2 cos θ dB T T Drawing from distributions – p. 20/33
The polar method � � � � 1 2 T − 1 1 1 2 cos θ dH − 1 ( b ) 2 cos θ 2 sin θ T − T H − 1 ( b ) = = 1 2 T − 1 1 1 2 sin θ 2 sin θ 2 cos θ db T T Therefore, � � dH − 1 ( b ) �� � = 1 � � � det 2 . � � db and f B ( T, θ ) = 1 1 2 π e − T/ 2 , 0 < T < + ∞ , 0 < θ < 2 π. 2 Product of • an exponential with mean 2: 1 2 e − T/ 2 • a uniform on [0 , 2 π [ : 1 / 2 π Drawing from distributions – p. 21/33
The polar method Therefore, • R 2 and θ are independent • R 2 is exponential with mean 2 • θ is uniform on (0 , 2 π ) Algorithm: 1. Let r 1 and r 2 be draws from U (0 , 1) . 2. Let R 2 = − 2 ln r 1 (draw from exponential of mean 2) 3. Let θ = 2 πr 2 (draw from U (0 , 2 π ) ) 4. Let √− 2 ln r 1 cos(2 πr 2 ) X = R cos θ = √− 2 ln r 1 sin(2 πr 2 ) Y = R sin θ = Drawing from distributions – p. 22/33
The polar method Issue: time consuming to compute sine and cosine Solution: generate directly the sine and the cosine • Draw a random point ( s 1 , s 2 ) in the circle of radius one centered at (0 , 0) . • How? Draw a random point in the square [ − 1 , 1] × [ − 1 , 1] and reject points outside the circle • Let ( R, θ ) be the polar coordinates of this point. • R 2 ∼ U (0 , 1) and θ ∼ U (0 , 2 π ) are independent R 2 s 2 1 + s 2 = 2 cos θ = s 1 /R sin θ = s 2 /R Drawing from distributions – p. 23/33
The polar method Original transformation: √− 2 ln r 1 cos(2 πr 2 ) X = R cos θ = √− 2 ln r 1 sin(2 πr 2 ) Y = R sin θ = Replace r 1 by t = R 2 ∼ U (0 , 1) , and the sine and cosine as described above s 2 1 + s 2 t = 2 √ � − 2 ln t s 1 − 2 ln t X = R cos θ = = s 1 √ t t √ � − 2 ln t s 2 − 2 ln t Y = R sin θ = = s 2 √ t t Drawing from distributions – p. 24/33
The polar method Algorithm: 1. Let r 1 and r 2 be draws from U (0 , 1) . 2. Define s 1 = 2 r 1 − 1 and s 2 = 2 r 2 − 1 (draws from U ( − 1 , 1) ). 3. Define t = s 2 1 + s 2 2 . 4. If t > 1 , reject the draws and go to step 1. 5. Return � � − 2 ln t − 2 ln t x = s 1 and y = s 2 . t t Drawing from distributions – p. 25/33
Transformations of standard normal • If r is a draw from N (0 , 1) , then s = br + a is a draw from N ( a, b 2 ) • If r is a draw from N ( a, b 2 ) , then e r is a draw from a log normal LN ( a, b 2 ) with mean e a +( b 2 / 2) and variance e 2 a + b 2 ( e b 2 − 1) Drawing from distributions – p. 26/33
Multivariate normal • If r 1 ,. . . , r n are independent draws from N (0 , 1) , and r 1 . . r = . r n • then s = a + Lr is a vector of draws from the n -variate normal N ( a, LL T ) , where • L is lower triangular, and • LL T is the Cholesky factorization of the variance-covariance matrix Drawing from distributions – p. 27/33
Multivariate normal Example: ℓ 11 0 0 L = ℓ 21 ℓ 22 0 ℓ 31 ℓ 32 ℓ 33 s 1 = ℓ 11 r 1 s 2 = ℓ 21 r 1 + ℓ 22 r 2 s 3 = ℓ 31 r 1 + ℓ 32 r 2 + ℓ 33 r 3 Drawing from distributions – p. 28/33
Recommend
More recommend
