Chapter 38: Bayesian Inference

The Mechanics of Belief Updating

Bayes' theorem is a direct consequence of the definition of conditional probability. If $A$ and $B$ are events with $P(B) > 0$:

In a statistical inference context, replace $A$ with $\theta$ (a parameter or hypothesis) and $B$ with $\mathcal{D}$ (observed data):

Each term has a name and a role:

• Prior <!--MATHBLOCK14-->: your distribution over <!--MATHBLOCK15--> before seeing data. Encodes domain knowledge, physical constraints, or deliberate ignorance.
• Likelihood <!--MATHBLOCK16-->: the probability of the observed data as a function of <!--MATHBLOCK17-->. This is identical to the likelihood used in maximum likelihood estimation (MLE).
• Posterior <!--MATHBLOCK18-->: your updated distribution over <!--MATHBLOCK19--> after observing data. This is the inferential target.
• Marginal likelihood (evidence) <!--MATHBLOCK20-->: a normalizing constant ensuring the posterior integrates to one. It is crucial for model comparison but irrelevant for single-model inference, so we often write the unnormalized form:

A Worked Medical Example

Consider a disease that affects 1% of the population. A diagnostic test has 90% sensitivity (true positive rate) and 95% specificity (true negative rate). You test positive. What is the probability you have the disease?

Let $D$ = "has disease", $T^+$ = "tests positive".

Despite a highly accurate test, a positive result only raises your probability from 1% to about 15%. This result — often called the base rate fallacy when ignored — is the canonical demonstration of why prior information matters profoundly. Frequentist hypothesis testing has no direct mechanism to incorporate this base rate.

# Bayesian updating for the medical screening example
prior_disease = 0.01
sensitivity = 0.90       # P(T+ | D)
specificity = 0.95       # P(T- | not D)
false_positive_rate = 1 - specificity

# P(T+)
p_positive = sensitivity * prior_disease + false_positive_rate * (1 - prior_disease)

# Posterior via Bayes
posterior = (sensitivity * prior_disease) / p_positive
print(f"Prior P(disease)     : {prior_disease:.3f}")
print(f"P(T+)                : {p_positive:.4f}")
print(f"Posterior P(D | T+)  : {posterior:.3f}")

Frequentist vs. Bayesian: The Core Distinction

The philosophical divide comes down to the interpretation of probability itself.

A frequentist defines probability as the limiting relative frequency of outcomes over infinitely many independent repetitions of an experiment. Parameters are fixed constants — asking for "the probability that $\mu = 5$" is meaningless because $\mu$ either is or is not 5. Inference is about constructing procedures (estimators, tests, confidence intervals) with guaranteed long-run properties.

A Bayesian defines probability as a measure of rational belief or uncertainty, applicable to any proposition including parameter values. Probabilities are personal in the sense that they encode a rational agent's state of knowledge given available information. Inference is about coherently updating beliefs via Bayes' theorem.

In practice, the differences surface in:

• Interval estimation: A 95% frequentist confidence interval is a random interval with the property that 95% of all such intervals constructed from repeated samples would contain the true parameter. A single interval either contains the parameter or it does not — no probability attaches to the individual interval. A 95% Bayesian credible interval directly states that, given the data and prior, there is a 95% posterior probability the parameter lies in that interval. The Bayesian statement is what most users mistakenly believe about confidence intervals.
• Sample size and stopping rules: Bayesian inference is valid regardless of when you stop collecting data. Frequentist <!--MATHBLOCK25-->-values are distorted by optional stopping, which creates practical problems in sequential experiments.
• Incorporating prior knowledge: Bayesian priors provide a formal channel for domain expertise. Frequentist methods have no analogous mechanism (though regularization methods like ridge regression have implicit Bayesian interpretations).

Choosing a Prior

The prior distribution $P(\theta)$ is the most distinctive and debated element of Bayesian inference. Critics argue it introduces subjectivity; practitioners counter that all statistical analyses embed assumptions, and Bayesian priors make those assumptions explicit and auditable. In practice, prior choice is guided by three considerations: (1) domain knowledge about plausible parameter ranges, (2) mathematical convenience, and (3) theoretical desiderata like invariance and non-informativeness.

Types of Priors

