Stochastic Simulation Generation of random variables Continuous - PowerPoint PPT Presentation

Stochastic Simulation Generation of random variables Continuous sample space Bo Friis Nielsen Institute of Mathematical Modelling Technical University of Denmark 2800 Kgs. Lyngby – Denmark Email: bfn@imm.dtu.dk

Plan W1.1-2 Plan W1.1-2 Random number generation Independent, uniformly distributed RN Distribution Independent Random variable Model Output DTU 02443 – lecture 4 2

Random variables Random variables Aim • The scope is the generation of independent random variables X 1 , X 2 , ... X n with a given distribution , F x ( x ) , (or probability density function [pdf]). • We assume we have access to a supply ( U i ) of random numbers, independent samples from the uniform distribution on ]0 , 1[ . • Our task is to transform U i into X i . DTU 02443 – lecture 4 3

Generation of (pseudo)random variates Generation of (pseudo)random variates • Inverse transformation techniques • Composition methods • Acceptance/rejection methods • Mathematical methods DTU 02443 – lecture 4 4

Uniform distribution I Uniform distribution I Our norm distribution or building brick, U i ∼ U (0 , 1) f ( x ) = 1 F ( x ) = x for 0 ≤ x ≤ 1 E ( U i ) = 1 1 Var ( U i ) = 2 12 1 0 1 DTU 02443 – lecture 4 5

Inverse transformation Inverse transformation The cumulative distribution function (CDF) F ( x ) = P ( X ≤ x ) F(x) F(x) = P(X<=x) x x x DTU 02443 – lecture 4 6

The random variable F ( X ) The random variable F ( X ) F(X) F(X) = U X X DTU 02443 – lecture 4 7

From U to X From U to X * U = F(X) * U = F(X) U * −1 X = F (U) DTU 02443 – lecture 4 8

Inversion method Inversion method The random variable U = F ( X ) U=F(X) X U = F ( X ) F ( x ) = P ( X ≤ x ) P ( U ≤ u ) = P ( F ( X ) ≤ u ) = P ( X ≤ F − 1 ( u )) = F ( F − 1 ( u )) = u I.e. F ( X ) is uniformly distributed. DTU 02443 – lecture 4 9

Inversion method Inversion method Assume g continuous and increasing and let: X = g ( Y ) Y ∼ F Y ( y ) = P ( Y ≤ y ) then F X ( x ) = P ( X ≤ x ) = P ( g ( Y ) ≤ x ) = P ( Y ≤ g − 1 ( x )) = F Y ( g − 1 ( x )) If Y = U then F U ( u ) = u , and F X ( x ) = g − 1 ( x ) . If X = F − 1 ( U ) then X will have the cdf F ( x ) ( ( F − 1 ( x )) − 1 = F ( x ) ). DTU 02443 – lecture 4 10

Uniform distribution II Uniform distribution II Now, focus on U ( a, b ) . 1 f ( x ) = a ≤ x ≤ b b − a F ( x ) = x − a F − 1 ( u ) = a + ( b − a ) u b − a X = a + ( b − a ) U ∼ U ( a, b ) DTU 02443 – lecture 4 11

Exponential distribution Exponential distribution The time between events in a Poisson process is exponentially distributed. (Arrival time) E ( X ) = 1 F − 1 ( u ) = − 1 F ( x ) = 1 − e ( − λx ) λ log(1 − u ) λ So (both 1 − U and U are uniformly distributed) X = − log( U ) ∼ exp( λ ) λ Exponential distribution 1 10 0.9 9 8 0.8 7 0.7 0.6 6 0.5 5 0.4 4 3 0.3 0.2 2 0.1 1 0 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 DTU 02443 – lecture 4 12

Pareto Pareto Is often used in connection to description of income (over a certain level). � β � k � � U − 1 X ∼ Pa ( k, β ) F ( x ) = 1 − x ≥ β X = β k x k k ( k − 1) 2 ( k − 2) β 2 E ( X ) = Var ( X ) = k > 1 , 2 k − 1 β Pareto distribution, Par(1,1) 1 Pareto with X ≥ 0 0.9 0.8 0.7 0.6 � − k � 1 + x 0.5 � � U − 1 k − 1 F ( x ) = 1 − X = β 0.4 0.3 β 0.2 0.1 0 1 2 3 4 5 6 7 8 9 10 DTU 02443 – lecture 4 13

Gaussian Gaussian X a result of many ( ∞ ) independent sources (Central limit theorem) X ∼ N ( µ, σ 2 ) Z = Φ − 1 ( U ) Z ∼ N (0 , 1) X = µ + σZ Gaussian distribution 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 −3 −2 −1 0 1 2 3 DTU 02443 – lecture 4 14

Mathematical Method Mathematical Method • By means of transformation and other techniques we can obtain a random variable with a certain distribution. The Box-Muller method A transformation from polar � ( θ = 2 πU 2 , r = − 2 log( U 1 ) ) into Cartesian coordinates ( X = Z 1 and Y = Z 2 ).      cos(2 πU 2 )  Z 1  = � − 2 log( U 1 ) Z 1 , Z 2 ∼ N (0 , 1)  sin(2 πU 2 ) Z 2 Central limit theorem n U i − n � X = eg. n = 6 2 i =1 DTU 02443 – lecture 4 15

Generation of cos and sin Generation of cos and sin Sine and cosine can be calculated by the following acceptance/rejection algorith m: 1. Generate V 1 , V 2 ∼ U ( − 1 , 1) 2. Generate R 2 = V 2 1 + V 2 2 3. If R 2 > 1 goto 1. 4. cos (2 πU 2 ) = V 1 R , sin (2 πU 2 ) = V 2 R DTU 02443 – lecture 4 16

LN ( α, β 2 ) LN ( α, β 2 ) Logarithmic Gaussian, LN ( α, β 2 ) Y ∼ LN ( α, β 2 ) log( Y ) ∼ N ( α, β 2 ) Y = e X X = α + β Z Z ∼ N (0 , 1) DTU 02443 – lecture 4 17

General and mulitvariate normal distribution General and mulitvariate normal distribution • Generate n independent values from an N (0 , 1) distribution, Z i ∼ N (0 , 1) . • X i = µ i + � i j =1 c ij Z j • Where c ij are the elements in the Cholesky factorisation of Σ , Σ = CC ′ DTU 02443 – lecture 4 18

Composition methods - hyperexponential Composition methods - hyperexponential distribution distribution m m p i e − λ i x = � � 1 − e − λ i x � � F ( x ) = 1 − p i i =1 i =1 Formally we can express Z = X I where I ∼ { 1 , 2 , . . . , m } with P ( I = i ) = p i and X I ∼ exp ( λ I ) 1. Choose I ∼ { 1 , 2 , . . . , m } with probabilities p i ’s 2. Z = − 1 λ I log ( U ) DTU 02443 – lecture 4 19

Composition methods - Erlang distribution Composition methods - Erlang distribution • The Erlang distribution is a special case of the Gamma distribution with integer valued shape parameter • An Erlang distributed random variable can be interpreted as a sum of independent exponential variables • We can generate an Erlang- n distributed random variate by adding n exponential random variates. E ( Y ) = n Var ( Y ) = n Y ∼ Erl n ( λ ) λ 2 λ with λ i = λ n n − 1 λ log ( U i ) = − 1 � � λ log (Π n Y = X i = i =1 U i ) i =1 i =1 DTU 02443 – lecture 4 20

Composition methods II Composition methods II Generalization: � f ( x ) = f ( x | y ) f ( y ) dy X given Y : f ( x | y ) Y : f ( y ) Y is typically a parameter (eg. the conditional distribution of X given Y = µ is N ( µ, σ 2 ) ) Generate: • Generate Y from f ( y ) . • Generate X from f ( x | y ) where Y is used. DTU 02443 – lecture 4 21

Acceptance/rejection Acceptance/rejection Problem: we would like to generate X from pdf f , but it is much faster to generate Y with pdf g . NB. X and Y have the same sample space. If f ( y ) g ( y ) ≤ c for all y and some c • Step 1. Generate Y having density g . • Step 2. Generate a random number U • If U ≤ f ( Y ) cg ( Y ) set X = Y .Otherwise return to step 1. g ( y ) dy f ( y ) cg ( y ) = f ( y ) dy c DTU 02443 – lecture 4 22

Statistics Toolbox Version 4.0 (R13) 20-Jun-2002 . . . Random Number Generators. betarnd - Beta random numbers. binornd - Binomial random numbers. chi2rnd - Chi square random numbers. exprnd - Exponential random numbers. frnd - F random numbers. gamrnd - Gamma random numbers. geornd - Geometric random numbers. hygernd - Hypergeometric random numbers. iwishrnd - Inverse Wishart random matrix. lognrnd - Lognormal random numbers. DTU 02443 – lecture 4 23

Exercise 3 Exercise 3 1. Generate simulated values from the following distributions (a) Exponential distribution (b) Normal distribution (at least with standard Box-Mueller) (c) Pareto distribution, with β = 1 and experiment with different values of k values: k = 2.05, k = 2.5, k = 3 and k = 4. Verify the results by comparing histograms with analytical results and perform tests for distribution type. 2. For the Pareto distribution with support on [ β, ∞ [ compare mean value and variance, with analytical results, which can be k calculated as E ( X ) = β k − 1 (for k > 1 ) and Var ( X ) = β 2 k ( k − 1) 2 ( k − 2) (for k > 2 ) 3. For the normal distribution generate 100 95% confidence intervals for the mean and variance, each based on 10 observations. Discuss the results.