Stochastic Simulation Generation of random variables Discrete - - PowerPoint PPT Presentation

Stochastic Simulation Generation of random variables Discrete sample space

Bo Friis Nielsen

Applied Mathematics and Computer Science Technical University of Denmark 2800 Kgs. Lyngby – Denmark Email: bfn@imm.dtu.dk

02443 – lecture 3 2

DTU

Plan W1.1-2 Plan W1.1-2

Random number generation Independent Random variable Output Independent, uniformly distributed RN Model Distribution

02443 – lecture 3 3

DTU

Random variables Random variables

Aim

• The scope is the generation of independent random variables

X1, X2, ... Xn with a given distribution, Fx(x), (or probability density function [pdf]).

• We assume we have access to a supply (Ui) of random numbers,

independent samples from the uniform distribution on ]0, 1[.

• Our task is to transform Ui into Xi.

02443 – lecture 3 4

DTU

Uniform distribution I Uniform distribution I

Our norm distribution or building block, U(0, 1) f(x) = 1 F(x) = x for 0 ≤ x ≤ 1

1 1

E(X) = 1

2

Var(X) =

1 12

02443 – lecture 3 5

DTU

Coin Coin

1 1 1/2 1 1/2

• r uniform distribution
• 1

1/2 1

X = 0, 1 P(X = i) = 1 2 X :=

• U > 1

2

• X = ⌊(2U)⌋

02443 – lecture 3 6

DTU

Bernoulli trial Bernoulli trial

Toss a coin with P(X = 1) = p and P(X = 0) = 1 − p. F(x) = P (X ≤ x) X = U > 1 − p X = 0 ≤ U ≤ 1 − p X = 1 1 − p < U ≤ 1

• 1

1 1−p

PDF CDF

1 p 1−p 1 1 1−p 1

02443 – lecture 3 7

DTU

A fair die A fair die

4 5 6 3 2 1 1/6 1/6 1 4 5 6 3 2 1

• r uniform distribution
• 1/6

1 1 3 5 2 4 6

X = 1, 2, ... 6 P(X = i) = 1/6 X = ⌊(6U)⌋ + 1 Can be generalized 6 → k.

02443 – lecture 3 8

DTU

Discrete distribution - direct (crude) method Discrete distribution - direct (crude) method

Suppose X can take k distinct values x1 < x2 < ... xk with pi = P(X = xi), i = 1, 2, ... , k Then X takes the value xi with probability pi if U falls in an interval with length pi. That is if

i−1

• j=1

pj < U ≤

i

• j=1

pj

• r

X = xi if F(xi−1) < U ≤ F(xi)

02443 – lecture 3 9

DTU

Geometric distribution, NB(1, p) Geometric distribution, NB(1, p)

The discrete time version of waiting time. Memory-less.

f(n) = P(X = n) = (1 − p)n−1p n = 1, 2, ... F(n) = P(X ≤ n) = 1 − (1 − p)n X = n if F(n−1) < U ≤ F(n) 1−(1−p)n−1 < U ≤ 1−(1−p)n n − 1 < log (1 − U) log (1 − p) ≤ n X =

• log (U)

log (1−p)

• + 1

02443 – lecture 3 10

DTU

Discrete distribution II Discrete distribution II

• 2

3 6 7 5 4 1 P1 P2 1

• 1. Generate U
• 2. Find the interval i which U belong to. Pi−1 < U ≤ Pi
• 3. output xi
• Linear search (E(X))
• Rearrangement of intervals
• Binary search
• Indexed search

02443 – lecture 3 11

DTU

Rejection Method Rejection Method

Simple rejection More optimistic: acceptance method.

Assume P(X = i) = pi for i = 1, 2, ... k.

• C
• 1

2 3 4 p(i)

Let c ≥ pi (then pi/c ≤ 1).

• 1. I = ⌊(k ∗ U1)⌋ + 1
• 2. if U2 ≤ pI/c output: I

Else goto 1. frequency for i :

1 k pi c

k

j=1 1 k pj c

= pi

02443 – lecture 3 12

02443

Alias method Alias method

• A method for generating discrete random variates of general

type

• From discrete uniform to general discrete
• Generate one random number
• One comparison
• Potentially one table look-up
• Drawback: Complex set-up procedure

02443 – lecture 3 13

DTU

A six-point distribution A six-point distribution

P(X = 1) = 17 96 P(X = 2) = 1 12 P(X = 3) = 1 3 P(X = 4) = 1 4 P(X = 5) = 1 24 P(X = 6) = 11 96

Alias method Alias method

• Setup procedure

⋄ Generate the table of F(I)-values, (which part of the mass belongs to I itself. ⋄ Generate the table of L(I)-values, (the alias of class I)

• Method at run time

⋄ Generate I:I = ⌊k ∗ U1⌋+1 ⋄ Test against F(I).If U2 ≤ F(I) then return X = I else return X = L(I). The L, F tables for the six-point distribution

F(1) = 1 F(2) = 1 2 F(3) = 15 16 F(4) = 1 F(5) = 1 4 F(6) = 11 16 L(1) = 1 L(2) = 4 L(3) = 1 L(4) = 4 L(5) = 3 L(6) = 3

02443 – lecture 3 15

DTU

Alias Method Alias Method

p(i)

• 1

2 3 4

On setup: generate F and L. Generation:

• 1. I = ⌊(k ∗ U1)⌋ + 1
• 2. if U2 ≤ F(I) output I

else output L(I).

The Alias tables The Alias tables

Generate F and L.

Pseudo code. p is a vector containing the probabilities.

• 1. L={1, . . . , k}
• 2. F=k*p (F = 1 is equivalent for the uniform dist.)
• 3. G=find(F>=1) and S=find(F<=1)
• 4. while ~isempty(S),

(a) i=G(1) and j=S(1) (b) L(j)=i and F(i)=F(i)-(1-F(j)) (c) if F(i)<1-eps then G(1)=[] and S=[S i] (d) S(1)=[]

02443 – lecture 3 17

DTU

Rejection Method Rejection Method

General method

Aim: We will generate X with probabilities pi = P (X = i). Assume Y with probabilities qi = P (Y = i) is easily generated and C ≥ pi

qi for all i = 1, . . . .

• 1. Generate Y with probability qi and let X⋆ = Y .
• 2. Generate U2.

If U2 ≤

pX⋆ CqX⋆ output X = X⋆

else reject.

02443 – lecture 3 18

DTU

Rejection Method: Probability for X = i: Rejection Method: Probability for X = i:

P(X = i) = P(X⋆ = i|accept) = P(X⋆ = i, accept) P(accept) = P(X⋆ = i)P(accept)|X⋆ = i) P(accept) = qi ·

pi Cqi

• j qj

pj Cqj

= pi

Excercise 2 Excercise 2

Discrete random variables

In the excercise you can use a build in procedure for generating random

• numbers. Compare the results obtained in simulations with expected results.

Use histograms (and tests).

• 1. Choose a value for the probability parameter p in the geometric di-

stribution and simulate 10,000 outcomes. You can experiment with a small, moderate and large value if you like.

• 2. Simulate the 6 point distribution with

X 1 2 3 4 5 6 pi 7/48 5/48 1/8 1/16 1/4 5/16 (a) by applying a direct (crude) method (b) by using the the rejction method (c) by using the Alias method