Informative priors encode substantive knowledge. If you are estimating the effect of a drug that is chemically similar to others with effect sizes of 0.1–0.3 standard deviations, a prior $\mathcal{N}(0.2, 0.05^2)$ is defensible and will shrink small-sample estimates toward the historically sensible region.

Weakly informative priors (the workhorse of applied Bayesian analysis) provide gentle regularization without strong claims. For a logistic regression coefficient, a $\mathcal{N}(0, 1)$ prior on the log-odds scale says effect sizes much beyond ±2 are unlikely — reasonable for most real-world problems without being dogmatic.

Non-informative / reference priors attempt to express ignorance. The uniform prior $P(\theta) \propto 1$ is the simplest, but is problematic because uniformity is not invariant under reparameterization (a uniform prior on $\sigma$ implies a non-uniform prior on $\sigma^2$). Jeffreys prior $P(\theta) \propto \sqrt{I(\theta)}$, where $I(\theta)$ is the Fisher information, is invariant under reparameterization and is the standard non-informative choice.

Conjugate Priors

A prior is conjugate to a likelihood if the posterior belongs to the same distributional family as the prior. Conjugacy yields closed-form posterior updates — no numerical integration required — making it invaluable for pen-and-paper analysis, fast computation, and building intuition.

Beta-Binomial (Proportion Estimation)

The canonical conjugate pair. If $X \mid \theta \sim \text{Binomial}(n, \theta)$ and $\theta \sim \text{Beta}(\alpha, \beta)$, then:

The parameters $\alpha$ and $\beta$ have a natural interpretation as pseudo-counts: $\alpha - 1$ prior successes and $\beta - 1$ prior failures. The posterior mean is:

This is a weighted average between the prior mean $\alpha/(\alpha+\beta)$ and the MLE $k/n$, with the prior's weight proportional to $\alpha + \beta$ (the prior's effective sample size). As $n \to \infty$, the MLE dominates — the prior is washed out by data.

from scipy.stats import beta as beta_dist
import numpy as np

# Three priors: uniform, skeptical (low base rate), optimistic
priors = [
    ("Uniform Beta(1,1)",     1,  1),
    ("Skeptical Beta(2,8)",   2,  8),
    ("Optimistic Beta(8,2)",  8,  2),
]

# Observed data: 6 successes out of 10 trials
k, n = 6, 10

print(f"{'Prior':<25} {'Post Mean':>10} {'95% CrI':>22}")
print("-" * 60)
for name, a, b in priors:
    a_post = a + k
    b_post = b + (n - k)
    rv = beta_dist(a_post, b_post)
    lo, hi = rv.ppf(0.025), rv.ppf(0.975)
    print(f"{name:<25} {rv.mean():>10.3f}   [{lo:.3f}, {hi:.3f}]")

Normal-Normal (Mean Estimation with Known Variance)

If $X_1, \ldots, X_n \mid \mu \sim \mathcal{N}(\mu, \sigma^2)$ with known $\sigma^2$, and $\mu \sim \mathcal{N}(\mu_0, \tau_0^2)$, then:

where the precision-weighted posterior mean and variance are:

Again, the posterior mean is a precision-weighted average of the prior mean and the sample mean. High prior precision (small $\tau_0^2$) anchors the estimate near the prior; high data precision (large $n/\sigma^2$) pulls it toward $\bar{x}$.

Gamma-Poisson (Rate Estimation)

If $X \mid \lambda \sim \text{Poisson}(\lambda)$ and $\lambda \sim \text{Gamma}(\alpha, \beta)$ (shape-rate parameterization), then:

The rate parameter $\beta$ acts as the prior's "exposure" — it is the pseudo-count of time intervals observed before any data.

Common Pitfalls in Prior Specification

• Diffuse priors on constrained parameters: A <!--MATHBLOCK32--> prior on a standard deviation seems vague, but since <!--MATHBLOCK33-->, half the prior mass is useless. Use half-Normal or Exponential priors for positive quantities.
• Prior-data conflict: If the data strongly contradicts the prior, the posterior can be bimodal or behave unexpectedly. Always run a prior predictive check — simulate data from the prior and verify it covers the plausible range of observations.
• Ignoring parameterization: A uniform prior on <!--MATHBLOCK34--> is a <!--MATHBLOCK35--> (Jeffreys) prior on <!--MATHBLOCK36--> itself. Know which scale you are placing your prior on.

