Stochastic Series - Pt 8 - Metropolis-Hastings algorithm

The Metropolis-Hastings algorithm was quite hard for me to understand. Not because the math was complicated, but because I couldn't understand why it even existed in the first place. In this guide we are going to look at why it's actually a helpful and interesting algorithm.

Main purpose

"The purpose of Metropolis-Hastings is to sample from a complex probability distribution without having to calculate every possibility."

There is A LOT to unpack here, so first let's break down these concepts in the right order.

Efficiently

Doing stuff efficiently in math is often about not doing unnecessary work, like calculating the probabilities of every possible state when we only need to evaluate the ones we visit.

MH can help us avoid that work, although it isn't automatically the fastest way to sample from every distribution.

Sample

A sample is one value or observation drawn from a probability distribution. Imagine you have a box with an unknown number of red and blue balls. You close your eyes, pull out one ball, write down the color, and put the ball back.

If you keep doing this, you end up with a collection of samples. For example, you might observe:

Blue balls: 42%
Red balls: 58%

The more samples you draw, the closer these proportions will generally get to the actual distribution of colors in the box.

Complex probability distribution

A complex probability distribution is a distribution where calculating or exploring the probabilities of all possible outcomes is computationally expensive or difficult.

An expensive operation can be defined as something that takes a lot of computational resources, like CPU and memory.

The important thing about MH is that we don't have to calculate the probability of every state. We do, however, need a way to calculate the relative probability of the states we visit.

Back to Metropolis-Hastings

So we use the MH algorithm to construct a Markov chain that follows a target probability distribution, giving us a sequence of samples.

We are not mimicking the transitions of another Markov chain. We are constructing our own chain so that its stationary distribution is the one we want.

How MH does it

Imagine we have two different states: A and B. We don't know the transition matrix (P); we only know the stationary distribution (π) we want.

P = ?

π = (\frac{2}{3}, \frac{1}{3})

So we want to construct a Markov chain where it spends approximately 66.7% of the time in state A and 33.3% in state B in the long run.

To make the numbers easier to follow, imagine 18 copies of the chain: 12 currently in A and 6 in B. This matches our target distribution.

Let's start by giving both states a 50% probability of proposing a move to the other state, and a 50% probability of proposing to stay where they are.

Direction	Starting	Probability	Expected volume
A -> B	12	0.5	6
B -> A	6	0.5	3

In math:

q (A, B) = 0.5

q (B, A) = 0.5

Looking at the table, we can see that each proposal to move between the states has the same probability: 50%.

But the expected volume moving from A to B is 6, while only 3 move from B to A. That doesn't preserve our desired distribution. We need a way to alter the actual transitions, so we add two concepts:

Proposal probability (q)
Acceptance probability (α)

The proposal probability tells us how likely we are to suggest moving from one state to another.

The acceptance probability tells us how likely we are to accept that suggestion.

Instead of accepting every proposed transition, MH gives us a chance of rejecting it. If we reject a proposal, we stay in the current state.

The Metropolis-Hastings formula

The full formula is:

α (x, y) = min (1, \frac{π (y) q (y, x)}{π (x) q (x, y)})

Where:

x = Current state
y = Proposed next state
π(x) = Target probability of state x
π(y) = Target probability of state y
q(x,y) = Probability of proposing y from x
q(y,x) = Probability of proposing x from y
α(x,y) = Probability of accepting the proposed transition

Let's calculate the acceptance probability from A to B:

α (A, B) = min (1, \frac{(1 / 3) \times 0.5}{(2 / 3) \times 0.5}) = 0.5

And from B to A:

α (B, A) = min (1, \frac{(2 / 3) \times 0.5}{(1 / 3) \times 0.5}) = 1

So we accept a proposed move from A to B with a probability of 50%, while a proposed move from B to A is always accepted.

Now our table becomes:

Direction	Starting	Proposal probability	Acceptance probability	Expected volume
A -> B	12	0.5	0.5	3
B -> A	6	0.5	1.0	3

Now the math becomes:

P (A, B) = q (A, B) \times α (A, B) = 0.5 \times 0.5 = 0.25

P (B, A) = q (B, A) \times α (B, A) = 0.5 \times 1 = 0.5

Notice that we now expect 3 moves in each direction. The flow of probability is balanced.

This is called detailed balance. In math:

π (A) P (A, B) = π (B) P (B, A)

For our example:

\frac{2}{3} \times 0.25 = \frac{1}{3} \times 0.5 = \frac{1}{6}

Detailed balance guarantees that our desired distribution is stationary: once the chain has that distribution, another step leaves it unchanged.

What happens when we reject a proposal?

Let's say we're in state A and propose moving to B. Our acceptance probability is 0.5.

We generate a random number between 0 and 1. If it's at most 0.5, we accept the move. Otherwise, we reject it and stay in A.

This is important because rejected proposals also contribute to the transition probabilities of staying in the same state.

Since the chain must either move or stay, we get:

P (A, A) = 1 - P (A, B) = 0.75

P (B, B) = 1 - P (B, A) = 0.5

So why do we even need MH?

We didn't know the transition matrix at first, but by choosing proposal probabilities and calculating the acceptance probabilities, we constructed it:

P = (\begin{matrix} 0.75 & 0.25 \\ 0.50 & 0.50 \end{matrix})

And this matrix has exactly the stationary distribution we wanted:

π P = π

Why not just use distribution directly?

In our example, we absolutely can. We know that A has a probability of 2/3 and B has a probability of 1/3. We could simply sample from those probabilities and skip MH entirely.

So why does MH exist?

Because in a more complex problem, we might not have a complete list of the probabilities of every possible state. Instead, we might have a mathematical expression that lets us compare the probabilities of two states.

And that's exactly what we're already doing in the MH formula:

\frac{π (y)}{π (x)}

For our example:

\frac{π (B)}{π (A)} = \frac{1 / 3}{2 / 3} = 0.5

In a much larger problem, we might be able to calculate this ratio without calculating the probability of every state or the total needed to normalize the distribution.

MH uses this comparison, together with our chosen proposal probabilities, to calculate how often we should accept a proposed move. This lets us generate samples without constructing the entire transition matrix.

After running the chain long enough, the proportion of time it spends in each state approaches our target distribution, provided the chain can properly explore its states. The samples are not independent, since each state depends on the previous one.

So knowing the complete probability distribution and knowing how to compare the probabilities of two states are not the same thing.

That's what makes Metropolis-Hastings useful.