Part 3: Markov Chain Monte Carlo

“It is a capital mistake to theorize before one has data. Insensibly one begins to twist facts to suit theories, instead of theories to suit facts.”
— Arthur Conan Doyle, A Scandal in Bohemia

“The purpose of computing is insight, not numbers.”
— Richard Hamming

Learning Outcomes¶

By the end of Part 3, you will be able to:

• Explain why a posterior distribution is scientifically valuable but computationally intractable to evaluate directly in moderate and high dimensions

• Distinguish optimization from posterior sampling, and explain why finding a peak is not the same as characterizing uncertainty

• Derive the Metropolis-Hastings acceptance rule from detailed balance and interpret what each term is doing

• Implement a Metropolis-Hastings sampler for a target density that can be evaluated up to a normalization constant

• Diagnose whether an MCMC run is trustworthy using trace plots, autocorrelation, effective sample size, and $\hat{R}$

• Connect MCMC back to the course-wide Observable $\to$ Model $\to$ Inference framework by adding the missing computational step

• Apply the full workflow to an astronomical inference problem and interpret the result scientifically

Prerequisites: Parts 1-2 of this module (measurement philosophy, Bayes’ theorem), Module 1 (sampling distributions, CLT, ergodicity), basic Python.

Roadmap: Where We’re Going¶

Part 2 established the epistemic object of Bayesian inference: not a best-fit parameter, but a posterior distribution over plausible parameter values given data and assumptions. Part 3 asks the next scientific question: how do we compute scientifically meaningful quantities from that posterior when direct integration is impossible?

In course language, the logic is now:

• Observable: the measurements we actually make

• Model: the forward relation and noise model that explain how those measurements arise from physical parameters

• Inference: the posterior distribution over unknown parameters

• Computation: the numerical machinery needed to represent that posterior when it cannot be written down or evaluated everywhere

MCMC does not replace inference. It is the computational bridge that makes Bayesian inference scientifically usable.

This lecture should be read through one computational chain:

The roadmap below is the same logic in flow-chart form:

The narrative arc: expectation as the scientific target, integration as the computational obstacle, sampling as the workaround, sampler construction as the algorithmic task, and diagnostics as the argument for trust. The goal is not just to memorize Metropolis-Hastings. The goal is to understand why posterior sampling is the natural computational continuation of Bayesian inference.

This isn’t just a numerical method. It’s a fundamentally different way of thinking about inference that emerged from mid-20th century physics and revolutionized statistics, machine learning, and computational science.

Part 1: The Computational Crisis — Why We Need MCMC¶

The Posterior is a Distribution, Not a Number¶

In Part 2, you learned Bayes’ theorem:

where

• $\theta$ are the model parameters (what we want to infer)

• $D$ is the observed data (what we have)

• $p(D | \theta)$ is the likelihood (how probable the data is given parameters)

• $p(\theta)$ is the prior (our beliefs — what we know about the parameters)

• $p(\theta | D)$ is the posterior (what we want)

• $p(D) = \int p(D|\theta)p(\theta)\,d\theta$ is the evidence (normalization constant)

The posterior $p(\theta \mid D)$ is a probability distribution over the parameter space. It tells you, for every possible value of $\theta$, how probable that value is given the data you observed.

From this point onward, we will use

as shorthand for the target distribution. This is only a change in notation, not a new object.

In the course language:

• Observable: the measured data $D$

• Model: the likelihood and prior you built in Part 2

• Inference: the posterior $\pi(\theta)$

• Computation: the missing question in this lecture is how to work with $\pi(\theta)$ without evaluating it everywhere

Let’s be concrete. Suppose you’re measuring the distance to a Cepheid variable star. The point here is not just to name the pieces, but to line them up in the exact order you would need for implementation.

A Complete Inference Pipeline You Could Reconstruct¶

Take a Cepheid with measured period $P = 10.0~\mathrm{days}$ and observed apparent magnitude $m_{\mathrm{obs}} = 18.50 \pm 0.15~\mathrm{mag}$.

1. Observable

The data are the quantities we actually measured:

• $P = 10.0~\mathrm{days}$

• $m_{\mathrm{obs}} = 18.50~\mathrm{mag}$

• measurement uncertainty $\sigma_m = 0.15~\mathrm{mag}$

2. Forward model

The Leavitt law gives the absolute magnitude:

and the distance modulus gives the predicted apparent magnitude at distance $d$:

This is the physics step: parameters $\to$ predicted observables.

3. Noise model

For this lecture’s 1D toy problem, we assume:

• independent Gaussian measurement noise,

• known standard deviation $\sigma_m = 0.15~\mathrm{mag}$,

• period measured precisely enough that we can hold $P$ fixed.

That assumption becomes

4. Likelihood

Those assumptions imply

5. Posterior

Once we choose a prior $p(d)$, Bayes gives

6. MCMC target

This posterior is exactly the quantity MCMC will target:

The posterior $p(d \mid m_{\mathrm{obs}}, P)$ tells you not just the “best” distance, but the full probability distribution over possible distances. This captures:

• Your uncertainty (how confident are you?)

• Asymmetries (is the uncertainty symmetric around the peak?)

• Correlations (if you had multiple parameters, how do they trade off?)

In one dimension, you could still evaluate $p(d \mid m_{\mathrm{obs}}, P)$ on a grid of distance values and plot it. That is useful pedagogically because it is the last case where brute force still feels easy.

But the moment the physics becomes slightly more realistic, the problem stops being one-dimensional.

If extinction $A_V$ is also unknown, then the forward model becomes

Now the same observed brightness can be explained by a closer, dustier Cepheid or a farther, cleaner one. That is a genuine physical degeneracy, not a numerical inconvenience.

{width=“100%” #fig-mcmc-degeneracy-2d}

Even this toy extension becomes expensive quickly. A grid with 400 points in $d$ and 400 points in $A_V$ already needs $1.6 \times 10^5$ likelihood evaluations. Add metallicity, intrinsic scatter, a calibration offset, and one nuisance parameter for each observed star, and brute force becomes the bottleneck.

The Curse of Dimensionality¶

Imagine you want to map out a posterior distribution by evaluating it on a grid. How many grid points do you need?

1D example: To resolve a 1D posterior to reasonable accuracy, you might need $N = 100$ grid points.

2D example: Now you have two parameters. You need $100 \times 100 = 10^4$ evaluations.

3D example: Three parameters $\to$ $100 \times 100 \times 100 = 10^6$ evaluations.

Scaling: For $d$ dimensions, you need $N^d$ evaluations. This is exponential scaling — the curse of dimensionality.

Let’s make this concrete with realistic numbers:

Even “small” problems with 5-10 parameters become computationally impossible. And many real scientific problems have hundreds or thousands of parameters:

• Climate models: ~100 parameters (ocean circulation, cloud physics, aerosols...)

• Exoplanet characterization: ~20 parameters per planet (mass, orbit, atmosphere...)

• Galaxy formation: ~50 parameters (star formation efficiency, black hole feedback...)

• Neural networks: millions to billions of parameters

Grid evaluation simply cannot work.

What About Optimization? Just Find the Peak¶

You might think: “Why do I need the full distribution? Just find the maximum of $p(\theta \mid D)$ and report that.”

That reasoning is appealing because optimization is familiar and computationally simpler. But sampling and optimization answer different epistemic questions.

• Optimization asks: Which parameter value maximizes the posterior?

• Posterior inference asks: Which parameter values remain plausible, and with what relative credibility, given the data and model assumptions?

The second question is scientifically richer because uncertainty, asymmetry, and parameter tradeoffs are often part of the result, not peripheral details.

This is maximum a posteriori (MAP) estimation, and it’s sometimes useful. But you lose crucial information:

1. Uncertainty quantification: The peak tells you nothing about how uncertain you are. Is it a sharp peak (confident) or broad (uncertain)? If you report $\theta = 5.0$ as the MAP estimate, can you distinguish whether your uncertainty is $\pm 0.1$ or $\pm 10$? The peak doesn’t tell you.

2. Asymmetries: The peak might be at $d = 500$ pc, but the distribution could be skewed — maybe $500 \pm 20$ pc toward smaller distances, but +100 pc toward larger distances (imagine a long tail from dust extinction uncertainty).

3. Correlations: With multiple parameters, knowing the peak doesn’t tell you about degeneracies. Maybe radius and temperature are correlated — high R and low T give the same luminosity as low R and high T.

4. Model comparison: The evidence $p(D)$ (the denominator in Bayes’ theorem) requires integrating over all parameter values. The peak tells you nothing about this.

5. Propagating uncertainty: If you use your inferred parameters to predict something else, you need the full distribution to get realistic error bars.

The scientific standard: Modern quantitative science does not stop at “best fit” values. It reports credible intervals and posterior structure. Bayesian posteriors do this directly.

{width=“100%” #fig-mcmc-posterior-full-object}

What Inference Actually Requires¶

Now we can state the computational target explicitly.

Everything in the rest of the lecture is in service of the same chain:

In science, we rarely want the posterior curve itself. We want derived quantities:

• the posterior mean distance,

• the probability that a parameter lies in some physically interesting range,

• a credible interval,

• a prediction under the forward model.

All of these can be written as posterior expectations or posterior integrals. For any function $f(\theta)$,

is the general quantity we want to compute.

For the Cepheid problem, examples include:

• $f(d) = d$ for the posterior mean distance,

• $f(d) = m_{\mathrm{model}}(d)$ for the mean predicted magnitude,

• indicator functions such as $f(d) = \mathbf{1}_{d < 500~\mathrm{kpc}}$ for probabilities of physically relevant events.

This is the conceptual hinge of the lecture:

• Inference gives us the posterior $\pi(\theta)$.

• Computation must turn that posterior into expectations we can report scientifically.

At this point, all of the standard strategies have failed:

• grid evaluation is intractable,

• analytic integration is impossible,

• and independent sampling is unavailable.

So MCMC is not merely a convenient tool we happen to like. Once direct integration, dense grids, and direct iid (independent and identically distributed) sampling are off the table, constructing a sampler becomes the natural general-purpose strategy if we still want posterior expectations. The natural general-purpose option is to construct a process whose long-run behavior reproduces the posterior. That conceptual move is what defines MCMC.

The Monte Carlo Solution: Sample, Don’t Integrate¶

Here’s where the inference becomes computational. If we had iid samples $\theta_1, \theta_2, \ldots, \theta_N$ drawn from the posterior, then the intractable integral above would become a sample average:

This is the Monte Carlo idea: replace integration with averaging over samples. It works because an expectation under $\pi(\theta)$ is, by definition, an average of $f(\theta)$ weighted by that distribution; Monte Carlo replaces the exact weighted integral with an empirical average over draws from the same distribution.

If you have $N$ iid samples $\{\theta_1, \theta_2, \ldots, \theta_N\}$ drawn from $p(\theta \mid D)$, then you can estimate many posterior quantities:

Mean:

Variance:

Credible intervals: Sort the samples and find the 16th and 84th percentiles for a 68% credible interval.

Predictions: For any function $f(\theta)$, estimate its expectation:

The key quantitative fact: For the Monte Carlo estimator

with $\theta_i \stackrel{\text{iid}}{\sim} \pi(\theta)$ and $\sigma_f^2 = \mathrm{Var}[f(\theta)]$, the estimator variance is

so the standard error — the standard deviation of the estimator itself — is

The distinction matters:

• $\mathrm{Var}(\hat{\mu}_f)$ tells you how the estimator’s variance shrinks with sample size,

• $\mathrm{SE}(\hat{\mu}_f)$ is the uncertainty scale for the estimator itself.

The Central Limit Theorem implies that the sampling distribution of the estimator has variance $\sigma_f^2/N$, so its standard deviation (the standard error) scales as $\sigma_f/\sqrt{N}$.

To gain one more decimal place of accuracy, you need about $100 \times$ more samples. That is slow convergence. But it is still fundamentally better than dense-grid integration in high dimensions, where the cost typically grows exponentially with the number of parameters.

This is revolutionary. Grid evaluation scales exponentially with dimension $(N^d)$. By contrast, iid Monte Carlo has an $O(N^{-1/2})$ convergence rate that does not itself worsen exponentially with dimension.

But here’s the catch: We do not have a machine that hands us independent samples from $\pi(\theta) = p(\theta \mid D)$. We can evaluate the posterior at points we choose, but evaluation is not the same as sampling.

So the problem changes form:

How do we generate samples from a distribution we can evaluate but cannot sample from directly?

This is where Markov Chain Monte Carlo enters. It’s a method to generate samples from a very broad class of target distributions by constructing a Markov chain that explores parameter space.

Part 2: Markov Chains — The Mathematical Foundation¶

The Core Idea: Build a Stochastic Process¶

Because direct independent sampling from $\pi(\theta) = p(\theta \mid D)$ is unavailable, we replace it with a stochastic process that produces dependent samples whose long-run frequencies still match the posterior.

This is the key conceptual relaxation. If independent samples are impossible, then we must allow dependence between successive samples. The simplest structure that introduces dependence without becoming hopelessly complicated is a Markov process: the next state depends only on the current one.

That is why Markov chains appear here. They are not arbitrary decoration. They are the minimal relaxation of independence that still gives us a tractable way to construct a sampler.

Instead of trying to sample directly from $p(\theta \mid D)$, we’ll build a stochastic process (a Markov chain) that:

1. Starts from any initial $\theta_0$

2. Takes random steps through parameter space

3. Gradually “forgets” where it started

4. Eventually produces samples distributed according to $p(\theta \mid D)$

Think of it like a random walker wandering around parameter space. We’ll design the rules of walking so that after sufficient wandering, the walker spends time in each region proportional to $p(\theta \mid D)$.

If $p(\theta \mid D)$ is high in some region, the walker visits often. If $p(\theta \mid D)$ is low, the walker rarely goes there. The stationary distribution of the random walk is the posterior.

Markov Chains: Formal Definition¶

A Markov chain is not yet an inference method. It is a type of stochastic process. We choose it because it is the simplest sequential mechanism rich enough to produce dependent samples with a controllable long-run distribution.

Its defining property is simple:

the next state depends only on the current state, not on the full past history.

In symbols,

This is the Markov property: the future depends on the present, not the past. Once that idea is clear, we can introduce the next crucial object: the transition kernel $T(\theta' \mid \theta)$, which tells us how the chain moves through parameter space.

Transition Kernels¶

The dynamics of a Markov chain are defined by the transition kernel (or transition probability):

This is the probability density of moving to $\theta'$ given that you’re currently at $\theta$. It must satisfy:

for all $\theta$ (probabilities must sum to 1).

Given a starting distribution $p_0(\theta)$, the distribution after one step is:

After $t$ steps:

Question: Can we design $T$ such that $p_t(\theta)$ converges to our target distribution $\pi(\theta) = p(\theta \mid D)$ as $t \to \infty$?

Answer: Yes! But we need two conditions:

1. Stationarity: The target distribution $\pi$ must be stationary with respect to $T$.

2. Ergodicity: The chain must be able to explore the entire parameter space.

Stationarity: The Distribution Must Be Preserved¶

A distribution $\pi(\theta)$ is stationary with respect to $T$ if applying the transition kernel doesn’t change it:

What it is: If the Markov chain is distributed according to $\pi$, then after one more step, it’s still distributed according to $\pi$. The distribution is an equilibrium of the dynamics.

Physical analogy: A gas in thermal equilibrium. Molecules are constantly moving and colliding, but the overall distribution (Boltzmann distribution) doesn’t change. Individual states change, but the probability distribution is stationary.

These two ideas play different roles:

• Stationarity tells us the target distribution is preserved by the dynamics.

• Ergodicity tells us the chain can actually find and explore that distribution from arbitrary starting points.

A chain can have the correct stationary distribution on paper and still be scientifically useless if it cannot explore the relevant regions of parameter space.

Ergodicity: The Chain Must Explore Everywhere¶

Stationarity means $\pi$ is preserved if you start from it. But how do we get there from an arbitrary starting point $\theta_0$?

We need ergodicity: The chain must be able to reach any state from any other state in finite time. Formally, ergodicity requires:

1. Irreducibility: Every state is accessible from every other state (there are no isolated regions).

2. Aperiodicity: The chain doesn’t get stuck in cycles. (E.g., don’t alternate deterministically: left, right, left, right forever.)

The Ergodic Theorem: If a Markov chain is ergodic and has stationary distribution $\pi$, then for any function $f(\theta)$:

Time averages along the chain converge to expectations under $\pi$. This is how MCMC works: run the chain for a long time, then use the samples to estimate expectations.

Detailed Balance: A Sufficient Condition for Stationarity¶

At this point we need a design test. A chain that merely “moves around” parameter space is not enough. A random walk could easily settle into the wrong distribution. We need a condition we can verify mathematically that guarantees the chain targets the right long-run distribution.

Here is the failure mode to keep in mind. Suppose we design a random walk that proposes steps and always accepts them. It will certainly explore parameter space, but it will not generally spend time in each region proportional to $\pi(\theta)$. In fact, for a symmetric proposal on a bounded space, it tends to behave much more like uniform exploration than posterior sampling.

So “motion” is not enough. We need a correction mechanism that reshapes the dynamics so the chain has the right equilibrium distribution. Detailed balance is the first clean way to express that requirement mathematically.

Detailed balance is not the only possible route to stationarity, but it is a powerful sufficient condition because it is local, checkable, and easy to build into an algorithm.

We need to design a transition kernel $T$ that has $\pi$ as its stationary distribution. There’s a powerful sufficient condition called detailed balance (also called microscopic reversibility):

In words: the probability flow from $\theta$ to $\theta'$ equals the flow from $\theta'$ back to $\theta$ when both are drawn from $\pi$.

Why does detailed balance imply stationarity? Here’s the key insight:

Detailed balance is a sufficient condition for stationarity. To see why, assume detailed balance holds and integrate both sides over $\theta$:

The right side factors since $\pi(\theta')$ doesn’t depend on the integration variable:

(The integral equals 1 because T is a probability density in its first argument.)

The left side is exactly the definition of how the distribution evolves forward one step:

This integral gives you the probability at $\theta'$ after one transition starting from $\pi$. Since $\text{LHS} = \text{RHS} = \pi(\theta')$, we’ve shown that starting from distribution $\pi$ produces distribution $\pi$ after one step. The distribution is stationary.

Now we have our blueprint:

1. Choose a target distribution $\pi(\theta) = p(\theta \mid D)$.

2. Design a transition kernel $T$ that satisfies detailed balance with respect to $\pi$.

3. The Markov chain will converge to $\pi$ (if ergodic).

4. Run the chain, collect samples, estimate expectations

The question: How do we actually construct such a $T$? That’s where Metropolis-Hastings comes in.

Part 3: The Metropolis-Hastings Algorithm¶

The Construction: Proposal + Accept/Reject¶

Here’s the brilliant idea due to Metropolis et al. (1953) and generalized by Hastings (1970):

Separate the transition kernel into two parts:

1. Proposal distribution $Q(\theta' \mid \theta)$: generates a candidate new state $\theta'$ given current state $\theta$

2. Acceptance probability $\alpha(\theta' \mid \theta)$: probability of accepting the proposal

For $\theta' \neq \theta$, the accepted-move component of the transition kernel is:

The full kernel also includes a self-transition at $\theta' = \theta$ that collects the rejected proposals. In other words, “reject and stay where you are” is part of the Markov dynamics itself, not an implementation afterthought.

The algorithm:

1. Start at some $\theta_0$ (anywhere in parameter space)

2. For $t = 0, 1, 2, ..., N-1$:

   • Propose: Draw $\theta'$ from $Q(\theta' \mid \theta_t)$

   • Evaluate: Compute acceptance probability $\alpha(\theta' \mid \theta_t)$

   • Decide:

     • With probability $\alpha$, accept: $\theta_{t+1} = \theta'$

     • With probability $1-\alpha$, reject: $\theta_{t+1} = \theta_t$ (stay put)

3. Return samples $\{\theta_1, \theta_2, \ldots, \theta_N\}$

The magic is in step 2: How do we choose $\alpha$ to satisfy detailed balance?

Deriving the Acceptance Probability¶

We want:

Substituting $T = Q\alpha$ where $Q$ is the proposal and $\alpha$ is the acceptance probability:

Rearrange:

Define the acceptance ratio:

When $\pi$ is the posterior: write

Then

so the evidence $p(D)$ cancels because it is constant in $\theta$.

We need $\alpha(\theta' \mid \theta)$ and $\alpha(\theta \mid \theta')$ to have ratio $r$. A simple choice that works:

Verification: Let

• If $R \geq 1$: $\alpha(\theta' \mid \theta) = 1$, $\alpha(\theta \mid \theta') = 1/R \leq 1$. Ratio: $1/(1/R) = R$.

• If $R < 1$: $\alpha(\theta' \mid \theta) = R$, $\alpha(\theta \mid \theta') = 1$. Ratio: $R/1 = R$.

Either way, the ratio condition is satisfied, so detailed balance holds!

Interpretation:

• If $\theta'$ has higher posterior probability than $\theta$ $(R > 1)$, always accept.

• If $\theta'$ has lower posterior probability $(R < 1)$, accept with probability $R$.

The accept/reject rule does something subtle and important:

• it favors moves toward higher posterior density,

• but it does not force greedy ascent.

That is exactly why MCMC samples a distribution instead of performing optimization. Occasional moves into lower-probability regions are not mistakes; they are part of correctly representing posterior uncertainty.

Proposal Distributions: Symmetric vs. Asymmetric¶

The acceptance probability depends on the proposal distribution $Q$. Two important cases:

Symmetric proposals: $Q(\theta' \mid \theta) = Q(\theta \mid \theta')$.

Example — Gaussian random walk (symmetric in $\theta$):

For symmetric $Q$, the $Q$-terms cancel and

Here you can evaluate $\pi$ using only log-likelihood + log-prior; any constant normalizer cancels. This is the Metropolis algorithm (the original 1953 version). You only need to evaluate the ratio of posterior probabilities.

Asymmetric proposals: $Q(\theta' \mid \theta) \neq Q(\theta \mid \theta')$

Example: Propose $\theta'$ by drawing from the prior $p(\theta)$.

Then:

You need the full ratio including the $Q$ terms. This is the Metropolis-Hastings algorithm (Hastings 1970 generalization).

Tuning the Proposal: The Art of Step Size¶

The proposal distribution Q determines how the chain explores. A crucial parameter is the step size (scale of proposals).

Too small: Proposals are always accepted $(\alpha \approx 1)$, but you take tiny steps. The chain explores very slowly. High autocorrelation between successive samples.

Too large: Proposals often land in low-probability regions. High rejection rate $(\alpha \approx 0)$. The chain stays stuck in one place most of the time. Again, slow exploration.

Just right: Moderate acceptance rate (typically 20-50%). The chain explores efficiently, taking substantial steps while still accepting frequently enough.

Adaptive tuning: In practice, you might run the chain for a burn-in period, monitor the acceptance rate, and adjust $\sigma$. Common strategy:

• Run 1000 steps (burn-in)

• If acceptance rate > 50%, multiply $\sigma$ by 1.2

• If acceptance rate < 20%, divide $\sigma$ by 1.2

• Repeat until acceptance rate is in target range

Once tuned, fix the proposal and run your production chain. (Don’t keep adapting during production — this violates the Markov property and detailed balance.)

From Likelihood to Log-Posterior Code¶

This is the bridge students usually need made explicit. For any Bayesian model, the implementation pipeline is:

1. Write the forward model that predicts observables from parameters.

2. Write the log-likelihood from your noise model.

3. Write the log-prior from your scientific assumptions and parameter bounds.

4. Add them to get the unnormalized log-posterior:

where the constant $C = -\log p(D)$ can be ignored.

In practice, this is the computational object nearly everything reduces to:

def log_posterior(theta):
    return log_likelihood(theta) + log_prior(theta)

MCMC algorithms work by comparing this quantity between the current and proposed states. Ratios in probability space become subtraction in log space, which is why log-space arithmetic is not optional numerical polish but part of the method itself.

5. Run the sampler on differences in $\log \pi(\theta)$ between the current and proposed states.

That is the whole computational bridge:

Purpose: the code below is a minimal reusable Metropolis-Hastings skeleton.

Inputs: a function that evaluates the unnormalized log-posterior, a proposal rule, an initial parameter vector, and the number of samples.

Outputs: a Markov chain of sampled parameter states plus the acceptance rate.

Common pitfall: on rejection, you must still record the current state again. Rejected proposals do not disappear from the chain; “staying put” is part of the Markov process.

Minimal Sampler Skeleton¶

import numpy as np  # for np.log, np.random.uniform, array handling

def metropolis_hastings(log_posterior, proposal, theta_init, n_samples):
    """
    Generic Metropolis-Hastings MCMC sampler.
    
    Parameters:
    -----------
    log_posterior : function
        Computes log pi(theta) (log posterior probability). Return the unnormalized log posterior: log p(D|theta) + log p(theta). Do not include log p(D); it is a constant and cancels in MH.
    proposal : function  
        Generates theta_prime given theta. Returns theta_prime and log_Q_ratio = log Q(theta|theta_prime) - log Q(theta_prime|theta) (order matters)
    theta_init : array
        Initial parameter values
    n_samples : int
        Number of MCMC samples to generate
    
    Returns:
    --------
    chain : array (n_samples, n_params)
        MCMC samples
    acceptance_rate : float
        Fraction of accepted proposals
    """
    
    theta = theta_init.copy()
    chain = np.zeros((n_samples, len(theta)))
    n_accepted = 0
    
    log_pi_current = log_posterior(theta)
    
    for i in range(n_samples):
        # 1. Propose new state
        theta_proposed, log_Q_ratio = proposal(theta)
        
        # 2. Evaluate posterior at proposed state
        log_pi_proposed = log_posterior(theta_proposed)
        
        # 3. Compute acceptance probability (in log space!)
        log_alpha = log_pi_proposed - log_pi_current + log_Q_ratio
        
        # 4. Accept or reject
        if np.log(np.random.uniform()) < log_alpha:
            # Accept
            theta = theta_proposed
            log_pi_current = log_pi_proposed
            n_accepted += 1
        # If rejected, theta stays the same (automatic in this code structure)
        
        # 5. Store sample
        chain[i] = theta
    
    acceptance_rate = n_accepted / n_samples
    return chain, acceptance_rate

Key implementation details:

1. Work in log space

   Purpose: show why a sampler evaluates log-probabilities instead of raw probabilities. Inputs: two extremely small probability values. Output: a finite log-ratio that can still be compared numerically. Pitfall: direct probabilities underflow to zero, so ratio calculations become undefined.

# Direct space: numerical disaster
p1 = 1e-1000  # Underflows to 0!
p2 = 1e-1001  # Also underflows to 0!
ratio = p2 / p1  # 0/0 = NaN. Game over.

# Log space: works perfectly  
log_p1 = -1000 * np.log(10)  # ~-2302.6
log_p2 = -1001 * np.log(10)  # ~-2305.0  
log_ratio = log_p2 - log_p1  # -2.3 (a perfectly fine number!)
# If log_p = -1000, then p = exp(-1000) = 0 (underflow!)
# But log_p = -1000 is a perfectly fine number to work with.

2. Proposal functions

   Purpose: generate a candidate move for a symmetric random-walk sampler. Input: current parameter vector and a proposal scale. Output: proposed parameter vector plus the proposal-density correction term. Pitfall: if the proposal is not actually symmetric, you cannot set log_Q_ratio = 0.0.

def gaussian_proposal(theta, sigma=1.0):
    theta_new = theta + np.random.normal(0, sigma, size=len(theta))
    log_Q_ratio = 0.0  # Symmetric, so Q(theta|theta_prime)/Q(theta_prime|theta) = 1
    return theta_new, log_Q_ratio

3. Acceptance decision: Compare log(uniform random) to $\log \alpha$ instead of comparing a uniform random draw to $\alpha$. This avoids exponentiating large negative numbers. (Use unnormalized log posterior: loglike + logprior; do not add log p(D).)

What About the Normalization Constant?¶

Recall: we only need $\pi(\theta)$ up to proportionality; this section shows explicitly how $p(D)$ cancels from the MH acceptance ratio.

Remember Bayes’ theorem:

The denominator $p(D)$ is often intractable (that’s why we’re doing MCMC in the first place!). But look at the acceptance probability:

The $p(D)$ terms cancel! You only need:

Profound implication: MCMC only requires evaluating the unnormalized posterior (likelihood $\times$ prior). You never need the evidence $p(D)$.

This is why MCMC works for arbitrarily complex models. As long as you can evaluate $p(D \mid \theta)$ (your likelihood/model) and $p(\theta)$ (your prior), you can sample from the posterior.

Part 4: Practical Diagnostics — Has Your Chain Converged?¶

You’ve run your MCMC algorithm for N steps. You have N samples. Now the crucial question: Are these samples from the posterior?

The ergodic theorem says the chain eventually converges to $\pi$. But how long is “eventually”? If you haven’t run long enough, your samples don’t represent the posterior — they’re still biased toward your starting point.

This is the burn-in problem: Early samples are not from $\pi$. You need to discard them.

Convergence is therefore not something you “prove” from a single plot. It is a judgment built from several pieces of evidence that all point in the same direction: stationarity, mixing, agreement across chains, and enough effectively independent samples.

Visual Diagnostic: Trace Plots¶

The first and most important diagnostic: Plot your samples over time.

{width=“100%” #fig-mcmc-trace-plot-comparison}

A trace plot shows $\theta_t$ vs. $t$. What to look for:

Good signs:

• Looks like random noise around a stable mean

• No trends (increasing or decreasing)

• No periodic oscillations

• Mixes well (explores the full range rapidly)

• Multiple chains from different starting points overlap

Bad signs:

• Systematic trend (chain hasn’t equilibrated)

• Stuck in one region for a long time (mixing poorly)

• Periodic oscillations (not ergodic, perhaps?)

• Multiple chains give different distributions (not converged)

Quantitative Diagnostic 1: Autocorrelation Function¶

The autocorrelation function (ACF) measures how correlated $\theta_t$ is with $\theta_{t+k}$:

where $\mu = \langle \theta \rangle$ is the chain mean.

Interpretation:

• $\rho_0 = 1$ (perfect correlation with itself)

• $\rho_k \to 0$ as $k \to \infty$ for an ergodic chain

• Larger $k$ needed for $\rho_k \approx 0$ means higher correlation (worse mixing)

Autocorrelation time $\tau$: the characteristic lag at which samples become independent:

In practice, truncate the infinite sum when $\rho_k$ becomes negligible (say, when $\rho_k < 0.05$ or when $\rho_k$ is no longer statistically significant).

Effective sample size (ESS):

If $\tau = 10$, then your $N = 10000$ samples are only as informative as $N_\text{eff} = 1000$ independent samples. Your error bars should be computed using $N_\text{eff}$, not $N$.

Quantitative Diagnostic 2: Gelman-Rubin Statistic¶

The Gelman-Rubin diagnostic (also called R-hat, $\hat{R}$) tests convergence by comparing multiple chains.

Idea: Run M chains from different starting points. If they’ve all converged to the same distribution, they should have:

• Similar means

• Similar variances

• Samples that are indistinguishable when pooled

The R-hat statistic compares variance between chains to variance within chains:

More precisely:

1. Run $M$ chains of length $N$ each (after burn-in)

2. Compute within-chain variance W (average variance across chains)

3. Compute between-chain variance B (variance of chain means)

4. Estimate total variance: $\hat{\text{Var}} = \frac{N-1}{N} W + \frac{1}{N} B$

5. Compute: $\hat{R} = \sqrt{\hat{\text{Var}} / W}$

Interpretation:

• $\hat{R} \approx 1$: Chains have converged (within-chain and between-chain variances agree)

• $\hat{R} > 1.1$: Chains have not converged (they give different answers)

• $\hat{R} > 1.2$: Serious convergence problem

Rule of thumb: $\hat{R} < 1.01$ for all parameters before trusting your results.

Burn-In: How Much to Discard?¶

The first part of each chain is burn-in: Samples before the chain has reached equilibrium. You must discard these before computing statistics.

How long is burn-in?

No universal answer. It depends on:

• How far your starting point is from the typical set

• How well your proposals are tuned

• The geometry of your posterior (simple unimodal vs. complicated multimodal)

• Dimensionality

Conservative approach:

1. Run the chain for $N$ total steps

2. Plot trace plots for several parameters

3. Identify where the chain stabilizes (stops drifting)

4. Discard everything before that point (typically 10-50% of samples)

5. Check that R-hat < 1.01 for the remaining samples

Practical tip: It’s safer to discard too much than too little. If you have 100,000 samples and discard the first 20,000, you still have plenty. But if you keep burn-in samples, your posterior estimates will be wrong.

Putting It All Together: A Convergence Checklist¶

Before trusting your MCMC results, verify:

Visual inspection:

• Trace plots look like stationary noise (no trends or drift)

• Trace plots from multiple chains overlap

• Posterior distributions from each chain look similar

Quantitative metrics:

• R-hat < 1.01 for all parameters

• Effective sample size ESS $> 400$ per parameter (for 95% credible intervals)

• Autocorrelation time $\tau < 100$ (or you’ve run long enough that $N/\tau > 400$)

Acceptance rate (for Metropolis-Hastings):

• Acceptance rate between 20% and 50%

• Not too high (>60% means proposals too timid)

• Not too low (<15% means proposals too aggressive)

Sensitivity checks:

• Results don’t change if you run longer

• Results don’t change if you discard more burn-in

• Results don’t change if you change starting points

If all these pass, you can trust that your samples approximate the true posterior.

Failure Modes to Watch For¶

Not every chain is scientifically useful, even if it produces a long list of numbers.

• Poor mixing: the chain moves so slowly that neighboring samples carry nearly the same information.

• Non-ergodic behavior: the chain cannot reach some relevant parts of parameter space, so time averages never represent the full posterior.

• Multimodal trapping: different chains stay in different modes and never communicate.

• Bad proposal scaling: proposals are so timid that exploration is glacial, or so aggressive that almost everything is rejected.

These are not cosmetic issues. Each one breaks the argument that sample averages represent posterior expectations.

Part 5: Real Example — Cepheid Variable Distance¶

Let’s apply everything you’ve learned to a concrete astrophysical problem: Inferring the distance to a Cepheid variable star.

The Science: Cepheids as Standard Candles¶

Cepheid variables are pulsating stars discovered by Henrietta Leavitt in 1908. She found a remarkable relationship: The period of pulsation correlates with intrinsic luminosity.

This makes Cepheids standard candles: If you measure the period $(P)$, you know the absolute magnitude (M). Compare to the observed apparent magnitude $(m)$, and you can infer the distance.

Notice the full chain here:

• We do not observe distance directly. We observe apparent magnitude.

• We introduce a model connecting distance to expected magnitude through the distance modulus.

• We encode measurement uncertainty with a likelihood and prior assumptions with a prior.

• The result is an inference: a posterior distribution over distance.

MCMC enters only after this reasoning is complete. Its job is not to create the science question, but to make the answer computable.

The distance modulus relation:

where d is the distance in parsecs.

The Period-Luminosity relation (Leavitt Law):

For Classical Cepheids in the V-band, the absolute magnitude depends on the logarithm of the period:

This gives $M_V$ for a reference metallicity. In practice, there are corrections for:

• Metallicity [Fe/H]

• Reddening/extinction

• Period-color-luminosity relations

For our example, we’ll use the idealized case: A Cepheid with well-measured period, negligible extinction, and solar metallicity.

The Problem: Inferring Distance from Photometry¶

What you observe:

• Period: P = 10.0 days (measured from light curve)

• Apparent magnitude: $m_V = 18.50 \pm 0.15$ mag (from photometry)

What you want to infer:

• Distance: $d$ in parsecs (or kiloparsecs)

• Uncertainty: Full posterior $p(d \mid m, P)$

Step 1: Forward Model¶

The physics runs forward from period and distance to apparent magnitude.

From the Leavitt Law with $P = 10.0$ days:

The distance modulus is:

So:

Solving for $d$:

That gives the familiar point estimate, but it does not yet tell us how uncertain we are.

Step 2: Noise Model and Likelihood Assumptions¶

To turn this into inference, we need to say how noise enters the measurement. For the basic 1D example, we assume:

• the period is known well enough to hold fixed,

• the magnitude uncertainty is Gaussian,

• the standard deviation is known: $\sigma_m = 0.15~\mathrm{mag}$,

• different error contributions are independent enough to combine into one Gaussian width.

The Gaussian choice is not arbitrary. It is a first-pass model justified by the Central Limit Theorem: photometric uncertainty is often the sum of many small effects such as photon-counting noise, detector read noise, and calibration scatter. If those assumptions failed badly, we would need a different likelihood.

Holding $P$ fixed also means that, within this simplified example,

because all period uncertainty has been neglected. In a more realistic model, uncertainty in $P$ would be part of the likelihood too.

That gives the observation model

Step 3: Likelihood¶

The corresponding likelihood is

with

• $m_{\text{model}}(d) = M_V + 5 \log_{10}(d \times 1000 / 10)$ (d in kpc, distance in pc for the log term)

• $\sigma_m = 0.15$ mag (measurement uncertainty)

or, in the exponential form students often implement first,

In log form, which is what the code should evaluate,

{width=“100%” #fig-mcmc-likelihood-1d}

Step 4: Prior¶

We now encode a simple scale prior. We will use a prior uniform in log-distance, which corresponds to the density

This is a common choice for positive scale parameters because it does not privilege one multiplicative scale over another.

We will actually sample in $\log d$, where the prior is uniform:

This corresponds to distances between 50 kpc and 10 Mpc, a reasonable range for extragalactic Cepheids (e.g., in M31, M33, or Local Group galaxies).

Step 5: Posterior¶

By Bayes’ theorem,

In log form, the quantity we actually hand to the sampler is

where the additive constant $C$ does not depend on $d$ and can be dropped.

This is the answer to the practical question “What function do I write?” You write the function that returns $\log p(m_{\mathrm{obs}} \mid d) + \log p(d)$.

Step 6: What Changes in More Than One Parameter?¶

Real Cepheid inference is rarely this clean. If extinction is also uncertain, then the model becomes

Now distance and extinction trade off against each other: a larger $A_V$ can mimic a larger distance. The posterior becomes a correlated ridge rather than a single narrow bump. That is exactly the kind of geometry MCMC is built for.

In that 2D version, the computational target becomes

This is the real conceptual bridge from Part 2 to Part 3: the physics and uncertainty model define the posterior, and MCMC supplies the machinery for exploring it when that posterior is no longer something you can exhaustively map on a grid.

Pause and Predict¶

Before seeing any numerical result, make a scientific prediction:

1. Roughly what distance should the posterior center on?

2. Should the posterior be nearly symmetric in $d$, or slightly skewed?

3. If the proposal scale is tuned reasonably, do you expect the chain to mix quickly in this one-parameter problem?

This matters because MCMC should not feel like a black box. You should have physical expectations before the chain runs.

What the Sampler Should Recover¶

If you connect the generic Metropolis-Hastings sampler above to this log-likelihood and prior, the chain should recover a posterior centered near $d \approx 323~\text{kpc}$ with a width of a few tens of kiloparsecs.

The important teaching point is not the exact number. It is the structure of the answer:

• a central distance estimate near the distance-modulus calculation

• a credible interval rather than a single “best” value

• a trace plot that reaches stationarity quickly in this simple one-dimensional problem

• an autocorrelation time small enough that a modest chain produces a useful effective sample size

For a well-tuned one-dimensional sampler here, an acceptance rate in the rough range of 25% to 50% is plausible, and an autocorrelation time of order 10 would indicate healthy mixing.

Scientific Interpretation¶

What changes when we switch from a point estimate to a posterior?

• We no longer say only “the star is about 323 kpc away.”

• We say “the data support a range of distances, with quantified uncertainty.”

• We can propagate that uncertainty forward into later inferences.

• We can compare this result consistently with other distance indicators.

That is the real payoff of MCMC. It turns an inverse problem into a defensible uncertainty statement.

Synthesis — The Conceptual Architecture¶

You’ve learned a lot. Let’s step back and see the structure of what you now know.

The Four Pillars of MCMC¶

Each pillar rests on the previous one. You can’t appreciate MCMC diagnostics without understanding Markov chains. You can’t design good samplers without knowing detailed balance. And you can’t trust your results without proper statistical foundations.

One way to summarize the entire lecture is:

• Part 1: why posterior computation is hard

• Part 2: what kind of stochastic process can solve it

• Part 3: how Metropolis-Hastings builds that process

• Part 4: how we decide whether to trust the run

• Part 5: what scientific interpretation looks like once the machinery works

Universal Patterns Across Modules¶

One of the profound lessons of this course: The same mathematical structures appear everywhere. Let’s make this explicit:

The profound unity: Whether you’re:

• Drawing random samples to estimate a population mean (Module 1)

• Simulating molecules reaching thermal equilibrium (Module 3)

• Sampling parameters to quantify uncertainty (Module 5)

You’re using the same mathematical structure: stochastic processes converging to equilibrium distributions, with time averages equaling ensemble averages.

This is why physics is so powerful and predictive. The same principles work everywhere, from atoms to stars to probability distributions.

From Theory to Practice: Your Toolbox¶

You now have:

1. Conceptual understanding:

• Why MCMC is necessary (curse of dimensionality)

• How MCMC works (Markov chains, detailed balance)

• When MCMC converges (ergodicity, stationarity)

2. Practical skills:

• Implementing Metropolis-Hastings

• Tuning proposal distributions

• Diagnosing convergence (trace plots, R-hat, ACF)

• Interpreting results (corner plots, credible intervals)

3. Scientific applications:

• Setting up Bayesian models (likelihood + prior)

• Running inference on real data (Cepheids)

• Quantifying uncertainty properly

4. Connections:

• To Module 1 (sampling, CLT, ergodicity)

• To Module 3 (phase space, thermalization, detailed balance)

• To Project 4 (Type Ia SNe, cosmological parameters)

This is a complete inference framework you can apply to any scientific problem.

Common Misconceptions About MCMC¶

Before moving to Project 4, let’s address some frequent misunderstandings that trip up even experienced practitioners:

Misconception 1: “High acceptance rate means good sampling”

Reality: High acceptance (>60%) usually means your proposals are too timid — you’re taking baby steps and exploring slowly. The chain mixes poorly and has high autocorrelation. Conversely, very low acceptance (<15%) means proposals are too aggressive. Optimal acceptance depends on dimension: aim for 20-40% in most cases.

Why it’s tempting: It feels good when proposals are accepted! But acceptance rate is about efficiency, not validity.

Misconception 2: “R-hat < 1.1 guarantees convergence”

Reality: R-hat < 1.1 is necessary but not sufficient for convergence. If your posterior has multiple well-separated modes and each chain gets stuck in one mode, R-hat will look fine but you haven’t converged to the full posterior. You’ve converged to individual modes.

How to catch it: Always visually inspect trace plots and 2D marginal distributions. Look for chains that never overlap. Run chains from diverse starting points spanning the prior range.

Misconception 3: “MCMC samples are independent”

Reality: MCMC samples are correlated. The effective sample size $N_\text{eff} = N/\tau$ accounts for this. If $\tau = 100$, your 10,000 samples carry as much information as 100 independent samples. Always report $N_\text{eff}$, not just $N$.

Practical impact: If you compute credible intervals assuming independence, they’ll be too narrow (overconfident). Use N_eff for uncertainty quantification.

Misconception 4: “Longer chains are always better”

Reality: Running a poorly-tuned chain longer doesn’t fix the problem — it just wastes time. If your proposals are bad, adding more samples doesn’t help. Fix the sampler first (tune step size, use better proposals, add constraints), then run long chains.

The exception: If your chain is well-tuned but exploring a complex posterior, then yes, longer is better. But “longer” means N_eff > 400 per parameter, not just more samples.

Misconception 5: “I can start MCMC anywhere”

Reality: While theory says chains converge from any starting point, starting very far from the typical set (the region of high probability) means long burn-in. Starting outside your prior bounds or at parameter values where the likelihood is undefined will cause immediate failure.

Best practice: Start from reasonable values. For well-specified priors, sample from the prior. For complex problems, run optimization first to find the posterior mode, then start MCMC nearby.

Misconception 6: “Burn-in is waste”

Reality: Burn-in is not waste — it’s how the chain finds the typical set! It’s like a spacecraft’s trajectory correction burns to reach the right orbit. The real waste is keeping burn-in samples and contaminating your posterior estimates.

How much to discard: Be conservative. If in doubt, discard more. Look at trace plots — when do they stabilize?

Misconception 7: “MCMC always works”

Reality: MCMC can fail in many ways:

• Multimodal posteriors where chains get trapped in one mode

• Extremely high-dimensional problems where mixing is glacially slow

• Posteriors with complex geometry (funnel shapes, banana-shaped ridges)

• Highly correlated parameters that random walk proposals can’t handle

• Likelihood functions with discontinuities or numerical instabilities

Recognition: That’s why diagnostics exist! If R-hat is bad, acceptance is extreme, or chains look weird, don’t trust the results. Fix the sampler or use a different method.

Misconception 8: “The posterior mean is the best parameter estimate”

Reality: Depends on your loss function! The posterior mean minimizes squared-error loss. The posterior median minimizes absolute-error loss. The posterior mode (MAP) minimizes 0-1 loss. Which one you report depends on what you’re trying to do. For asymmetric posteriors, they can differ substantially.

Best practice: Report the full posterior (credible intervals) rather than a single point estimate. Let readers make their own decisions about point estimates.

Looking Ahead: Project 4 and Beyond¶

What’s Next¶

In Project 4, you’ll implement everything from this module to measure cosmological parameters from Type Ia supernovae:

Your tasks:

1. Build the forward model (cosmological distance-redshift relation)

2. Implement Metropolis-Hastings MCMC

3. Apply diagnostics (trace plots, R-hat, ACF, corner plots)

4. Run inference on real SNe data (JLA sample)

5. Measure $\Omega_m$ and $h$ (the contents of the universe!)

6. Extend to Hamiltonian Monte Carlo (using your leapfrog integrator from Project 2)

This is the same data and methods that earned the 2011 Nobel Prize in Physics for the discovery of dark energy.

Advanced Topics (Module 5.4 Preview)¶

After mastering Metropolis-Hastings, you’ll learn:

Hamiltonian Monte Carlo (HMC):

• Uses gradient information for efficient exploration

• Based on Hamiltonian dynamics (Module 3!)

• Much lower autocorrelation than random walk

• The leapfrog integrator from Project 2 is exactly what HMC uses

Affine-invariant ensemble samplers (like emcee):

• Multiple “walkers” evolve together

• Self-tuning, no manual proposal tuning needed

• Handles strong correlations automatically

No-U-Turn Sampler (NUTS):

• Adaptive HMC that automatically tunes step size and trajectory length

• The algorithm behind Stan and PyMC

• State of the art for general-purpose MCMC

But these are all variations on the core principles you learned here. Metropolis-Hastings is the foundation.

Where is MCMC Heading? The Frontier of Computational Inference¶

MCMC is not a settled field — it’s actively evolving. Here are some exciting research directions that may transform how we do inference in the next decade:

Neural network-assisted MCMC: Recent work uses machine learning to learn optimal proposal distributions from data. Instead of hand-tuning proposals, train a neural network to predict good proposals based on the posterior landscape. Early results show dramatic speedups for complex problems. This combines classical statistics with modern deep learning.

Adaptive and self-tuning methods: Modern samplers can adapt their behavior during sampling without violating detailed balance (which traditional theory said was impossible!). Methods like adaptive Metropolis (Haario et al. 2001) and delayed rejection algorithms show that we can be smarter about exploration while maintaining theoretical guarantees.

Non-reversible MCMC: Traditional MCMC satisfies detailed balance (microscopic reversibility). But what if we intentionally break reversibility to explore faster? Non-reversible Langevin samplers and lifted particle filters can have much lower autocorrelation than reversible methods. This is cutting-edge research connecting to non-equilibrium statistical mechanics.

Variational inference and hybrid methods: Not everything requires sampling. Variational methods approximate posteriors with simpler distributions (like Gaussians) by optimization rather than sampling. Modern hybrid approaches combine MCMC’s accuracy with variational inference’s speed. The future may be algorithms that switch between sampling and optimization adaptively.

Quantum MCMC: As quantum computers mature, researchers are developing quantum algorithms for Bayesian inference. Quantum walks could explore posterior landscapes exponentially faster than classical MCMC — though practical implementations are years away.

Amortized inference: Instead of running MCMC separately for each dataset, train a model once that can instantly produce posteriors for new data. This “simulation-based inference” approach (also called likelihood-free inference or neural posterior estimation) is revolutionizing fields where likelihood evaluation is expensive.

The common thread: All these advances build on the foundations you learned — Markov chains, stationarity, convergence theory. Understanding Metropolis-Hastings deeply prepares you to understand and contribute to these cutting-edge methods. The principles don’t change; the implementations get ever more sophisticated.

MCMC has come a long way from simulating nuclear reactions on MANIAC I in 1953. Where will it go next? Perhaps you’ll help decide.

The Professional Path¶

After building MCMC from scratch, you’ll be able to use professional tools intelligently:

emcee (Python): Affine-invariant ensemble sampler
PyMC (Python): Full Bayesian modeling framework with HMC/NUTS
Stan (C++): Gold-standard Bayesian inference (autodiff, NUTS)
JAX/numpyro (Python): HMC with Just-In-Time (JIT) compilation and GPU acceleration

You won’t be a black-box user clicking buttons. You’ll understand:

• Why these tools work (detailed balance, ergodicity)

• How to diagnose problems (R-hat, trace plots, divergences)

• When methods fail (multimodal posteriors, stiff likelihoods)

• What the algorithms are actually doing under the hood

That’s the glass-box philosophy in action.

Congratulations! You’ve learned one of the most important computational techniques in modern science. MCMC underpins:

• Bayesian inference across all sciences

• Machine learning (Bayesian neural networks, variational inference)

• Computational physics (lattice QCD, molecular dynamics)

• Statistics (everywhere!)

More importantly, you’ve seen how measurement, probability, physics, and computation all connect into a unified framework for learning from data.

Now go build your inference engine. The universe is waiting to be measured.

Self-Assessment Rubric¶

Before moving forward, honestly assess your mastery of this module. This rubric helps you identify areas for review and confirms readiness for Project 4.

Level 1: Conceptual Understanding¶

Can you explain the ideas to someone else?

• Basic (Developing): I can explain why MCMC is needed in general terms but struggle with the curse of dimensionality details.

• Proficient: I can clearly explain the curse of dimensionality with concrete examples (e.g., 10^20 evaluations for d=10).

• Advanced: I can explain the conceptual progression: crisis (curse) $\to$ solution (Monte Carlo) $\to$ problem (how to sample?) $\to$ answer (Markov chains) $\to$ implementation (Metropolis-Hastings), and justify each step.

Level 2: Mathematical Foundations¶

Can you derive key results?

• Basic (Developing): I can state the detailed balance equation but cannot derive how it implies stationarity.

• Proficient: I can derive that detailed balance implies stationarity by integrating both sides, and I understand why the normalization constant cancels in the acceptance probability.

• Advanced: I can derive the Metropolis-Hastings acceptance criterion from detailed balance, verify it satisfies the ratio condition for both R > 1 and R < 1 cases, and explain why min(1, r) is optimal.

Level 3: Implementation Skills¶

Can you code MCMC from scratch?

• Basic (Developing): I can follow the pseudocode but need significant help implementing it in Python.

• Proficient: I can implement Metropolis-Hastings for a simple 1D problem (like the Cepheid example) with working proposal and acceptance logic, using log-space arithmetic correctly.

• Advanced: I can implement M-H for multidimensional problems, tune proposals adaptively, handle constraints (bounded parameters), and debug common issues (underflow, stuck chains, poor mixing).

Level 4: Diagnostic Expertise¶

Can you tell if your sampler is working?

• Basic (Developing): I know I should check trace plots but I’m not sure what “good” looks like.

• Proficient: I can create and interpret trace plots, compute R-hat and ACF, determine when burn-in is sufficient, and recognize clear convergence failures (stuck chains, high R-hat).

• Advanced: I can diagnose subtle problems (multimodal posteriors where chains are in different modes but R-hat looks okay, systematic biases in proposals), adjust sampling strategy accordingly, and justify my burn-in and thinning choices quantitatively.

Level 5: Scientific Application¶

Can you apply MCMC to real problems?

• Basic (Developing): I can run provided MCMC code but struggle to set up likelihoods and priors for new problems.

• Proficient: I can formulate a Bayesian inference problem from scratch (define likelihood from measurement model, choose appropriate priors, implement log-posterior), run MCMC, and interpret credible intervals correctly.

• Advanced: I can handle complex scientific problems (multiple parameters with constraints, multimodal posteriors, expensive likelihoods requiring optimization), choose appropriate diagnostics for the specific problem, and present results in publication-ready format with proper uncertainty quantification.

Level 6: Connections and Transfer¶

Can you see how MCMC fits into the bigger picture?

• Basic (Developing): I see MCMC as an isolated technique for Bayesian inference.

• Proficient: I understand how MCMC connects to Module 1 (CLT, LLN, ergodicity) and Module 3 (detailed balance = microscopic reversibility, Boltzmann distribution analogy).

• Advanced: I can articulate the deep mathematical unity across modules (same stochastic process theory governs sampling, molecular dynamics, and MCMC), explain how different MCMC variants (HMC, ensemble samplers) relate to the core principles, and identify when MCMC is the right tool versus alternatives (optimization, variational inference, grid methods).

Remember: This is self-assessment for learning, not grading. Be honest with yourself. Identifying gaps now makes Project 4 more productive.

Next Module: Hamiltonian Monte Carlo (HMC) and advanced sampling methods