From Posterior Distribution to Actionable Summaries

The posterior distribution $P(\theta \mid \mathcal{D})$ is a complete description of uncertainty about $\theta$ after observing data. In practice, you often need to reduce this distribution to a few interpretable numbers. This section covers the main posterior summaries and their uses.

Point Estimates

Three natural point estimates emerge from the posterior:

• Posterior mean <!--MATHBLOCK4-->: minimizes expected squared error loss. The workhorse for continuous parameters.
• Posterior median: minimizes expected absolute error loss. More robust to skewed posteriors and heavy tails.
• Maximum a posteriori (MAP) estimate <!--MATHBLOCK5-->: the posterior mode. Equivalent to regularized MLE (e.g., ridge regression uses a Gaussian prior implicitly). Not invariant under reparameterization, and can be highly misleading for multimodal or skewed posteriors — the mode of a distribution is not generally representative of it.

For symmetric unimodal posteriors (like a well-identified Gaussian), all three coincide. For skewed posteriors (common with small samples, bounded parameters, or hierarchical models), they can differ substantially.

Types of Credible Intervals

A credible interval (CrI) $[L, U]$ satisfies:

Unlike a frequentist confidence interval, this statement has a direct probability interpretation conditional on the observed data. Two common constructions:

Equal-Tailed Interval (ETI)

Place $\alpha/2$ probability in each tail. For a 95% ETI, compute the 2.5th and 97.5th percentiles of the posterior. This is easy to compute from MCMC samples:

import numpy as np
from scipy.stats import beta as beta_dist

# Posterior: Beta(16, 8) from Beta-Binomial update
rv = beta_dist(16, 8)

# 95% Equal-Tailed Interval
eti_low, eti_high = rv.ppf(0.025), rv.ppf(0.975)
print(f"95% ETI: [{eti_low:.3f}, {eti_high:.3f}]")
print(f"Width: {eti_high - eti_low:.3f}")

Highest Density Interval (HDI)

The HDI (sometimes called HPD — Highest Posterior Density) is the shortest interval containing the specified probability mass. For a unimodal posterior, it is the interval $[L, U]$ such that every point inside has higher posterior density than every point outside:

For symmetric posteriors, ETI and HDI coincide. For skewed posteriors, the HDI is shorter and better captures where the parameter is most plausibly located. The ArviZ library provides az.hdi() which computes HDIs directly from MCMC samples.

When the Distinction Matters

Consider an exponentially distributed posterior (common for rate parameters near zero). The ETI includes a lower tail far from the bulk of the distribution, while the HDI correctly identifies that values near the mode are most credible. For reporting to practitioners, the HDI is usually the more honest summary.

import numpy as np
import arviz as az
from scipy.stats import expon

# Simulate samples from a skewed posterior (e.g., Gamma(2, rate=5))
rng = np.random.default_rng(7)
samples = rng.gamma(shape=2, scale=1/5, size=50000)

# Equal-tailed interval
eti = np.percentile(samples, [2.5, 97.5])

# HDI via ArviZ
hdi = az.hdi(samples, hdi_prob=0.95)

print(f"Posterior mean  : {samples.mean():.4f}")
print(f"Posterior median: {np.median(samples):.4f}")
print(f"95% ETI         : [{eti[0]:.4f}, {eti[1]:.4f}]  width={eti[1]-eti[0]:.4f}")
print(f"95% HDI         : [{hdi[0]:.4f}, {hdi[1]:.4f}]  width={hdi[1]-hdi[0]:.4f}")

Posterior Probability of Hypotheses

One of the most practically useful features of Bayesian inference is the ability to directly compute the probability that a hypothesis is true given the data. This requires no notion of $p$-values or null hypothesis significance testing.

For example, in an A/B test with conversion rates $\theta_A$ and $\theta_B$, you can compute $P(\theta_B > \theta_A \mid \mathcal{D})$ directly from posterior samples:

from scipy.stats import beta
import numpy as np

