Bayesian Inference: From Bayes' Theorem to MCMC

Bayesian inference is a framework for updating beliefs about the world as new evidence arrives. At its core is one equation — Bayes' theorem — but the implications span everything from A/B testing to large-scale neural network training.

This post walks through the theory step by step: from the theorem itself, through conjugate priors (where the math is clean), and into MCMC (where math meets computation).

1. The Theorem

Bayes' theorem follows directly from the definition of conditional probability:

$P(A \mid B) = \frac{P(B \mid A) \, P(A)}{P(B)}$

In parameter estimation, we rewrite it with $\theta$ (the parameters) and $D$ (the data):

$P(\theta \mid D) = \frac{P(D \mid \theta) \, P(\theta)}{P(D)}$

The denominator expands via the law of total probability:

$P(D) = \int P(D \mid \theta) \, P(\theta) \, d\theta$

This integral is the normalising constant — and it's often the hardest part.

Key insight: The posterior is proportional to likelihood × prior. The denominator is just a scale factor: $P(\theta \mid D) \propto P(D \mid \theta) \, P(\theta)$

Deeper: This proportionality is why MCMC methods work — we can sample from the posterior without ever computing $P(D)$, as long as we can evaluate the numerator up to a constant.

2. Conjugate Priors

A prior is conjugate to a likelihood when the posterior has the same functional form as the prior.

Beta-Binomial

The classic example: binomial likelihood + beta prior = beta posterior.

$y \sim \text{Binomial}(n, \theta)$

$\theta \sim \text{Beta}(\alpha, \beta)$

$\theta \mid y \sim \text{Beta}(\alpha + y, \;\beta + n - y)$

The beta distribution's PDF:

$f(\theta; \alpha, \beta) = \frac{\theta^{\alpha - 1} (1 - \theta)^{\beta - 1}}{B(\alpha, \beta)}$

where $B(\alpha, \beta)$ is the beta function:

$B(\alpha, \beta) = \int_0^1 t^{\alpha - 1} (1 - t)^{\beta - 1} \, dt$

Common conjugate pairs

The Gaussian-normal conjugate posterior demonstrates precision-weighted averaging. The posterior mean:

$\mu_n = \frac{\frac{\mu_0}{\sigma_0^2} + \frac{n \bar{y}}{\sigma^2}}{\frac{1}{\sigma_0^2} + \frac{n}{\sigma^2}}$

This is a precision-weighted combination of the prior mean $\mu_0$ and the data mean $\bar{y}$. As $n \to \infty$, the data dominates and $\mu_n \to \bar{y}$.

3. When Conjugacy Breaks

For most real-world models, no conjugate prior exists and the evidence integral becomes intractable:

$P(D) = \int P(D \mid \theta) \, P(\theta) \, d\theta \qquad \text{(no closed form)}$

This is where we turn to Markov Chain Monte Carlo.

4. Markov Chain Monte Carlo

MCMC constructs a Markov chain whose stationary distribution is the posterior $P(\theta \mid D)$. After enough steps, samples from the chain approximate samples from the posterior.

Metropolis-Hastings

The algorithm:

1. Start at initial $\theta_0$
2. For $t = 1$ to $T$:
   • Propose $\theta^* \sim q(\theta^* \mid \theta_t)$
   • Compute the acceptance ratio: $\alpha = \min\!\left(1,\; \frac{P(\theta^* \mid D) \; q(\theta_t \mid \theta^*)}{P(\theta_t \mid D) \; q(\theta^* \mid \theta_t)}\right)$
   • Accept $\theta^*$ with probability $\alpha$

Intuition: Always move to a more probable state; sometimes move to a less probable state (proportional to the ratio). This prevents the chain from getting stuck in local modes.

A minimal Python implementation:

import numpy as np
 
def metropolis_hastings(
    log_posterior,         # function: log P(theta | D)
    proposal_sampler,      # function: draw theta* ~ q(· | theta)
    proposal_log_prob,     # function: log q(theta* | theta)
    theta_0,               # initial position
    n_steps: int = 10_000,
):
    samples = [theta_0]
    current = theta_0
 
    for _ in range(n_steps):
        proposed = proposal_sampler(current)
 
        # log acceptance ratio
        log_alpha = (
            log_posterior(proposed)
            + proposal_log_prob(current, proposed)
            - log_posterior(current)
            - proposal_log_prob(proposed, current)
        )
 
        if np.log(np.random.uniform()) < log_alpha:
            current = proposed   # accept
 
        samples.append(current)
 
    return np.array(samples)

For a symmetric proposal (e.g. random-walk: $q(\theta^* \mid \theta) = q(\theta \mid \theta^*)$), the ratio simplifies because the proposal terms cancel:

$\alpha = \min\!\left(1,\; \frac{P(\theta^* \mid D)}{P(\theta_t \mid D)}\right)$

Burn-in and thinning

Two practical concerns when sampling:

1. Burn-in — discard the first $B$ samples while the chain converges
2. Thinning — keep every $k$-th sample to reduce autocorrelation

The number of effective samples is:

$n_{\text{eff}} = \frac{T}{1 + 2 \sum_{k=1}^\infty \rho_k}$

where $\rho_k$ is the autocorrelation at lag $k$.

A healthy chain:

Hamiltonian Monte Carlo

For high-dimensional posteriors, Metropolis-Hastings explores too slowly. HMC augments the parameter space with momentum variables and simulates Hamiltonian dynamics to propose distant states with high acceptance.

The Hamiltonian:

$H(\theta, r) = -\log P(\theta \mid D) + \frac{1}{2} r^T M^{-1} r$

where $r$ is the momentum and $M$ is the mass matrix. The leapfrog integrator discretises the dynamics:

def leapfrog(theta, r, grad_log_p, step_size, n_steps):
    """Simulate Hamiltonian dynamics."""
    r += 0.5 * step_size * grad_log_p(theta)
 
    for _ in range(n_steps - 1):
        theta += step_size * r
        r += step_size * grad_log_p(theta)
 
    theta += step_size * r
    r += 0.5 * step_size * grad_log_p(theta)
 
    return theta, r

HMC is the backbone of modern probabilistic programming languages like Stan and PyMC.

5. Bayesian vs Frequentist

A concrete example

A coin flipped 10 times lands heads 8 times. Estimate $\theta$ (probability of heads).

Bayesian with a $\text{Beta}(1, 1)$ (uniform) prior:

$\theta \mid y \sim \text{Beta}(1 + 8,\; 1 + 2) = \text{Beta}(9, 3)$

Posterior mean: $\dfrac{9}{9 + 3} = 0.75$

95% credible interval: $[0.55, 0.89]$

Frequentist maximum likelihood estimate:

$\hat{\theta} = \frac{8}{10} = 0.80$

95% Wald confidence interval:

$0.80 \pm 1.96 \times \sqrt{\frac{0.80 \times 0.20}{10}} = [0.55, 1.05]$

Both give similar numeric answers, but the interpretation is radically different — the Bayesian can directly state "there's a 95% probability $\theta$ lies in $[0.55, 0.89]$", while the frequentist must say "95% of such intervals would contain the true $\theta$".

Wrapping Up

Bayesian inference is beautiful theory — but it's also practical. With MCMC, you can fit models that would be impossible with closed-form statistics.

What to explore next

• Derive Bayes' theorem from the product rule of probability
• Compute a Beta-Binomial posterior by hand
• Implement Metropolis-Hastings for a simple problem
• Read Statistical Rethinking by Richard McElreath
• Explore Hamiltonian Monte Carlo for higher dimensions
• Compare credible vs confidence intervals on real data

Further reading

"All models are wrong, but some are useful." — George Box (1976)