Some notes on Markov Chain Monte Carlo: Metropolis-Hastings and Gibbs Steps for Exploring Posterior Distributions

1 Practical Bayesian Inference

In the Bayesian paradigm we build a joint model for data \underline{y} = (y_{1}, …, y_{N}), and parameters τ:

p(\underline{y}, τ) = p(\underline{y}|τ) · p(τ)

and we always want to know something about the posterior density

\begin{aligned}p(τ|\underline{y}) & = & \frac{p(\underline{y}|τ) · p(τ)}{∫ p(\underline{y}|t) · p(t) dt} \\ & ∝ & p(\underline{y}|τ) · p(τ) = p(\underline{y}, τ) \end{aligned}

For example,

• Posterior mode: argmax_{τ} p(τ|\underline{y})

• Posterior mean: E[τ|\underline{y}] = ∫ τ p(τ|\underline{y}) dτ (or some other risk-minimizer)

• Credible interval: Find a set A such that p(τ ∈ A|\underline{y}) = 0.95

• Graph or “shape” of p(τ|\underline{y})

\underline{All} of these operations require that you know a lot about the shape of p(τ|\underline{y}). How do you learn enough about the shape of p(τ|\underline{y})?

• If p(τ) is conjugate to p(\underline{y}|τ), then p(τ|\underline{y}) will belong to the same family as p(τ).

• If not, then you have a numerical integration problem. These problems are difficult because they depend on the “interesting” features of p(τ|\underline{y}).

– Difficult but do-able if p(τ|\underline{y}) is smooth and low-dimensional.

– Even if p(τ|\underline{y}) is smooth, as d = dim(τ) grows, space becomes too large for deterministic methods. Curse of dimensionality: For n “evaluations” of p(τ|\underline{y}), learn at rate O(1/n^{k/d}).

• Monte Carlo (MC) methods have long been used to help with the curse. With a random sample τ^{(1)}, τ^{(2)}, …, τ^{(n)}, can plot histograms (shape), compute averages (integration), etc.

– Improves convergence to O(1/\sqrt{n})

– Less smoothing of interesting features as n grows

– Curse returns: how do you sample from p(τ|\underline{y})?

2 MCMC: The Cliff’s Notes

The problem:

• Learn about π(τ) = π(θ, β, λ), where τ = (θ, β, λ) is some high-dimensional set of variables (parameters).

The essential idea:

• Define a (stationary) Markov chain M_{0}, M_{1}, M_{2}, … with states M_{k} = (θ^{k}, β^{k}, λ^{k}).

• Simulate the M_{k}’s; under regularity conditions (e.g., Tierney, 1994), M_{k} will converge in distribution to a stationary distribution π(θ, β, λ):

π(θ^{1}, β^{1}, λ^{1}) = ∫ P [M_{k+1} = (θ^{1}, β^{1}, λ^{1}) | M_{k} = (θ^{0}, β^{0}, λ^{0})] × π(θ^{0}, β^{0}, λ^{0}) d(θ^{0}, β^{0}, λ^{0})

• For Bayes, design chain so that π(θ, β, λ) turns out to be the posterior, p(θ, β, λ|\underline{y}).

Typical MCMC paradigm:

• Write p(θ, β, λ|\underline{y}) = p(τ_{1}, τ_{2}, …, τ_{M}|\underline{y}), where (τ_{1}, …, τ_{M}) is a disjoint partition of the original model parameters (θ, β, λ).

• General theory of MCMC (e.g., Tierney, 1994; Chib and Greenberg, 1995): construct state M_{k} = (τ_{1}^{(k)}, …, τ_{M}^{(k)}), by sampling each τ_{m} from its “complete conditional” distribution:

To step from M_{k−1} = (τ_{1}^{(k−1)}, …, τ_{M}^{(k−1)}) to M_{k} = (τ_{1}^{(k)}, …, τ_{M}^{(k)}):

