Part 3: Markov Chain Monte Carlo — From Theory to Practice

“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:

1. Explain why direct evaluation of posteriors fails in high dimensions (curse of dimensionality)

2. Derive the Metropolis-Hastings acceptance criterion from detailed balance

3. Recognize the connection between MCMC sampling and ergodic exploration from Module 1

4. Implement a Markov chain that provably converges to any target distribution

5. Diagnose convergence using trace plots, autocorrelation, and the Gelman-Rubin statistic

6. Apply MCMC to real astronomical inference problems (Cepheid distance moduli)

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

Roadmap: Where We’re Going¶

This is a long module because we’re building something profound: a universal inference engine that works for any scientific model, any data, any dimensionality. Here’s the journey:

The narrative arc: We’ll start with a crisis (can’t evaluate high-dimensional integrals), discover a profound solution (sampling via Markov chains), derive the mathematical machinery (detailed balance), build the algorithm (Metropolis-Hastings), learn how to use it (diagnostics), and apply it to real science (Cepheid standard candles preparing for Type Ia SNe).

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|θ)p(θ)dθ$ is the evidence (normalization constant)

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

Let’s be concrete. Suppose you’re measuring the distance to a Cepheid variable star. Your parameter is $θ = d$ (distance in parsecs). You have:

• Data D: Observed apparent magnitude $m = 15.2 ± 0.1$ mag

• Model: Absolute magnitude $M = -4.5$ (from period-luminosity relation)

• Likelihood: Gaussian, $p(m|d) = N(m_\text{theory}(d), σ²)$ where $m_\text{theory(d)} = M + 5\log_{10}(d) - 5$

The posterior $p(d|m)$ 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 low dimensions (1-2 parameters), you could just evaluate $p(d|m)$ on a grid of distance values and plot it. Done! But what if you have 10 parameters? Or 100?

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 × 100 = 10^4$ evaluations.

3D example: Three parameters → $100 × 100 × 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(θ|D)$ and report that!”

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 $θ = 5.0$ as the MAP estimate, can you distinguish whether your uncertainty is $±0.1$ or $±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 ± 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 doesn’t report “best fit” values. It reports credible intervals from the full posterior distribution. This is what p-values tried to do (badly). Bayesian posteriors do it right.

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

Here’s a profound insight from the early days of computing (1940s-1950s): You don’t need to evaluate the posterior everywhere. You need to sample from it.

If you have N samples ${θ₁, θ₂, ..., θ_N}$ drawn from $p(θ|D)$, then you can estimate anything:

Mean:

Variance:

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

Predictions: For any function f(θ), estimate its expectation:

The magic: The accuracy of Monte Carlo estimates scales as $\sigma/\sqrt{N}$, where $\sigma$ is the standard deviation. The variance of Monte Carlo estimates scales as $1/\sqrt{N}$ regardless of dimensionality — this follows directly from the Central Limit Theorem you studied in Module 1.

To get one more decimal place of accuracy, you need $100×$ more samples. But crucially: this doesn’t depend on dimensionality! Going from 10 to 100 dimensions doesn’t change this. It’s the same $\sqrt{N}$ convergence.

This is revolutionary. Grid evaluation scales exponentially with dimension $(N^d)$. Monte Carlo scales with accuracy $(1/\sqrt{N})$, independent of dimension.

But here’s the catch: How do you draw samples from $p(θ|D)$ when you don’t know what $p(θ|D)$ looks like? You can evaluate it at any point (numerator of Bayes’ theorem), but sampling requires knowing the shape.

This is where Markov Chain Monte Carlo enters. It’s a method to generate samples from any target distribution, no matter how complex, by constructing a Markov chain that explores parameter space.

Part 2: Markov Chains — The Mathematical Foundation¶

The Core Idea: Build a Stochastic Process¶

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

1. Starts from any initial $θ₀$

2. Takes random steps through parameter space

3. Gradually “forgets” where it started

