Today's Topics

• Markov chain Monte Carlo (MCMC)

• Metropolis-Hastings algorithm

• Gibbs sampler

Problem

Given unnormalized probability density/mass function u(x)∝π(x), estimate high-dimensional integral

I=EX[f(X)]=∫f(x)π(x)dx,X∼π(x).

We have already introduced inverse transform algorithm, rejection sampling, and importance sampling.

But they are not "general" enough.

MCMC

Markov chain Monte Carlo (MCMC) tackles this problem by constructing a Markov chain

X0→X1→⋯→XN,

and approximate the integral as

ˆI=1NN∑i=1f(Xi).

We hope that ˆI→I=EX[f(X)] almost surely as N→∞.

MCMC

Unlike previously introduced Monte Carlo methods, the Markov chain X0→⋯→XN in general contains dependent sample points {Xi}.

The nature of Markov chains implies that Xt+1 should be only based on Xt (with other independent sources of randomness) and not on X0,…,Xt−1.

Therefore, with a given starting point X0, the whole Markov chain can be specified via a conditional distribution T(y|x):

X1∼T(⋅|X0),X2∼T(⋅|X1),…

Transition Kernel

• Such a conditional distribution is typically called a transition kernel

• Denoted as T(x→y)≡T(y|x)

• Think of x as the current position/state, and y the next position/state in the Markov chain

• Our target is to find such a T(x→y) such that ˆI→I  a.s. holds

Regularity Conditions

We first require T(x→y) to satisfy some mild regularity conditions:

• Irreducible: it is possible to get to any state from any state

• Aperiodic: the chain cannot get trapped in cycles

• Positive recurrent: revists every neighborhood infinitely often

These conditions are rather weak and can be easily satisfied.

Invariant/Stationary Distribution

The "real" condition that matters is that π(x) is an invariant/stationary distribution of T(x→y), meaning that

π(y)=∫π(x)T(x→y)dx.

• Combined with the regularity conditions, we have ˆI→I  a.s.

• This property is called the ergodicity of Markov chains

• X0 can be chosen arbitrarily

Detailed Balance

• How to satisfy the stationary distribution condition?

• We typically consider a stronger, but easier-to-verify condition called detailed balance: π(x)T(x→y)=π(y)T(y→x)

• Easy to verify that it implies the stationary distribution condition

• Sufficient but not necessary

Important Algorithms

• Metropolis-Hastings algorithm

• Gibbs sampler

• Metropolis-adjusted Langevin algorithm

• Hamiltonian Monte Carlo

Metropolis-Hastings Algorithm

• Listed as one of the "Top 10 Algorithms" in the 20th century ^[1]

• Foundation of many other MCMC algorithms

• Simple and widely used

Metropolis-Hastings Algorithm

• Choose an initial value X0

• Choose an arbitrary proposal distribution q(y|x)

• Given Xt, propose a new point Y∼q(⋅|Xt)

• Define α=min(1,π(Y)q(Xt|Y)π(Xt)q(Y|Xt))

• With probability α set Xt+1=Y; otherwise set Xt+1=Xt

Some (Important) Remarks

• We do not really need to know π(x)

• The unnormalized density u(x)∝π(x) suffices, since α=min(1,π(Y)q(Xt|Y)π(Xt)q(Y|Xt))=min(1,u(Y)q(Xt|Y)u(Xt)q(Y|Xt))

• q(y|x) can be chosen almost completely freely (but some would be better)

• These properties show the real power of M-H

Some (Important) Remarks

• If the "acceptance" step Xt+1=Y fails, we still need to make Xt+1=Xt (i.e., repeat the previous position once)

• This is different from the rejection sampling

Transition Kernel

• Let α(x,y)=min(1,u(y)q(x|y)u(x)q(y|x))

• Probability of proposing xt+1 and accepting it is q(xt+1|xt)α(xt,xt+1)

• Probability of accepting a move is α(xt)=∫q(y|xt)α(xt,y)dy

• The transition kernel is T(xt→xt+1)=q(xt+1|xt)α(xt,xt+1)+(1−α(xt))δ(xt+1=xt), where δ(⋅) is a Dirac measure

Verifying Detailed Balance

• We need to verify π(xt)T(xt→xt+1)=π(xt+1)T(xt+1→xt)

• First term of LHS is π(xt)q(xt+1|xt)α(xt,xt+1)=π(xt)q(xt+1|xt)⋅min(1,π(xt+1)q(xt|xt+1)π(xt)q(xt+1|xt))=min(π(xt)q(xt+1|xt),π(xt+1)q(xt|xt+1))

• Second term of LHS is π(xt)(1−α(xt))δ(xt=xt+1)=π(xt+1)(1−α(xt+1))δ(xt+1=xt)

Random Walk M-H

• The Metropolis-Hastings algorithm gives full freedom in choosing the proposal distribution q(y|x)

• One commonly-used implementation of the M-H algorithm is the random walk M-H

• It makes q(y|x) a normal distribution N(x,σ2I)

• The acceptance probability simplifies to α=min(1,u(Y)u(Xt))

• σ2 is a hyperparameter

Visualization

https://chi-feng.github.io/mcmc-demo/app.html?algorithm=RandomWalkMH&target=banana

Extensions

• Do we have a better choice than random walk M-H?

• How to design new MCMC algorithm based on M-H?

• We will see an excellent example based on the Langevin algorithm

• A good review of the history of M-H can be found in [2]

• Also introduces a lot of extensions of M-H

Gibbs Sampler

• Suppose the variable X we want to sample can be partitioned as X=(X1,…,Xd)

• Each component Xi can be a scalar or a vector

• The Gibbs sampler has the following iterations:

  • Select an index i from {1,…,d} randomly or deterministically

  • Sample X(k+1)i∼π(xi|X(k)−i)

  • Set X(k+1)=(X(k)1,…,X(k+1)i,…,X(k)d)

Transition Kernel

• Suppose index i is selected with probability ρi

• X(k+1) and X(k) differ only in component Xi

• T(x(k)→x(k+1))=ρiπ(x(k+1)i|x(k)−i)

Detailed Balance

• We need to verify π(x(k))T(x(k)→x(k+1))=π(x(k+1))T(x(k+1)→x(k))

• Note that x(k)−i=x(k+1)−i

• LHS is π(x(k))ρiπ(x(k+1)i|x(k)−i)=π(x(k)−i)π(x(k)i|x(k)−i)ρiπ(x(k+1)i|x(k)−i)=π(x(k+1)−i)π(x(k)i|x(k+1)−i)ρiπ(x(k+1)i|x(k+1)−i)

• By symmetry, equal to RHS

Visualization

https://chi-feng.github.io/mcmc-demo/app.html?algorithm=GibbsSampling&target=banana

References