rng = np.random.default_rng(42)
# A: 45 conversions / 200 trials  => posterior Beta(46, 156)
# B: 55 conversions / 200 trials  => posterior Beta(56, 146)
theta_A = rng.beta(46, 156, size=100_000)
theta_B = rng.beta(56, 146, size=100_000)

prob_B_better = (theta_B > theta_A).mean()
lift = ((theta_B - theta_A) / theta_A).mean()
print(f"P(B > A | data)     : {prob_B_better:.3f}")
print(f"Expected lift of B  : {lift*100:.1f}%")

This output — a direct probability statement — is far more actionable for a product manager than "reject $H_0$ at $p = 0.04$."

The Credible Interval vs. Confidence Interval Revisited

The two constructions often yield numerically similar results under flat priors, which explains why practitioners often inadvertently apply the Bayesian interpretation to a frequentist interval. The conceptual difference is genuine:

• A 95% CI is a property of the procedure: if you ran the experiment 100 times and computed a CI each time, roughly 95 would contain the true parameter. A single CI either contains or does not contain the parameter — no probabilistic statement applies to it individually.
• A 95% CrI is a statement about the current parameter given current data: there is a 95% posterior probability that <!--MATHBLOCK15--> lies in this interval, conditional on the observed data and the prior.

For a practitioner making a decision based on a single experiment, the Bayesian credible interval is the more relevant object.

Why MCMC?

Conjugate priors yield beautiful closed-form posteriors, but they cover only a tiny fraction of real-world models. As soon as you add hierarchy, non-standard likelihoods, or correlated parameters, the normalizing constant $P(\mathcal{D}) = \int P(\mathcal{D} \mid \theta) P(\theta) d\theta$ becomes intractable. For a model with $d$ parameters, this integral has dimensionality $d$ — and even modest models can have dozens or hundreds of parameters.

Markov Chain Monte Carlo sidesteps the integral entirely. Instead of computing the posterior analytically, MCMC constructs a Markov chain whose stationary distribution is the posterior. After a warm-up (burn-in) period, samples from the chain are approximately independent draws from $P(\theta \mid \mathcal{D})$, which can be used to estimate posterior means, standard deviations, credible intervals, and any other posterior quantity.

Metropolis-Hastings Algorithm

The Metropolis-Hastings (MH) algorithm is the foundation of MCMC. Given a current state $\theta^{(t)}$:

1. Propose a new state <!--MATHBLOCK8--> from a symmetric proposal distribution <!--MATHBLOCK9--> (e.g., <!--MATHBLOCK10-->, <!--MATHBLOCK11-->).
2. Compute the acceptance ratio:

Note the marginal likelihood $P(\mathcal{D})$ cancels in the ratio — this is why MCMC only requires the unnormalized posterior.

1. Accept <!--MATHBLOCK13--> with probability <!--MATHBLOCK14-->; otherwise stay at <!--MATHBLOCK15-->.

The chain is guaranteed to converge to the posterior under mild regularity conditions (irreducibility and aperiodicity). The proposal scale $\sigma_p$ is a critical tuning parameter: too small and the chain explores slowly (high autocorrelation); too large and most proposals are rejected.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, beta

rng = np.random.default_rng(42)

# Target: Beta(5, 3) — we pretend we can only evaluate it up to a constant
def log_posterior(theta):
    if not (0 < theta < 1): return -np.inf
    return 4*np.log(theta) + 2*np.log(1 - theta)  # log Beta(5,3) kernel

n_samples = 20_000
proposal_std = 0.1
chain = np.empty(n_samples)
chain[0] = 0.5

accepted = 0
for t in range(1, n_samples):
    current = chain[t-1]
    proposal = current + rng.normal(0, proposal_std)
    log_alpha = log_posterior(proposal) - log_posterior(current)
    if np.log(rng.uniform()) < log_alpha:
        chain[t] = proposal
        accepted += 1
    else:
        chain[t] = current

print(f"Acceptance rate  : {accepted/(n_samples-1):.3f}  (target: ~0.23-0.44)")
print(f"Posterior mean   : {chain[2000:].mean():.4f}  (true: {5/8:.4f})")
print(f"Posterior std    : {chain[2000:].std():.4f}")

Hamiltonian Monte Carlo and NUTS

