Lecture III

An Introduction to Markov Chain Monte Carlo

Lecture III

Abstract

Connection with Statistics

To sample from an arbitrary random vector X Î R^p with pdf f just turn f into a canonical distribution by choosing T and Z(T) and defining an energy function by,

Plain Metropolis with Symmetric Proposal Distribution

1. Propose a candidate state y from a proposal distribution, q(y|x) that may depend on the current x. This proposal distribution is assumed to be symmetric in x and y, i.e., for plain vanilla Metropolis,

   q(y|x) = q(x|y)

   this can be easily implemented, e.g., take

   y = x + e

   with e a random vector radially symmetric about 0. Typical choices for the distribution of e are multivariate Cauchy, Gaussian, or T distribution with a few degrees of freedom.

2. If the new proposed state y uses less energy than the current state x then go there with probability one. If the new state is more expensive, in terms of energy, than the one we are currently on, then test your luck and go there with a probability exponentially decreasing in the difference of energy. More formally the acceptance probability is,

   a(x,y) = exp ì í î - (E(y) - E(x))+ T ü ý þ = min { 1, f(y) f(x) }

The Plain Vanilla Metropolis

{

X₀ ¬ x₀ (arbitrary)
for t = 0,1,2,¼

{

y ¬  sample from q(·|x_t)
u ¬ unif(0,1)

if u £ min{1,f(y)/f(x_t)} then X_t+1 ¬ y
else X_t+1 ¬ x_t

}

Example: Javascript implementation of plain vanilla Metropolis for sampling from N(0,1) with

Solution by Haihong Li

Metropolis-Hastings

which is true, since when the min. on one side is 1, the min. on the other side isn't and the denominator cancels with the term outside the parenthesis producing the equality of both sides.

As always, detailed balance is enough to assure that f is the stationary distribution for the chain. But detailed balance is not necessary. It is possible for a chain to have f as its stationary distribution without p(y|x)f(x) = p(x|y)f(y) for all x,y. All that it is needed for stationarity is, (its definition)

Knowing that f is a stationary distribution for a Markov Chain with transition kernel p(x|y) is not enough to warrant that the chain will eventually sample states according to the density f. It only says that if at some time t the density of X_t is f then it will remain f for all subsequent times since the law of X_t+1 is Tf = f. From the classical fix point theorem of functional analysis we know that if the f's belong to a complete metric space, with distance function d, and the operator T is contractive, i.e.

then the iterates f₀, T(f₀), T(T(f₀)),¼ converge geometrically in the metric of the space to the unique fix point for T, provided that the initial f₀ is not too far from the fix point. Notice that,

Thus, the Markov Chain is expected to have the fix point as its long-run distribution. As we will see later the conditions of the classic fix point theorem are too restrictive. Under mild regularity conditions on the kernel (e.g. when T is only contractive on the average) a large class of Markov Chains will be geometrically ergodic. More on this theme when we look at convergence theorems...

Component-Wise Metropolis

If J Ì I we write x^-J the set of all coordinates except those in the set J, i.e. x^-J = x^I\J. Let I₁,I₂,¼,I_m form a partition of I so that x gets separated into m groups,

Start the Metropolis-Hastings algorithm from somewhere and suppose that at time t we are visiting state x, then at time t+1 we either remain at x or we go to y by modifying only one of the m components of x, i.e.,

The value of k can be chosen at random among {1,2,¼,m} or one after another cyclically, in which case to have a homogeneous Markov chain we need to define one iteration only after a full sweep over the set of m components is completed. Choosing the components at random eliminates this problem but may end up, by chance, neglecting some of the components for a long time. An intermediate solution is to randomize the order and then update the components in that order. In all these cases the acceptance probability is the usual Metropolis-Hastings formula but when x differs from y in the k component we have,

where q_k is the proposal function q for updating component k and ratio f(y)/f(x) is written as the ratio of full conditionals, where,

By writing the acceptance probability in this way we can see that the Markov Chain associated to the sequence of updates of component k, with all the other components fix at x^-I_k, has the full conditional (defined above) as stationary distribution. Eventually all the components end up sampling from the correct conditionals.

Gibbs Sampler

In this case the acceptance probability becomes 1 and the Gibbs sampler always accepts the proposal candidates. Gibbs sampling is just a special case of Metropolis. This is the justification to the popular statement: to sample from a joint distribution just sample repeatedly from its one dimensional conditionals given whatever you've seen at the time.

The Plain Vanilla Gibbs

{

X₀ ¬ x₀ (arbitrary)
for t = 0,1,2,¼

{

k ¬ unif. on {1,2,¼,p}
y^k ¬ sample from f(x^k|x_t^-{k})
X_t+1 ¬ {y^k,x_t^-{k}}

}

Example: Javascript implementation of Gibbs sampler for generating samples from (X,Y) with one dimensional conditionals given by: