Understanding Metropolis-Hastings and Gibbs Sampling from First Principles

Learning Objectives

By the end of this lesson, you will:

Understand why MCMC was invented.
Understand what a Markov Chain is.
Understand why sampling solves Bayesian inference.
Understand the logic behind the Metropolis-Hastings acceptance ratio.
Learn Gibbs Sampling from first principles.
Implement both algorithms in Python.
Apply them to healthcare, epidemiology, and supply chain problems.

Why Do We Need MCMC?

Recall Bayes’ theorem:

p(\theta|D) \frac {p(D|\theta)p(\theta)} {\int p(D|\theta)p(\theta)d\theta}

The challenge is usually the denominator:

\int p(D|\theta)p(\theta)d\theta

For most realistic Bayesian models, this integral cannot be evaluated analytically.

Examples include:

Healthcare:

\theta= \text{treatment effect}

Epidemiology:

\theta= \text{disease transmission rate}

Supply Chain:

\theta=\text{future demand parameter}

Instead of solving the integral directly, MCMC asks a different question:

Can we generate samples from the posterior distribution?

If we can generate:

[ \theta^{(1)}, \theta^{(2)}, \dots, \theta^{(N)} ]

from

p(\theta|D)

then almost every Bayesian quantity can be estimated using simple averages.

The Big Idea Behind MCMC

Suppose we want the posterior mean.

Mathematically:

E(\theta|D)=\int \theta p(\theta|D)d\theta

If we have posterior samples:

\theta^{(1)}, \theta^{(2)}, \dots, \theta^{(N)}

then:

E(\theta|D) \approx \frac1N \sum_{i=1}^{N} \theta^{(i)}

Similarly:

Posterior variance:

Var(\theta|D) \approx \frac1N \sum_{i=1}^{N} (\theta^{(i)}-\bar{\theta})^2

Posterior probability:

