MCMC and Bayesian Data Analysis(PPT在文件模块)

How to generate a sample from $p(x)$?

Let's first see how Matlab samples from a $p(x)$. In Matlab, there are several common probability distributions.

Try univariate Gaussian distribution

p= normpdf(linspace(xmin , xmax , num_of_points) , mean, standard_deviation);%PDF

c= normcdf(linspace(xmin , xmax , num_of_points) , mean, standard_deviation);%CDF

y=normrnd(mean, standard_deviation,number_of_samples, 1);%Random Number Generating Method

Try PDF:

x=linespace(-1,1,1000);

p=normpdf(x ,0, 1);

plot(x,p);

Note: linespace returns a vector which is usually accessed like this

x(1)%the first elem, not x(0)

x(:)

x(1,:)

Try RSM:

y=normrnd(0, 1,100, 1);%试试采10000个样本

hist( y , 20 );%20 bars

Try univariate uniform distribution

p= unifpdf(linspace(xmin , xmax , num_of_points) , a,b);%PDF

c= unifcdf(linspace(xmin , xmax , num_of_points) , a,b);%CDF

y=unifrnd(a,b,number_of_samples, 1);%RNG

Try PDF:

x=linespace(-10,10,1000);

p= unifpdf(x ,-5,5);

plot(x,p);

Try RSM:

y=unifrnd(-5, 5,100, 1);%试试采10000个样本

hist( y , 20 );%20 bars

Matlab provides random number generating functions for some standard $p(x)$, it doesn't provide us sampling functions for a general $p(x)$. Here I show some common sampling methods.

Inverse Transform Sampling(ITS)

with descret variables

This method generates random numbers from any probability distribution given the inverse of its cumulative distribution function . The idea is to sample uniformly distributed random numbers (between 0 and 1) and then transform these values using the inverse cumulative distribution function(InvCDF)(which can be descret or continous). If the InvCDF is descrete, then the ITS method just requires a table lookup, like shown in Table 1.

Table 1. Probability of digits observed in human random digit generation experiment

There is a method called randsample in Matlab that can implement the sampling process using the Table 1. See the code below.

%Note: The randsample doesn't defaultly exist in Octave-core package, install statistic package from http://octave.sourceforge.net/statistics/ before using randsample.

% probabilities for each digit

theta=[0.000; ... % digit 0

0.100; ... % digit 1

0.090; ... % digit 2

0.095; ... % digit 3

0.200; ... % digit 4

0.175; ... % digit 5

0.190; ... % digit 6

0.050; ... % digit 7

0.100; ... % digit 8

0.000];

seed = 1; rand( 'state' , seed );% fix the random number generator

K = 10000;% let's say we draw K random values

digitset = 0:9;

Y = randsample(digitset,K,true,theta);

figure( 1 ); clf;

counts = hist( Y , digitset );

bar( digitset , counts , 'k' );

xlim([-0.5 9.5]);

xlabel( 'Digit' );

ylabel( 'Frequency' );

title( 'Distribution of simulated draws of human digit generator' );

pause;

Instead of using the built-in functions such as randsample or mnrnd, it is helpful to consider how to implement the underlying sampling algorithm using the inverse transform method which is:

(1) Calculate $F(X)$.

(2) Sample u from Uniform(0,1).

(3) Get a sample $x^{i}$ of $P(X)$, which is $F(u)^{-1}$.

(4) Repeat (2) and (3) until we get enough samples.

Note: For discrete distributions, $F(X)^{-1}$ is discrete, the way to get a sample $x^{i}$ is illustrated below where $u=0.8,~x^{i}=6$ .

with continuous variables

This can be done with the following procedure:

(1) Draw U ∼ Uniform(0, 1).

(2) Set $X=F(U)^{-1}$

(3) Repeat

For example, we want to sample random numbers from the exponential distribution where  its CDF is F (x|λ) = 1 − exp(−x/λ) . Then $F(u|\gamma)^{-1}=-log(1-u)\gamma$. Therefore replace $F(U)^{-1}$ with $F(u|\gamma)^{-1}$.

p=-log(1-unifrnd(0,1,10000,1))*2;

hist(p,30);

Reject Sampling

Applied situation: impossible/difficult to compute CDF of $P(X)$.

Advantage: unlike MCMC, it doesn't require of any “burn-in” period, i.e., all samples obtained during sampling can immediately be used as samples from the target distribution $p(\theta)$.

Based on the Figure above, the method is:

(1) Choose a proposal distribution q(θ) that is easy to sample from.

(2) Find a constant c such that cq(θ) ≥ p(θ) for all θ.

(3) Draw a proposal θ from q(θ).

(4) Draw a u from Uniform[0, cq(θ)].

(5) Reject the proposal if u > p(θ), accept otherwise. Actually, since u is sampled from Uniform[0, cq(θ)], it is equal to state like this " Reject if $u\in[p(\theta),cq(\theta)]$, accept otherwise".

(6) Repeat steps 3, 4, and 5 until desired number of samples is reached; each accepted sample $\theta$ is a draw from p(θ).

For example

then the code is

k=100000;%draw k samples

c=2;

theta_vec=unifrnd(0,1,k,1)%gen a proposal vector from q($\theta$)

cq_vec=c*unifpdf(theta_vec);%cq(theta) vector

p_vec=2*theta_vec;%p(theta) vector

u_vec=[];

for cq=cq_vec

u_vec=[u_vec;unifrnd(0,cq)];

end

r=theta_vec.*(u_vec<p_vec);

r(r==0)=[];%remove the “0” elements

hist(r,20);

MCMC Sampling

Before getting to know MCMC sampling, we first get to know Monte Carlo Integration and Markov Chain.

For example:

%Implement the Markov Chain involving x under Beta(200(0.9x^((t-1))+0.05),200(1-0.9x^((t-1)-0.05))

fa=inline('','x')%parameter a for beta

fb=inline('200*(1-0.9*x-0.05)','x');%parameter b for beta

no4mc=4;%4 markove chains

states=unifrnd(0,1,1,no4mc);%initial states

N=1000;%200 samples drawn from 4 chains

X=states;

for i=1:N

states=betarnd(fa(states),fb(states));

X=[X;states];

end;

plot(X);

pause;

Metroplis Sampling