Random walk Metropolis-Hastings suffers from slow exploration in high-dimensional problems: to double effective sample size, you need roughly four times the computation. Hamiltonian Monte Carlo (HMC) uses gradient information to propose long-range moves that maintain high acceptance rates.

HMC augments the parameter space with auxiliary momentum variables $\mathbf{p}$ and simulates Hamiltonian dynamics:

The gradient $\nabla_\theta \log P(\theta \mid \mathcal{D})$ steers proposals toward high-probability regions. The No-U-Turn Sampler (NUTS) (Hoffman & Gelman 2014) adapts the trajectory length automatically, eliminating the difficult-to-tune leapfrog step count. NUTS is the default sampler in PyMC, Stan, and most modern probabilistic programming frameworks.

Convergence Diagnostics

MCMC samples are not i.i.d. — they are autocorrelated, and early samples depend on the starting point. Diagnostics assess whether chains have converged to the target distribution.

<!--MATHBLOCK19--> (Potential Scale Reduction Factor)

Run $m \geq 4$ chains from different starting points. $\hat{R}$ compares between-chain variance to within-chain variance:

where $W$ is the mean within-chain variance and $\hat{V}$ is an estimate of the marginal posterior variance. $\hat{R} \approx 1$ indicates convergence; values above 1.01 warrant concern; values above 1.1 indicate serious non-convergence.

Effective Sample Size (ESS)

Due to autocorrelation, $N$ MCMC samples contain less information than $N$ i.i.d. samples. The bulk ESS estimates the effective number of independent samples for central estimates, and the tail ESS assesses reliability of tail quantiles (used for credible interval bounds). ArviZ recommends $\text{ESS}_{\text{bulk}} \geq 100$ per chain and $\text{ESS}_{\text{tail}} \geq 100$ per chain as minimum thresholds.

Trace Plots and Divergences

Trace plots (parameter value vs. iteration) should look like stationary "fuzzy caterpillars" — all chains mixing rapidly across the same region. Funnel-shaped traces, slow drift, or chains stuck in local regions indicate problems. In HMC/NUTS, divergent transitions signal that the sampler encountered regions of high curvature and could not follow the Hamiltonian dynamics faithfully — they are strong indicators of model misspecification, requiring reparameterization or prior adjustment.

The Probabilistic Programming Paradigm

Probabilistic programming languages (PPLs) allow practitioners to specify generative models declaratively — in code that mirrors the mathematical model — and automatically derive posterior inference algorithms. PyMC (formerly PyMC3) is the leading Python PPL, using the Aesara/PyTensor computational backend for automatic differentiation (enabling HMC/NUTS) and providing a rich library of distributions, model components, and diagnostics through tight ArviZ integration.

The PyMC workflow follows a principled Bayesian pipeline:

1. Specify model: define priors and likelihood inside a pm.Model() context.
2. Prior predictive check: sample from the model before observing data to verify the prior is sensible.
3. Fit (sample): run NUTS to obtain posterior samples.
4. Diagnostics: check <!--MATHBLOCK0-->, ESS, trace plots, and divergences.
5. Posterior predictive check: compare data simulated from the posterior to the observed data.

Example: Bayesian Linear Regression

Consider a linear model $y_i = \alpha + \beta x_i + \epsilon_i$, $\epsilon_i \sim \mathcal{N}(0, \sigma^2)$. A Bayesian treatment places priors on all three parameters and computes their joint posterior.

import numpy as np
import pymc as pm
import arviz as az
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)
n = 100
true_alpha, true_beta, true_sigma = 1.5, 2.3, 0.8
x = rng.uniform(-2, 2, n)
y = true_alpha + true_beta * x + rng.normal(0, true_sigma, n)

with pm.Model() as linear_model:
    # Priors
    alpha = pm.Normal("alpha", mu=0, sigma=10)
    beta  = pm.Normal("beta",  mu=0, sigma=10)
    sigma = pm.HalfNormal("sigma", sigma=1)

    # Likelihood
    mu = alpha + beta * x
    likelihood = pm.Normal("y_obs", mu=mu, sigma=sigma, observed=y)

    # Prior predictive check
    prior_pred = pm.sample_prior_predictive(samples=200, random_seed=42)

    # Sample posterior
    trace = pm.sample(2000, tune=1000, chains=4,
                      target_accept=0.9, random_seed=42, progressbar=False)

    # Posterior predictive check
    post_pred = pm.sample_posterior_predictive(trace, random_seed=42)