1. τ_{1}^{(k)} ∼ p(τ_{1}|τ_{2}^{(k−1)}, …, τ_{M}^{(k−1)}, \underline{y});

2. τ_{2}^{(k)} ∼ p(τ_{2}|τ_{1}^{(k)}, τ_{3}^{(k−1)}, …, τ_{M}^{(k−1)}, \underline{y});

3. τ_{3}^{(k)} ∼ p(τ_{3}|τ_{1}^{(k)}, τ_{2}^{(k)}, τ_{4}^{(k−1)}, …, τ_{M}^{(k−1)}, \underline{y});

⋮  ⋮  ⋮  ⋮

M. τ_{M}^{(k)} ∼ p(τ_{M}|τ_{1}^{(k)}, τ_{2}^{(k)}, τ_{3}^{(k)}, …, τ_{M−1}^{(k)}, \underline{y},);

This is the basic idea behind all so-called “Gibbs samplers.”

Markov Chain Monte Carlo (MCMC) vs. MC:

• Replace draws from p(τ|\underline{y}) with a sequence of low-dimensional draws from p_{1}(τ_{1}), p_{2}(τ_{2}), etc.

• For Bayes, often can exploit “partial” conjugacy to simplify problem of drawing from p(τ_{(k)}).

• Costs: dependent rather than iid samples, approximate rather than exact sampling of p(τ|\underline{y}).

3 Implementation Details

Sampling schemes

Once a set of complete conditionals is decided upon:

• If you can sample τ_{m} directly from its complete conditional distribution, do so—this is a \underline{Gibbs step};

• Otherwise perform a \underline{Metropolis step}:

– Sample a “proposal value” τ_{m}^{∗} from any convenient proposal distribution q_{k}(τ_{m}|τ_{m}^{(k−1)})

– Compute the “acceptance probability” α_{m}^{∗} from q_{k}(τ_{m}|τ_{m}^{(k−1)}) and the complete conditional p(τ_{m}|rest) ≡ p(τ_{m} | τ_{1}^{(k)}, …, τ_{m−1}^{(k)}, τ_{m+1}^{(k−1)}, …, τ_{M}^{(k−1)}, \underline{y}):

α_{n}^{∗} = min \{\frac{p(τ_{n}^{∗}|rest)q_{k}(τ_{n}^{(k−1)}|τ_{n}^{∗})}{p(τ_{n}^{(k−1)}|rest)q_{k}(τ_{n}^{∗}|τ_{n}^{(k−1)})}, 1\};

– Accept τ_{n}^{(k)} = τ_{n}^{∗} with probability α_{n}^{∗}; otherwise set τ_{n}^{(k)} = τ_{n}^{(k−1)}.

This is often called “Metropolis-Hastings within Gibbs”.

Identifying the complete conditionals

Suppose the θ’s, β’s and λ’s have prior distributions p(θ|λ), p(β) and p(λ). Then the “complete conditional” for θ, for example, is

p(θ|rest) = p(θ|\underline{y}, β, λ) = \frac{p(\underline{y}, θ, β, λ)}{∫ p(\underline{y}, t, β, λ) dt}
= \frac{p(\underline{y}|θ, β) p(θ|λ) p(β) p(λ)}{∫ p(\underline{y}|t, β) p(t|λ) p(β) p(λ) dt}
∝ p(\underline{y}|θ, β) p(θ|λ)

Key observation: The shape of p(θ|\underline{y}, β, λ) is determined by just the parts of the likelihood that depend explicitly on θ.

Convergence

M_{0}, M_{1}, M_{2}, …, M_{k}, “Burn-in” segment

\underbrace{M_{k+1}, …, M_{k+n}} Usable MCMC sample

How large should k be?

• Use acf() in Splus/R to check when Corr(M_{ℓ}, M_{ℓ+k}) small;

• Use time series plots to see when the chain has “stabilized”.

• Run several chains from different starting points; let k be so large that (within-chain variation) < (between chain variation).

