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 parameters: integers [math]\displaystyle{ \,a, b, m }[/math], and an initial value [math]\displaystyle{ \,x_0 }[/math] we call the seed
• a sequence of integers is defined as

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

Example: [math]\displaystyle{ \,a=13, b=0, m=31, x_0=1 }[/math] creates a uniform histogram.

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

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

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

MATLAB code for displaying the values of x generated:

x;

MATLAB code for plotting the histogram of x:

hist(x);

Facts about this algorithm:

• In this example, 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 [math]\displaystyle{ m-1 }[/math].
• Dividing the numbers by [math]\displaystyle{ m-1 }[/math] yields numbers in the interval [math]\displaystyle{ [0,1) }[/math].
• MATLAB's rand function once used this algorithm with [math]\displaystyle{ a=7^5, b=0, m=2^{31}-1 }[/math], for reasons described in Park and Miller's 1988 paper "Random Number Generators: Good Ones are Hard to Find" (available online).

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]