# Diagnostics
print(az.summary(trace, var_names=["alpha", "beta", "sigma"]).to_string())
print(f"\nTrue values: alpha={true_alpha}, beta={true_beta}, sigma={true_sigma}")

The az.summary() output reports posterior mean, standard deviation, HDI bounds, $\hat{R}$, and ESS for each parameter. Well-specified models should show $\hat{R} < 1.01$ and ESS in the thousands.

Example: Hierarchical Model for Grouped Data

Hierarchical (multilevel) models are one of Bayesian inference's greatest practical advantages. Suppose you have test scores from students in $J = 8$ schools, with varying sample sizes per school. A pooled model ignores school-level variation; a fully unpooled model estimates each school independently and overfits small groups. A hierarchical model finds the optimal partial pooling automatically.

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

rng = np.random.default_rng(0)
J = 8
true_school_means = rng.normal(70, 10, J)
groups_sizes = rng.integers(5, 30, J)

school_idx = np.concatenate([np.full(n, j) for j, n in enumerate(groups_sizes)])
scores = np.concatenate([
    rng.normal(m, 5, n)
    for m, n in zip(true_school_means, groups_sizes)
])

with pm.Model() as hierarchical_model:
    # Hyperpriors
    mu_school    = pm.Normal("mu_school", mu=70, sigma=20)
    sigma_school = pm.HalfNormal("sigma_school", sigma=10)

    # School-level priors (non-centered parameterization for better geometry)
    z_school = pm.Normal("z_school", mu=0, sigma=1, shape=J)
    school_means = pm.Deterministic("school_means",
                                    mu_school + z_school * sigma_school)

    # Likelihood
    sigma_obs = pm.HalfNormal("sigma_obs", sigma=5)
    obs = pm.Normal("scores", mu=school_means[school_idx],
                    sigma=sigma_obs, observed=scores)

    trace_hier = pm.sample(2000, tune=1000, chains=4,
                           target_accept=0.9, random_seed=42, progressbar=False)

print(az.summary(trace_hier, var_names=["mu_school", "sigma_school",
                                         "school_means"]).to_string())

Non-Centered Parameterization

The non-centered parameterization (NCP) above — writing school<em>means = mu + z <em> sigma instead of directly sampling school</em>means ~ Normal(mu, sigma) — is a critical practical trick. In the centered parameterization, the school means and hyperparameters are highly correlated when group sample sizes are small, causing the sampler to struggle in a "funnel" geometry. The NCP decorrelates the geometry, enabling efficient sampling. PyMC's documentation flags this as the most important reparameterization in hierarchical modeling.

Posterior Predictive Checks

A posterior predictive check (PPC) assesses model fit by simulating new data from the posterior and comparing the distribution of simulated data to observed data. In PyMC, pm.sample<em>posterior</em>predictive() draws parameter samples from the posterior and then simulates observations.

Key diagnostics:

• Do simulated datasets match the spread and shape of observed data?
• Are there systematic biases (e.g., the model never generates extreme values that appear in the data)?
• Does the model reproduce key test statistics (mean, variance, skewness, proportion of zeros)?

PPCs are more informative than $p$-values because they directly show how* the model fails — whether it underpredicts tails, fails to capture multimodality, or generates impossible values.

Comparing Models the Bayesian Way

In frequentist statistics, model comparison often relies on information criteria (AIC, BIC) or likelihood ratio tests. Bayesian inference provides richer alternatives that naturally account for model complexity through the prior and propagate uncertainty across models.

Bayes Factors

The Bayes factor $\text{BF}_{12}$ is the ratio of marginal likelihoods (evidences) for two models $\mathcal{M}_1$ and $\mathcal{M}_2$:

