Optimization and Simulation Drawing from distributions Michel - PowerPoint PPT Presentation

Optimization and Simulation Drawing from distributions Michel Bierlaire Transport and Mobility Laboratory School of Architecture, Civil and Environmental Engineering Ecole Polytechnique F´ ed´ erale de Lausanne M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 1 / 51

Discrete distributions Outline Discrete distributions 1 Continuous distributions 2 Transforming draws 3 Monte-Carlo integration 4 Summary 5 Appendix 6 M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 2 / 51

Discrete distributions 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 Inverse transform method 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. M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 3 / 51

Discrete distributions 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 M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 4 / 51

Discrete distributions Discrete distributions Acceptance-rejection 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. M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 5 / 51

Discrete distributions Acceptance-rejection: analysis Probability to be accepted during a given iteration P ( Y = y , accepted) = P ( Y = y ) P (accepted | 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 (1 − 1 c ) n − 1 p x P ( X = x | n ) = c M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 6 / 51

Discrete distributions Acceptance-rejection: analysis + ∞ � P ( X = x ) = P ( X = x | n ) n =1 � � n − 1 p x + ∞ � 1 − 1 = c c n =1 c p x = c = p x . Reminder: geometric series + ∞ � 1 x n = 1 − x n =0 M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 7 / 51

Discrete distributions 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 . M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 8 / 51

Continuous distributions Outline Discrete distributions 1 Continuous distributions 2 Transforming draws 3 Monte-Carlo integration 4 Summary 5 Appendix 6 M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 9 / 51

Continuous distributions Continuous distributions Inverse Transform Method 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 ). M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 10 / 51

Continuous distributions 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 M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 11 / 51

Continuous distributions 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 ( ε ) M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 12 / 51

Continuous distributions 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). M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 13 / 51

Continuous distributions 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. M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 14 / 51

Continuous distributions 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 . √ e e ε − ε 2 / 2 = e ε − ε 2 2 = e − ( ε − 1)2 2 − 1 f X ( ε ) 1 Rejection method, with cf Y ( ε ) = 2 M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 15 / 51

Continuous distributions Rejection Method: example Algorithm: 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. Otherwise, go to step 2 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, 2006, Chapter 5]. M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 16 / 51

Continuous distributions Draws from the exponential 9000 Exponential 8000 7000 6000 5000 4000 3000 2000 1000 0 M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 17 / 51

Continuous distributions Rejected draws 9000 Rejected draws 8000 7000 6000 5000 4000 3000 2000 1000 0 M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 18 / 51

Continuous distributions Accepted draws 9000 Accepted draws 8000 7000 6000 5000 4000 3000 2000 1000 0 M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 19 / 51

Continuous distributions Rejected and accepted draws 9000 Accepted draws Rejected draws 8000 7000 6000 5000 4000 3000 2000 1000 0 M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 20 / 51

Continuous distributions Drawing from the standard normal distribution Accept/reject algorithm is not e ffi cient Polar method: no rejection (see appendix) M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 21 / 51

Continuous distributions 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) M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 22 / 51

Continuous distributions 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 M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 23 / 51

Continuous distributions Multivariate normal Example: ⎛ ⎞ 0 0 ℓ 11 ⎝ ⎠ 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 M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 24 / 51

Transforming draws Outline Discrete distributions 1 Continuous distributions 2 Transforming draws 3 Monte-Carlo integration 4 Summary 5 Appendix 6 M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 25 / 51

Transforming draws Transforming draws Method Consider draws from the following distributions: normal: N (0 , 1) (draws denoted by ξ below) uniform: U (0 , 1) (draws denoted by r below) Draws R from other distributions are obtained from nonlinear transforms. Lognormal(a,b) � − (ln x − a ) 2 � 1 R = e a + b ξ f ( x ) = √ exp 2 b 2 xb 2 π M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 26 / 51

Transforming draws Transforming draws Cauchy(a,b) � � � 2 �� − 1 � x − a � � π ( r − 1 f ( x ) = 1 + R = a + b tan 2) π b b χ 2 ( a ) ( a integer) a � f ( x ) = x ( a − 2) / 2 e − x / 2 ξ 2 R = j 2 a / 2 Γ ( a / 2) j =1 Erlang(a,b) ( b integer) b � f ( x ) = ( x / a ) b − 1 e − x / a R = − a ln r i a ( b − 1)! j =1 M. Bierlaire (TRANSP-OR ENAC EPFL) Optimization and Simulation 27 / 51