4. Eventually produces samples distributed according to $p(θ|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(θ|D)$.

If $p(θ|D)$ is high in some region, the walker visits often. If $p(θ|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 a sequence of random variables $\{θ₀, θ₁, θ₂, ...\}$ where the probability of the next state depends only on the current state:

This is the Markov property: The future depends on the present, not the past. The chain has “no memory” beyond its current location.

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 $θ'$ given that you’re currently at $θ$. It must satisfy:

for all $θ$ (probabilities must sum to 1).

Given a starting distribution $p₀(θ)$, the distribution after one step is:

After $t$ steps:

Question: Can we design $T$ such that p_t(θ) converges to our target distribution π(θ) = p(θ|D) as t → ∞?

Answer: Yes! But we need two conditions:

1. Stationarity: The target distribution $π$ 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 $π(θ)$ 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 $π$, then after one more step, it’s still distributed according to $π$. 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.

Ergodicity: The Chain Must Explore Everywhere¶

Stationarity means $π$ is preserved if you start from it. But how do we get there from an arbitrary starting point $θ₀$?

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 π, then for any function f(θ):

Time averages along the chain converge to expectations under $π$. 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¶

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

In words: The probability flow from θ to θ’ equals the flow from θ’ back to θ when both are drawn from π.

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 θ:

The right side factors since π(θ’) 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 θ’ after one transition starting from π. Since LHS = RHS = π(θ’), we’ve shown that starting from distribution π produces distribution π after one step. The distribution is stationary.

Now we have our blueprint:

1. Choose a target distribution $π(θ) = p(θ|D)$.

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

3. The Markov chain will converge to $π$ (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(θ'|θ)$: Generates a candidate new state θ’ given current state θ

2. Acceptance probability $α(θ'|θ)$: Probability of accepting the proposal

The full transition kernel is:

The algorithm:

1. Start at some $θ₀$ (anywhere in parameter space)

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

   • Propose: Draw $θ'$ from $Q(θ'|θ_t)$

   • Evaluate: Compute acceptance probability $α(θ'|θ_t)$

   • Decide:

     • With probability $α$, accept: $θ_{t+1} = θ'$

     • With probability $1-α$, reject: $θ_{t+1} = θ_t$ (stay put)

3. Return samples ${θ₁, θ₂, ..., θ_N}$

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

Deriving the Acceptance Probability¶

We want:

Substituting $T = Qα$ where $Q$ is the proposal and $α$ 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 $α(θ'|θ)$ and $α(θ|θ')$ to have ratio $r$. A simple choice that works:

Verification: Let

• If $R ≥ 1: α(θ'|θ) = 1, α(θ|θ') = 1/R ≤ 1$. Ratio: $1/(1/R) = R$ ✓

• If $R < 1: α(θ'|θ) = R, α(θ|θ') = 1$. Ratio: $R/1 = R$ ✓

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

Interpretation:

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

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

This means uphill moves are always accepted, downhill moves are sometimes accepted. The chain preferentially explores high-probability regions but can escape local optima by occasionally accepting downward moves.

Proposal Distributions: Symmetric vs. Asymmetric¶

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

Symmetric proposals: $Q(θ'|θ) = Q(θ|θ')$

Example: Gaussian random walk:

Propose by adding Gaussian noise: $θ' = θ + ε \text{ where } ε \sim N(0, Σ)$. Here $\mathcal{N}(0, Σ)$ is the normal (Gaussian) distribution with mean 0 and covariance $\Sigma$.

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

Example — Gaussian random walk (symmetric in $θ$):

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

Here you can evaluate $π$ 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(θ'|θ) ≠ Q(θ|θ')$

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

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 $(α ≈ 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 $(α ≈ 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 $σ$ by 1.2

• If acceptance rate < 20%, divide $σ$ 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.)

Pseudocode: The Complete Algorithm¶

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 π(θ) (log posterior probability). Return the unnormalized log posterior: log p(D|θ) + log p(θ). Do not include log p(D); it is a constant and cancels in MH.
    proposal : function  
        Generates θ' given θ. Returns θ’ and log_Q_ratio = log Q(θ|θ') − log Q(θ'|θ) (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: Posteriors often involve tiny numbers (like 10⁻¹⁰⁰⁰). Directly computing exp(-1000) causes underflow — it becomes exactly 0 in floating point arithmetic. But log-space arithmetic handles this gracefully:

# 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 * log(10)  # ~-2302.6
log_p2 = -1001 * 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: For symmetric Gaussian:

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(θ|θ')/Q(θ'|θ) = 1
    return theta_new, log_Q_ratio

3. Acceptance decision: Compare log(uniform random) to log $α$ instead of comparing uniform random to $α$. 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 × 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|θ)$ (your likelihood/model) and $p(θ)$ (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 π. 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 π. You need to discard them.

Visual Diagnostic: Trace Plots¶

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

A trace plot shows $θ_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 θ_t is with θ_{t+k}:

where $μ = ⟨θ⟩$ is the chain mean.

Interpretation:

• $ρ₀ = 1$ (perfect correlation with itself)

• $ρ_k → 0$ as $k → ∞$ for an ergodic chain

• Larger $k$ needed for $ρ_k ≈ 0$ means higher correlation (worse mixing)

Autocorrelation time τ: The characteristic lag at which samples become independent:

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

Effective sample size (ESS):

If $τ = 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 $τ < 100$ (or you’ve run long enough that $N/τ > 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.

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.

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 ± 0.15 mag (from photometry)

What you want to infer:

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

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

The forward model:

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

The distance modulus is:

So:

Solving for $d$:

But this is just a point estimate. What’s the uncertainty? Let’s use Bayesian inference.

Bayesian Setup¶

Parameters: $θ =$ distance $d$ (in kpc for numerical convenience)

Data: Observed apparent magnitude$m_\text{obs} = 18.50$ mag

Likelihood: Assuming Gaussian measurement errors:

where:

• $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)

In log form:

Prior: Uniform in log-space (Jeffrey’s prior for a scale parameter):

We’ll actually sample in log(d) to make this 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).

Posterior: By Bayes’ theorem:

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.

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_eff = N/τ accounts for this. If τ = 100, your 10,000 samples carry as much information as 100 independent samples. Always report N_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 $Ωₘ$ 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 (training deep networks is MCMC-like)

• 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) → solution (Monte Carlo) → problem (how to sample?) → answer (Markov chains) → 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