The marginal likelihood integrates over all parameter values, automatically penalizing model complexity (Occam's razor): a model with more parameters that spreads its prior mass over a large space will have a lower marginal likelihood unless it fits the data substantially better. Jeffreys' scale provides rough guidelines: $\text{BF} > 3$ is moderate evidence, $\text{BF} > 10$ is strong, $\text{BF} > 100$ is decisive.

However, Bayes factors are sensitive to the prior specification on parameters — even diffuse priors can strongly influence them. Computing them requires the intractable normalizing constant, usually via bridge sampling or other specialized methods.

WAIC and LOO-CV

The Widely Applicable Information Criterion (WAIC) and Leave-One-Out Cross-Validation (LOO-CV) are practical, computationally tractable alternatives to Bayes factors for comparing models fitted to the same data.

LOO-CV estimates the expected out-of-sample predictive accuracy using Pareto-smoothed importance sampling (PSIS), which reuses the posterior samples from the full model without refitting. ArviZ's az.compare() directly reports LOO scores with standard errors:

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

rng = np.random.default_rng(0)
# True data-generating process: quadratic
x = rng.uniform(-3, 3, 80)
y = 0.5 * x**2 - x + 1 + rng.normal(0, 1.5, 80)

def fit_model(degree, x, y):
    X = np.column_stack([x**k for k in range(1, degree+1)])
    with pm.Model():
        intercept = pm.Normal("intercept", 0, 10)
        betas = pm.Normal("betas", 0, 5, shape=degree)
        sigma = pm.HalfNormal("sigma", 3)
        mu = intercept + pm.math.dot(X, betas)
        pm.Normal("y", mu=mu, sigma=sigma, observed=y)
        idata = pm.sample(1000, tune=500, chains=2,
                          idata_kwargs={"log_likelihood": True},
                          random_seed=42, progressbar=False)
    return idata

idata_linear = fit_model(1, x, y)
idata_quad   = fit_model(2, x, y)
idata_cubic  = fit_model(3, x, y)

comparison = az.compare(
    {"linear": idata_linear, "quadratic": idata_quad, "cubic": idata_cubic},
    ic="loo"
)
print(comparison[["elpd_loo", "p_loo", "d_loo", "weight"]].to_string())

The elpd<em>loo column reports the expected log predictive density. The model with the highest elpd</em>loo is preferred. d_loo shows the difference from the best model, and weight gives the stacking weights — the optimal convex combination of models for prediction.

Prior Sensitivity Analysis

In any serious Bayesian analysis, you should verify that your conclusions do not change substantially under reasonable alternative priors. This is called prior sensitivity analysis:

1. Fit the model with your chosen prior.
2. Fit it again with 2-3 alternative priors spanning a reasonable range.
3. Compare posterior summaries (means, credible intervals, probabilities of hypotheses).

If the conclusions are stable, you can report them with confidence. If the posterior changes dramatically, you need either more data or a more careful prior justification. Instability is not a failure — it is useful information about what the data alone cannot resolve.

Workflow Summary and Common Pitfalls

The Recommended Bayesian Workflow (Gelman et al.)

1. Choose model and priors with domain knowledge; consult literature.
2. Prior predictive check: simulate data from the prior. Does it cover the plausible observation space? Is it over-concentrated?
3. Fit to simulated data (fake-data simulation): verify the sampler recovers the known true parameters before fitting real data.
4. Fit to real data and run diagnostics (<!--MATHBLOCK7-->, ESS, trace plots, divergences).
5. Posterior predictive check: does the model generate data similar to observed data?
6. Prior sensitivity analysis: do conclusions change under alternative reasonable priors?
7. Model comparison: if multiple candidate models, use LOO-CV.

Pitfalls to Avoid

• Centered parameterizations in hierarchies: Use non-centered parameterizations (see previous section) to avoid funnel geometry.
• Ignoring divergent transitions: Divergences in NUTS are not statistical noise — they indicate the sampler failed to explore correctly. Fix the model (often via reparameterization or more informative priors) before trusting any posterior summaries.
• Confusing prior and posterior predictive checks: Prior predictive checks run before fitting and check prior reasonableness. Posterior predictive checks run after fitting and check model fit. Both are essential.
• Multiple comparison issues in Bayesian A/B tests: Even Bayesian tests can be gamed by sequential peeking if you are using posterior probability as a decision threshold without a proper decision-theoretic framework. Use expected loss or a proper stopping rule.
• Treating MAP as a full Bayesian inference: MAP estimation discards all uncertainty quantification. It is a useful point estimate for optimization but is not a substitute for full posterior inference.