• CODA subroutine package for Splus/R offers a menu of such checks (www.mrc-bsu.cam.ac.uk/bugs/).

4 Example: The hierarchical beta-binomial model

Non-hierarchical version

Basic beta-binomial model, for examinees i = 1, …, N, each getting y_{i} of n_{i} questions right. [\underline{y} = (y_{1}, …, y_{N})]

1. y_{i}|θ_{i} ∼ Bin(θ_{i}, n_{i});

2. θ_{i}|α, β ∼ Beta(α, β) Full model p(\underline{y}, \underline{θ}|α, β) is

∏_{i=1}^{N} (\begin{matrix}n_{i} \\ y_{i} \end{matrix}) θ_{i}^{y_{i}}(1 − θ_{i})^{n_{i}−y_{i}} ∏_{i=1}^{N} \frac{Γ(α + β)}{Γ(α)Γ(β)}θ_{i}^{α−1}(1 − θ_{i})^{β−1}

Use the “key observation” to calculate

p(θ_{i}|y_{i}, α, β) ∝ θ_{i}^{α+y_{i}−1}(1 − θ_{i})^{β+n_{i}−y_{i}−1}

i.e. θ_{i}|y_{i} ∼ Beta(α + y_{i}, β + n_{i} − y_{i}). Choose fixed values for α, β, then use Beta(α + y_{i}, β + n_{i} − y_{i}) for posterior mode, posterior mean, credible interval for θ, etc. No need for anything fancy yet.

Hierarchical version

What if we want to estimate α, β (estimate the shape of the population / latent distribution)? Add a third modeling assumption:

3. α, β ∼ p(α, β)

Now the full model p(\underline{y}, \underline{θ}, α, β) is

∏_{i=1}^{N} (\begin{matrix}n_{i} \\ y_{i} \end{matrix}) θ_{i}^{y_{i}}(1 − θ_{i})^{n_{i}−y_{i}} ∏_{i=1}^{N} \frac{Γ(α + β)}{Γ(α)Γ(β)}θ_{i}^{α−1}(1 − θ_{i})^{β−1}p(α, β)

Can again use the “key observation” to calculate

\begin{aligned}p(θ_{i}|y_{i}, α, β) & = & Beta(θ_{i}|α + y_{i}, β + n_{i} − y_{i}) \\ p(α, β|\underline{y}) & ∝ & p(α, β) ∏_{i=1}^{N} \frac{Γ(α + β)}{Γ(α)Γ(β)} \frac{Γ(α + y_{i})Γ(β + n_{i} − y_{i})}{Γ(α + β + n_{i})} \end{aligned}

An MCMC solution

From the full model p(\underline{y}, \underline{θ}, α, β) is

∏_{i=1}^{N} (\begin{matrix}n_{i} \\ y_{i} \end{matrix}) θ_{i}^{y_{i}}(1 − θ_{i})^{n_{i}−y_{i}} ∏_{i=1}^{N} \frac{Γ(α + β)}{Γ(α)Γ(β)}θ_{i}^{α−1}(1 − θ_{i})^{β−1}p(α, β)

the “key observation” gives the complete conditionals

p(θ_{i}|rest) = Beta(θ_{i}|α + y_{i}, β + n_{i} − y_{i})
p(α|rest) ∝ [\frac{Γ(α + β)}{Γ(α)}]^{N} ∏_{i=1}^{N} θ_{i}^{α}p(α, β)

p(β|rest) ∝ [\frac{Γ(α + β)}{Γ(β)}]^{N} ∏_{i=1}^{N}(1 − θ_{i})^{β}p(α, β)

• Usually take p(α, β) = p(α)p(β); in this case we could take independent Gamma or Inv-Gamma densities, e.g.

• Canned routines sample θ_{i} ∼ Beta(⋯)

• Can set up Metropolis steps for α and β using Normal proposal draws (“random walk M-H”). Normal variances are “tuning parameters”.