P(\theta>0\mid D) \approx \frac{\#(\theta^{(i)}>0)}{N}

The entire Bayesian problem becomes:

Generate samples from the posterior distribution.

What Is a Markov Chain?

A Markov Chain is a sequence:

[ X_1,X_2,X_3,\dots ]

satisfying:

P(X_{t+1}|X_t,X_{t-1},\dots)= P(X_{t+1}|X_t)

The future depends only on the present.

It does not depend on the entire history.

Healthcare example:

Suppose a patient can be:

Healthy
Sick

The patient’s health tomorrow depends mostly on today’s condition.

This is a Markov process.

The Goal of MCMC

Instead of sampling directly from

p(\theta|D)

which is difficult,

we construct a Markov chain whose long-run distribution is

p(\theta|D)

Then after running the chain long enough:

\theta^{(t)} \sim p(\theta|D)

This is the central idea behind MCMC.

Metropolis-Hastings Algorithm

Metropolis-Hastings is the most widely used MCMC algorithm.

The algorithm repeatedly:

Proposes a move.
Decides whether to accept it.
Repeats.

Over time the chain spends the correct proportion of time in each region of the posterior.

Step 1: Start Somewhere

Suppose our current value is:

\theta_t

For example:

\theta_t=2

The starting point is usually arbitrary.

Step 2: Generate a Proposal

We generate a candidate value:

\theta^* \sim q(\theta^*|\theta_t)

A common choice is:

\theta^* \sim N(\theta_t,\sigma^2)

This means we randomly explore nearby values.

Suppose:

\theta^*=2.4

Step 3: Decide Whether To Move

This is where the brilliance of Metropolis-Hastings appears.

We compute:

r=\frac {p(\theta^*|D)} {p(\theta_t|D)}

for symmetric proposal distributions.

Accept with probability:

\min(1,r)

Why Do We Use This Ratio?

This is one of the most important questions in Bayesian computation.

Suppose:

Current location:

p(\theta_t|D)=0.10

Proposed location:

p(\theta^*|D)=0.50

Then:

[ r= \frac{0.50}{0.10}= 5 ]

Since:

r>1

we automatically accept.

Why?

Because the proposal lies in a region that is five times more probable under the posterior.

The chain should spend more time there.

Why Not Always Reject Worse Points?

Suppose:

Current location: 0.50
Proposed location: 0.25

Then:

r= \frac{0.25}{0.50} =0.50

Accept with probability: 50%

At first this seems strange.

Why move to a less likely location?

Because if we only accepted better moves, the chain could become trapped.

Imagine a posterior with two peaks.

If the chain reaches one peak, it may never leave.

Allowing occasional downhill moves enables exploration of the entire posterior distribution.

The Real Metropolis-Hastings Ratio

For asymmetric proposal distributions:

r=\frac {p(\theta^*|D) q(\theta_t|\theta^*)} {p(\theta_t|D) q(\theta^*|\theta_t)}

The extra terms compensate for proposal bias.

For symmetric proposals:

q(\theta_t|\theta^*)=q(\theta^*|\theta_t)

and they cancel.

Why Is MCMC So Powerful?

Notice:

p(\theta|D)=\frac {p(D|\theta)p(\theta)} {p(D)}

Substituting into the acceptance ratio:

r=\frac {\frac{p(D|\theta^*)p(\theta^*)}{p(D)}} {\frac{p(D|\theta_t)p(\theta_t)}{p(D)}}

The difficult term: $p(D)$ cancels completely.

This is the key reason MCMC became so successful.

The hardest integral in Bayesian statistics never needs to be computed.

Example 1: Healthcare — Drug Effectiveness

Suppose:

\theta=\text{average reduction in blood pressure}

Posterior:

p(\theta|D)

Goal:

P(\theta>0|D)

Generate 50,000 Metropolis-Hastings samples.

Estimate:

P(\theta>0\mid D) \approx \frac{\#(\theta^{(i)}>0)}{50000}

This gives the probability that the treatment truly improves patient outcomes.

Example 2: Epidemiology

Suppose:

[ \beta=\text{infection transmission rate} ]

Posterior:

p(\beta|D)

Goal:

P(\beta>1)

Values above one indicate an expanding epidemic.

Metropolis-Hastings samples provide the answer immediately.

Example 3: Supply Chain Forecasting

Suppose:

\lambda=\text{monthly demand rate}

Posterior:

p(\lambda|Sales)

Generate posterior samples:

\lambda^{(1)}, \dots, \lambda^{(50000)}

Then estimate:

P(Demand>Inventory)

using simulation.

Python Example: Metropolis-Hastings

import numpy as np

def target(x):
    return np.exp(-x**2/2)

n = 10000

theta = 0

samples = []

for i in range(n):

    proposal = np.random.normal(
        theta,
        1
    )

    r = target(proposal) / target(theta)

    if np.random.rand() < min(1, r):
        theta = proposal

    samples.append(theta)

samples = np.array(samples)

import numpy as np

def target(x):
    return np.exp(-x**2/2)

n = 10000

theta = 0

samples = []

for i in range(n):

    proposal = np.random.normal(
        theta,
        1
    )

    r = target(proposal) / target(theta)

    if np.random.rand() < min(1, r):
        theta = proposal

    samples.append(theta)

samples = np.array(samples)

Understanding Gibbs Sampling

Metropolis-Hastings is very general.

Gibbs Sampling is a special case that becomes extremely efficient when conditional distributions are available.

The Core Idea

Suppose we have two parameters $\theta_1$ and $\theta_2$

The joint posterior:

p(\theta_1,\theta_2|D)

may be difficult.

However, suppose we know:

p(\theta_1|\theta_2,D)

and

p(\theta_2|\theta_1,D)

These conditional distributions may be easy to sample from.

Gibbs Sampling Algorithm

Start with:

(\theta_1^{(0)},\theta_2^{(0)})

Update:

\theta_1^{(1)} \sim p(\theta_1|\theta_2^{(0)},D)

Then:

\theta_2^{(1)} \sim p(\theta_2|\theta_1^{(1)},D)

Repeat indefinitely.

Why Does Gibbs Sampling Work?

Imagine moving around a mountain.

Metropolis-Hastings jumps randomly.

Gibbs moves one coordinate at a time.

Update

\theta_1

Then update:

\theta_2

Then update:

\theta_1

again.

Eventually the chain explores the entire posterior.

Why Is There No Acceptance Ratio?

Because every draw already comes from the correct conditional distribution.

Every proposed value is automatically valid.

Mathematically:

r=1

for every iteration.

Nothing is rejected.

Example 1: Blood Pressure Model

Parameters:

\mu=\text{population mean blood pressure} , \\ \sigma^2=\text{population variance}

Posterior:

p(\mu,\sigma^2|D)

Update:

p(\mu|\sigma^2,D)

then

p(\sigma^2|\mu,D)

Example 2: Hospital Readmission Model

Parameters:

\beta_0

and

\beta_1

Posterior:

p(\beta_0,\beta_1|D)

Update:

\beta_0|\beta_1,D

then

\beta_1|\beta_0,D

Example 3: Supply Chain Forecasting

Parameters:

\lambda=\text{demand rate}

\alpha=\text{seasonality effect}

Posterior:

p(\lambda,\alpha|Sales)

Update each parameter conditionally.

Generate future demand forecasts using posterior samples.

Python Example: Gibbs Sampling

import numpy as np

n = 10000

x = 0
y = 0

samples = []

for i in range(n):

    x = np.random.normal(
        y,
        1
    )

    y = np.random.normal(
        x,
        1
    )

    samples.append([x, y])

samples = np.array(samples)

import numpy as np

n = 10000

x = 0
y = 0

samples = []

for i in range(n):

    x = np.random.normal(
        y,
        1
    )

    y = np.random.normal(
        x,
        1
    )

    samples.append([x, y])

samples = np.array(samples)

Metropolis-Hastings vs Gibbs Sampling

Feature	Metropolis-Hastings	Gibbs Sampling
Requires conditional distributions	No	Yes
Acceptance ratio	Yes	No
Every proposal accepted	No	Yes
General-purpose	Yes	No
Easier to apply	Yes	Sometimes
More efficient when conditionals known	No	Yes

Key Takeaways

MCMC solves Bayesian inference by generating samples from the posterior distribution.

Metropolis-Hastings uses the acceptance ratio:

r=\frac {p(\theta^*|D)} {p(\theta_t|D)}

to ensure the chain visits regions according to their posterior probability.

The normalizing constant:

p(D)

cancels from the acceptance ratio, making Bayesian inference feasible.

Gibbs Sampling is a special case of MCMC where we sample directly from conditional distributions.

Because each sample already comes from the correct conditional distribution:

r=1

and every proposal is accepted.

Together, Metropolis-Hastings and Gibbs Sampling form the foundation of modern Bayesian computation and lead directly to advanced methods such as Hamiltonian Monte Carlo, NUTS, Sequential Monte Carlo, Bayesian Nonparametrics, and Bayesian Machine Learning.

References

Gelman, Carlin, Stern, Dunson, Vehtari, and Rubin. Bayesian Data Analysis.
Robert and Casella. Monte Carlo Statistical Methods.
Brooks, Gelman, Jones, and Meng. Handbook of Markov Chain Monte Carlo.
Murphy. Machine Learning: A Probabilistic Perspective.

nerd-ish

Leave a ReplyCancel reply

Measure Theory Lesson 30: Radon Measures

Measure Theory Lesson 52: Bott Periodicity — The Miracle Behind K-Theory

Measure Theory Lesson 51: K-Theory for Operator Algebras — Measuring the Shape of a Space Without Points

Lesson 2: Markov Chain Monte Carlo (MCMC)

Understanding Metropolis-Hastings and Gibbs Sampling from First Principles

Learning Objectives

Why Do We Need MCMC?

The Big Idea Behind MCMC

What Is a Markov Chain?

The Goal of MCMC

Metropolis-Hastings Algorithm

Step 1: Start Somewhere

Step 2: Generate a Proposal

Step 3: Decide Whether To Move

Why Do We Use This Ratio?

Why Not Always Reject Worse Points?

The Real Metropolis-Hastings Ratio

Why Is MCMC So Powerful?

Example 1: Healthcare — Drug Effectiveness

Example 2: Epidemiology

Example 3: Supply Chain Forecasting

Python Example: Metropolis-Hastings

Understanding Gibbs Sampling

The Core Idea

Gibbs Sampling Algorithm

Why Does Gibbs Sampling Work?

Why Is There No Acceptance Ratio?

Example 1: Blood Pressure Model

Example 2: Hospital Readmission Model

Example 3: Supply Chain Forecasting

Python Example: Gibbs Sampling

Metropolis-Hastings vs Gibbs Sampling

Key Takeaways

References

Share this:

Like this:

Related posts:

Leave a ReplyCancel reply

Discover more from nerd-ish