32.3 Model Comparison and Decision Making

1. The Widely Applicable Information Criterion (WAIC)

Two climate models fit last century's temperature record equally well, yet one predicts catastrophic warming by 2100 and the other predicts gentle drift; which should a policymaker trust, and how would you know? That question sits at the heart of model comparison: we want to choose the model that will make the best predictions on future data drawn from the same process. The expected log pointwise predictive density (ELPD) measures exactly this:

WAIC is a scoring rule that estimates how well a Bayesian model will predict new, unseen observations, using only the posterior samples you already have. Choosing between models based on training-data fit alone is misleading. A complex model can memorize noise and perform terribly on fresh data. WAIC corrects for this by subtracting a data-driven penalty that grows with effective model complexity. The mechanism works as follows: average each observation's likelihood across all posterior draws to get a raw fit score, then subtract the variance of that log-likelihood across draws. Observations whose likelihood swings wildly across parameter values signal overfitting, so they incur a larger penalty. Use WAIC when you need a quick, posterior-sample-based estimate of predictive performance and your model is reasonably well-specified. Prefer LOO-CV (Section 2) when you need per-observation diagnostics or suspect influential outliers. Prefer Bayes factors (Section 3) when the scientific question concerns which hypothesis is true rather than which model predicts best.

We cannot compute ELPD directly because we do not know the true data-generating distribution \(p_{\text{true}}\). WAIC (Watanabe, 2010) estimates ELPD from the posterior predictive distribution:

The first term is the log pointwise predictive density (lppd), where \(\hat{p}(y_i \mid D) = \frac{1}{S} \sum_{s=1}^{S} p(y_i \mid \theta^{(s)})\) averages the likelihood of each observation over posterior samples. The second term \(V_i = \text{Var}_{s}[\log p(y_i \mid \theta^{(s)})]\) is a bias correction that penalizes model complexity: observations that are highly sensitive to the particular parameter values (high variance across posterior samples) contribute more to the penalty. This penalty is the Bayesian analogue of the parameter count in AIC.

Computing WAIC in Practice

WAIC is computed entirely from posterior samples, making it straightforward to calculate after any MCMC or VI run. Higher WAIC (less negative ELPD) means better expected out-of-sample prediction. In short: the model that best predicts data it has not yet seen is the one that has learned the pattern, not memorized the noise.

2. LOO-CV with Pareto Smoothed Importance Sampling

LOO-CV is the gold standard for estimating out-of-sample prediction: fit the model \(n\) times, each time leaving out one observation, and measure how well the model predicts the left-out point. The LOO-CV estimate of ELPD is:

where \(D_{-i}\) is the dataset with observation \(i\) removed. Naive LOO-CV requires refitting the model \(n\) times, which is computationally prohibitive for MCMC-based inference. Vehtari, Gelman, and Gabry (2017) showed that Pareto Smoothed Importance Sampling (PSIS) can approximate LOO-CV from a single posterior fit. PSIS reweights the existing posterior samples to approximate the leave-one-out posteriors, then stabilizes the importance weights by fitting a generalized Pareto distribution (a flexible family for modeling the shape of distribution tails) to the largest weights. (For a dataset of 1,000 observations, that is a 1,000x speedup: one posterior fit instead of one thousand.)

The quality of the PSIS approximation is diagnosed by the Pareto \(\hat{k}\) parameter for each observation. If \(\hat{k} < 0.7\), the approximation is reliable. If \(\hat{k} > 0.7\), that observation has too much influence on the posterior, and you need to actually refit the model without that observation (or consider whether the observation is an outlier that warrants investigation, connecting to the anomaly detection methods of Chapter 30).

import pymc as pm
import arviz as az
import numpy as np

# Two competing models for the same data
rng = np.random.default_rng(42)
n = 40
x = rng.uniform(0, 10, size=n)
y_true = 3.0 * np.sin(0.5 * x) + rng.normal(0, 0.8, size=n)

# Model 1: Linear relationship
with pm.Model() as linear_model:
    a = pm.Normal("a", mu=0, sigma=5)
    b = pm.Normal("b", mu=0, sigma=5)
    sigma = pm.HalfNormal("sigma", sigma=3)
    mu = a + b * x
    pm.Normal("y", mu=mu, sigma=sigma, observed=y_true)
    trace_linear = pm.sample(2000, tune=1000, chains=4, random_seed=42)

