stat341f11

Sampling - Sept 20, 2011

From [math]\displaystyle{ x \sim~ f(x) }[/math] sample [math]\displaystyle{ \,x_{1}, x_{2}, ..., x_{1000} }[/math]

Sample from uniform distribution.

Computers can't generate random numbers as they are deterministic but they can produce pseudo random numbers.

Multiplicative Congruential

• involves three integers
• a, b, m
• an initial value x₀, we call the seed
• a sequence of integers is defined as

[math]\displaystyle{ x_{k+1} = (ax_{k} + b) \mod{m} }[/math]

Example: a=13 b=0 m=31 x₀=1 creates a uniform histogram

Matlab code for generating 1000 randoms using the multiplicative congruential method:

a=13;
b=0;
m=31;
x=1;

for ii = 1:1000
    x(ii+1) = mod(a*x(ii)+b,m);
end

Facts about this algorithm:

• the first 30 terms in the sequence are a permutation of integers from 1 to 30 then the sequence repeats itself
• values are between 0 and m-1
• if you divide it by m-1 the result is numbers uniformly distributed in the interval of 0 and 1
• in matlab the values choosen are a=7⁵ b=0 m=2³¹-1 due to a 1988 paper showing they are optimal

Inverse Transform Method

[math]\displaystyle{ P(a\lt x\lt b)=\int_a^{b} f(x) dx }[/math]

[math]\displaystyle{ cdf=F(x)=P(X\lt =x)=\int_{-\infty}^{x} f(x) dx }[/math]

Assume cdf & cdf ^-1 can be found

Theorem:

Take [math]\displaystyle{ U \sim~ \mathrm{Unif}[0, 1] }[/math] and let [math]\displaystyle{ x=F^{-1}(u) }[/math]. Then x has distribution function [math]\displaystyle{ F() }[/math], where [math]\displaystyle{ F(x)=P(X\lt =x) }[/math].

Let [math]\displaystyle{ F^{-1}() }[/math] denote the inverse of [math]\displaystyle{ F() }[/math]. Therefore [math]\displaystyle{ F(x)=u \implies x=F^{-1}(u) }[/math]

Take the the exponential distribution for example

[math]\displaystyle{ \,f(x)={\lambda}e^{-{\lambda}x} }[/math]

[math]\displaystyle{ \,F(x)=\int_0^x {\lambda}e^{-{\lambda}u} du }[/math]

[math]\displaystyle{ \,F(x)=1-e^{-{\lambda}x} }[/math]

Let: [math]\displaystyle{ \,F(x)=y }[/math]

[math]\displaystyle{ \,y=1-e^{-{\lambda}x} }[/math]

[math]\displaystyle{ \,ln(1-y)={-{\lambda}x} }[/math]

[math]\displaystyle{ \,x=\frac{ln(1-y)}{-\lambda} }[/math]

[math]\displaystyle{ \,F^{-1}(x)=\frac{-ln(1-x)}{\lambda} }[/math]

Therefore, to get a exponential distribution from a uniform distribution takes 2 steps.

• Step 1. Draw [math]\displaystyle{ U \sim~ \mathrm{Unif}[0, 1] }[/math]
• Step 2. [math]\displaystyle{ x=\frac{-ln(1-F(u))}{\lambda} }[/math]

Now we just have to show the generated points have a cdf of F(x)

[math]\displaystyle{ \,p(F^{-1}(u)\lt =x) }[/math]

[math]\displaystyle{ \,p(F(F^{-1}(u))\lt =F(x)) }[/math]

[math]\displaystyle{ \,p(u\lt =F(x)) }[/math]

[math]\displaystyle{ \,=F(x) }[/math]

QED

Discrete Case

This same technique can be applied to the discrete case. Generate a discrete random variable [math]\displaystyle{ \,x }[/math] that has probability mass function [math]\displaystyle{ \,p(x=x_i)=p_i }[/math] where [math]\displaystyle{ \,x_0\lt x_1\lt x_2... }[/math] and [math]\displaystyle{ \,\sum p_i=1 }[/math]

• Step 1. Draw [math]\displaystyle{ u \sim~ \mathrm{Unif}[0, 1] }[/math]
• Step 2. [math]\displaystyle{ \,x=x_i }[/math] if [math]\displaystyle{ \,F(x_{i-1})\lt u\lt =F(x_i) }[/math]