# Model 2: Sinusoidal relationship (closer to truth)
with pm.Model() as sine_model:
    amp = pm.Normal("amp", mu=0, sigma=5)
    freq = pm.HalfNormal("freq", sigma=2)
    sigma = pm.HalfNormal("sigma", sigma=3)
    mu = amp * pm.math.sin(freq * x)
    pm.Normal("y", mu=mu, sigma=sigma, observed=y_true)
    trace_sine = pm.sample(2000, tune=1000, chains=4, random_seed=42)

# LOO-CV comparison
loo_linear = az.loo(trace_linear, linear_model)
loo_sine = az.loo(trace_sine, sine_model)

# Compare models: ranks by ELPD, reports differences with standard errors
comparison = az.compare(
    {"linear": trace_linear, "sinusoidal": trace_sine},
    ic="loo",
)
print(comparison)

# Check Pareto k diagnostics for the winning model
print(f"\nMax Pareto k (sinusoidal): {loo_sine.pareto_k.values.max():.3f}")
print("Values > 0.7 indicate influential observations needing investigation.")

3. Bayes Factors

The Bayes factor directly compares the marginal likelihoods (where the marginal likelihood of a model is the probability of the observed data averaged over all parameter values under that model's prior) of two models:

A Bayes factor of 10 means the data are 10 times more likely under Model 1 than Model 2. Kass and Raftery (1995) provide an interpretation scale: \(\text{BF} > 3.2\) is "substantial" evidence, \(> 10\) is "strong," and \(> 100\) is "decisive." The Bayes factor embodies the Bayesian Occam's razor: a simpler model (narrower prior) that explains the data gets a higher marginal likelihood than a complex model (broad prior) that can explain anything but concentrates less prior mass on the actual observations.

Computing Bayes factors is harder than computing LOO-CV because it requires the marginal likelihood, which is an integral over the entire parameter space. Several approximation methods exist:

• Bridge sampling (Meng and Wong, 1996): uses an auxiliary "bridge" distribution to connect posterior samples from both models. Available in the bridgesampling R package and Python implementations.
• Savage-Dickey ratio: for nested models where one model fixes a parameter \(\psi = \psi_0\) that the other estimates freely, the Bayes factor equals the ratio of the posterior to the prior density at \(\psi_0\). Elegant when applicable.
• Harmonic mean estimator: tempting but numerically unstable; the variance is often infinite. Avoid in practice.

import pymc as pm
import numpy as np
from scipy.stats import norm, gaussian_kde

# Simulated experiment: small effect with moderate noise
rng = np.random.default_rng(42)
effect_true = 0.35     # small but real effect
data = rng.normal(effect_true, 1.0, size=25)

with pm.Model() as effect_model:
    delta = pm.Normal("delta", mu=0, sigma=1)     # prior centered at 0
    sigma = pm.HalfNormal("sigma", sigma=2)
    pm.Normal("obs", mu=delta, sigma=sigma, observed=data)
    trace = pm.sample(4000, tune=2000, chains=4, random_seed=42)

# Savage-Dickey ratio: BF for delta = 0 (null) vs delta != 0
posterior_samples = trace.posterior["delta"].values.flatten()
prior_at_zero = norm.pdf(0, loc=0, scale=1)  # Normal(0,1) density at 0

kde = gaussian_kde(posterior_samples)
posterior_at_zero = kde.evaluate(0.0)[0]

bf_null = posterior_at_zero / prior_at_zero
bf_alt = 1.0 / bf_null

print(f"Prior density at delta=0:     {prior_at_zero:.4f}")
print(f"Posterior density at delta=0:  {posterior_at_zero:.4f}")
print(f"BF (null vs alternative):     {bf_null:.3f}")
print(f"BF (alternative vs null):     {bf_alt:.3f}")
# BF_alt > 3 suggests substantial evidence for a nonzero effect

4. WAIC vs. LOO-CV vs. Bayes Factors

These three approaches answer subtly different questions, and choosing among them depends on the scientific context, as summarized in Table 32.3:

In practice, LOO-CV is the default recommendation (Vehtari et al., 2017). It is easy to compute from posterior samples, provides per-observation diagnostics, and answers the practical question of which model predicts better. Bayes factors are appropriate when the scientific question is genuinely about comparing theories (does this effect exist? is this physical law the correct one?) rather than about prediction.

5. Calibration: Does Your Uncertainty Match Reality?

Selecting the best model is only half the job; you also need to verify that the winning model's stated uncertainty is trustworthy.

A model that reports 95% credible intervals should contain the true value 95% of the time. If it does, the model is well-calibrated. If 95% intervals contain the truth only 80% of the time, the model is overconfident (reporting narrower uncertainty than it should). If they contain the truth 99% of the time, the model is underconfident (wasting information by reporting wider uncertainty than necessary).

The standard tool for checking calibration is the probability integral transform (PIT). For each observation \(y_i\), compute the posterior predictive cumulative distribution function (CDF), where the CDF at a value is the probability that a new draw from the model would fall at or below that value, evaluated at \(y_i\): \(u_i = P(Y_i \leq y_i \mid D)\). If the model is well-calibrated, these PIT values should be uniformly distributed on \([0, 1]\) (because under a correct model, the probability of an observation falling below any given quantile is, by definition, equal to that quantile). Deviations from uniformity diagnose specific problems: a U-shaped PIT histogram indicates that the model's predictive intervals are too narrow (observations land in the tails more often than predicted), and an inverse U-shape indicates that the model's predictive intervals are too wide (observations cluster near the center).

import arviz as az
import numpy as np

def check_calibration(trace, model, observed, var_name="y"):
    """Check calibration via PIT histogram and coverage analysis.

    Args:
        trace: InferenceData from MCMC sampling
        model: PyMC model context
        observed: array of observed values
        var_name: name of the observed variable in the model
    Returns:
        pit_values: array of PIT values for each observation
    """
    with model:
        ppc = pm.sample_posterior_predictive(trace, random_seed=42)

    # Compute PIT values: fraction of posterior predictive samples <= observed
    ppc_samples = ppc.posterior_predictive[var_name].values
    ppc_flat = ppc_samples.reshape(-1, len(observed))  # (total_draws, n_obs)

    pit_values = np.mean(ppc_flat <= observed[np.newaxis, :], axis=0)

    # Check coverage at multiple nominal levels
    print("Calibration report:")
    for level in [0.50, 0.80, 0.90, 0.95]:
        alpha = 1 - level
        lower = np.quantile(ppc_flat, alpha / 2, axis=0)
        upper = np.quantile(ppc_flat, 1 - alpha / 2, axis=0)
        coverage = np.mean((observed >= lower) & (observed <= upper))
        gap = abs(coverage - level)
        status = "GOOD" if gap < 0.05 else "WARN" if gap < 0.10 else "BAD"
        print(f"  {level*100:.0f}% interval: coverage = {coverage:.3f} [{status}]")

    # ArviZ provides the PIT plot as a one-liner
    az.plot_loo_pit(idata=trace, y=var_name)

    return pit_values

6. Value of Information: Which Experiment Should We Run Next?

Once calibration confirms that the model's uncertainty reflects reality, a natural next question arises: where should we look to reduce that uncertainty most efficiently?

The expected information gain (EIG) answers a question that is central to scientific discovery: given our current knowledge (the posterior), which possible experiment would teach us the most? The EIG of an experiment \(\xi\) is the expected reduction in posterior entropy:

This is the expected KL divergence (a non-negative measure of how much one probability distribution differs from another; zero only when the two distributions are identical) between the posterior we would have after observing the outcome \(y\) and the posterior we have now. Experiments with high EIG are those whose outcomes would substantially change our beliefs, regardless of what those outcomes turn out to be.

Computing EIG exactly requires marginalizing over all possible experimental outcomes, an integral over the data space. For models that support sampling from the posterior predictive distribution (the distribution of new observations implied by the posterior over parameters), Monte Carlo provides a practical estimator:

import numpy as np
from scipy.stats import norm

def expected_information_gain(
    posterior_samples, likelihood_fn, design_points, n_mc=200, rng=None
):
    """Estimate expected information gain for candidate experimental designs.

    Uses the nested Monte Carlo estimator of Myung et al. (2013).

    Args:
        posterior_samples: (n_post, n_params) draws from current posterior
        likelihood_fn: callable (theta, x_design) -> (mean, std) of observation
        design_points: array of candidate experimental conditions
        n_mc: number of synthetic observations per design point
        rng: numpy random generator
    Returns:
        eig: array of EIG estimates, one per design point
    """
    if rng is None:
        rng = np.random.default_rng(0)
    n_post = len(posterior_samples)
    eig = np.zeros(len(design_points))

    for d_idx, x_d in enumerate(design_points):
        mi_samples = []

        for _ in range(n_mc):
            # 1. Draw a "true" theta from the posterior
            idx = rng.integers(n_post)
            theta_star = posterior_samples[idx]
            mu_star, sigma_star = likelihood_fn(theta_star, x_d)

            # 2. Simulate an observation under theta_star
            y_sim = rng.normal(mu_star, sigma_star)

            # 3. Evaluate log p(y_sim | theta) for all posterior samples
            log_likes = np.array([
                norm.logpdf(y_sim, *likelihood_fn(theta, x_d))
                for theta in posterior_samples
            ])

            # 4. Log marginal: log (1/S) sum_s p(y_sim | theta_s)
            log_marginal = np.logaddexp.reduce(log_likes) - np.log(n_post)

            # 5. Log conditional: log p(y_sim | theta_star)
            log_cond = norm.logpdf(y_sim, mu_star, sigma_star)

            mi_samples.append(log_cond - log_marginal)

        eig[d_idx] = np.mean(mi_samples)

    return eig

# Example: which temperature to measure next for a reaction rate study?
candidate_temps = np.linspace(250, 500, 50)  # Kelvin
# eig = expected_information_gain(posterior_samples, arrhenius_likelihood, candidate_temps)
# optimal_temp = candidate_temps[np.argmax(eig)]
# print(f"Most informative next measurement: {optimal_temp:.0f} K")

7. The Decision-Theoretic Discovery Cycle

The complete Bayesian discovery cycle integrates model building, comparison, calibration, and experimental design into a coherent loop. Figure 32.4 illustrates the six stages of this cycle, showing how each stage feeds into the next. Figure 32.3.1 illustrates Bayesian decision-theoretic discovery cycle.

1. Formulate competing hypotheses as Bayesian models (Section 32.1).
2. Compute posteriors for each model (Section 32.2).
3. Compare models using LOO-CV to identify the best-supported hypothesis.
4. Check calibration to verify that the winning model's uncertainty is trustworthy.
5. Compute EIG to determine the most informative next experiment.
6. Run the experiment, update the data, and return to step 2.

This loop is the quantitative backbone of the automated experiment design systems in Chapter 46 and the self-driving laboratories of Chapter 55, where it runs autonomously: a robot executes experiments while a Bayesian model decides what to try next.

Exercises

1. (Conceptual) A Bayes factor of 0.5 for Model A versus Model B means the data moderately favor Model B. But a LOO-CV comparison shows Model A has higher ELPD. Explain how these two criteria can disagree. Which should you trust, and under what circumstances? Hint: consider the role of priors in each criterion.
2. (Coding) Generate data from a quadratic model \(y = 2 + 0.5x + 0.1x^2 + \epsilon\) with \(n = 30\) points and \(\epsilon \sim \text{Normal}(0, 1)\). Fit three PyMC models (linear, quadratic, cubic) and compare them using az.compare with LOO-CV. How does the ranking change as you increase \(n\) to 100 and 500? At what sample size does LOO-CV reliably identify the quadratic model as best?
3. (Analysis) Using the EIG estimation function from this section, design an optimal 5-point experiment for a simple linear regression model \(y = a + bx + \epsilon\). Start with a prior \(a \sim \text{Normal}(0, 10)\), \(b \sim \text{Normal}(0, 10)\), \(\sigma = 1\), and candidate design points \(x \in \{-5, -4, \ldots, 4, 5\}\). After each observation (simulated from the true model with \(a=1, b=2\)), update the posterior and recompute EIG. Where does the optimal design place its points, and how does this compare to classical D-optimal design (where D-optimality selects design points that maximize the determinant of the Fisher information matrix, concentrating measurements at the extremes of the design space)?