Bayesian Model Comparison & BMA

Topic 25 §25.5 introduced two distinct quantities that both begin with the word “marginal”:

• Marginal likelihood $m(\mathbf{y}) = \int L(\theta)\pi(\theta)\,d\theta$: integrates the data likelihood against the prior. A single scalar summarizing the fit of the model.
• Posterior predictive $p(\tilde y \mid \mathbf{y}) = \int p(\tilde y \mid \theta) p(\theta \mid \mathbf{y})\,d\theta$: integrates a new-data likelihood against the posterior. A full density over new observations.

The marginal likelihood is the evidence supporting the model; the posterior predictive is the model’s prediction under that evidence. §27.2 formalizes the first, §27.7 develops the second.

MCMC works because it cancels $m(\mathbf{y})$. That cancellation is exactly what makes $m(\mathbf{y})$ hard to estimate from MCMC draws — the draws carry information about the posterior’s shape but not its scale. Every technique in §27.6 is a way to recover the lost scale: Laplace uses the analytic Gaussian-integral formula, bridge and path sampling tie posterior draws to a calibrated proposal, nested sampling transforms the integral into a 1-D quadrature over prior-mass levels. Each has failure modes; the one-size-fits-all estimator does not exist.

Given a candidate set $\mathcal{M}$, prior model probabilities $\\{p(\mathcal{M}_k)\\}_{k=1}^M$, and marginal likelihoods $\\{m_k(\mathbf{y})\\}_{k=1}^M$, the posterior probability of model $k$ is

$p(\mathcal{M}_k \mid \mathbf{y}) \;=\; \frac{p(\mathcal{M}_k)\, m_k(\mathbf{y})}{\sum_{j=1}^M p(\mathcal{M}_j)\, m_j(\mathbf{y})}.$

The denominator is the total evidence $p(\mathbf{y}) = \sum_j p(\mathcal{M}_j) m_j(\mathbf{y})$ — a number independent of $k$, acting purely as a normalizer.

Consider $y \mid \theta \sim \text{Binomial}(n, \theta)$ with $\theta \sim \text{Beta}(\alpha, \beta)$. The marginal likelihood is

$m(y) \;=\; \int_0^1 \binom{n}{y}\, \theta^y (1-\theta)^{n-y} \cdot \frac{\theta^{\alpha-1}(1-\theta)^{\beta-1}}{B(\alpha, \beta)}\, d\theta.$

The integrand is proportional to $\text{Beta}(\alpha + y, \beta + n - y)$, so the integral evaluates to

$m(y) \;=\; \binom{n}{y} \cdot \frac{B(\alpha + y, \beta + n - y)}{B(\alpha, \beta)}.$

Taking logs: $\log m(y) = \log \binom{n}{y} + \log B(\alpha + y, \beta + n - y) - \log B(\alpha, \beta)$. For $(n, y, \alpha, \beta) = (20, 12, 1, 1)$ the answer is $\log(1/21) \approx -3.04452$ — the closed-form reference we use in §27.6 to verify the bridge / path / Laplace estimators (brief §6.2 T27.4). The function betaBinomialLogMarginal in bayes.ts computes this with stable lnGamma arithmetic.

Bernardo–Smith (1994) classify model-comparison settings by the status of truth. M-closed: the true data-generating process is in $\mathcal{M}$. M-open: the truth is outside $\mathcal{M}$; we’re choosing the best approximating model. M-complete: truth lies in a (possibly unenumerated) extension of $\mathcal{M}$, but we treat $\mathcal{M}$ as a working set. BMA is principled in M-closed and M-complete (it averages over the posterior within $\mathcal{M}$); for M-open problems, stacking (formalml) or predictive-weight optimization is often preferred because BMA weights concentrate on a single “best” model that may still be far from truth. Most applied settings are M-open; the distinction matters when interpreting BMA weights.

A common pitfall: calling $m(\mathbf{y})$ the “likelihood of the model.” It is not. The likelihood of the model at $\theta$ is $L(\theta) = p(\mathbf{y} \mid \theta)$ — a function of $\theta$. The marginal likelihood $m(\mathbf{y})$ is a scalar — the likelihood $L$ averaged against the prior. Two models with identical MLE fits (same peak of $L$) can have very different $m$ if one has a concentrated prior that matches the data and the other has a diffuse prior that allocates probability mass to regions where $L$ is low.

$p(\mathcal{M}_k \mid \mathbf{y})$ depends on the entire set $\mathcal{M}$ through the denominator. Adding or removing candidate models — even ones with negligible weight — can shift the reported posterior probabilities. This is not a bug; it’s a statement about what Bayesian model comparison actually asks. The question “which model explains the data best?” is inherently relative: best compared to what? A model with $p(\mathcal{M}_1 \mid \mathbf{y}) = 0.95$ in a pair and $p(\mathcal{M}_1 \mid \mathbf{y}) = 0.03$ in an eight-model set is not contradictory — it answers different questions.

Independence and conditional independence continue to be written with the Topic 16 §16.9 convention $\perp\!\!\!\perp$. The M-open / M-closed / M-complete trichotomy is used narratively (§27.2 Rem 3) without a dedicated symbol. When $\mathcal{M}_k$ and $\mathcal{M}_{k'}$ share parameters (nested models), we’ll say so explicitly — the shared parameter is integrated against the joint prior, not multiplied.

For any two models $\mathcal{M}_0$ and $\mathcal{M}_1$ with positive prior probabilities,

$\underbrace{\frac{p(\mathcal{M}_1 \mid \mathbf{y})}{p(\mathcal{M}_0 \mid \mathbf{y})}}_{\text{posterior odds}} \;=\; \underbrace{\frac{p(\mathcal{M}_1)}{p(\mathcal{M}_0)}}_{\text{prior odds}} \;\times\; \underbrace{\mathrm{BF}_{10}(\mathbf{y})}_{\text{data's update}}.$

The Bayes factor is thus prior-free in one sense — it’s the part of the comparison due to the data, independent of prior model probabilities — and prior-dependent in another: each $m_k$ depends on the within-model prior $\pi_k$. Rem 7 below elaborates the second sense; Rem 12–14 follow the consequences.

Let $y_1, \ldots, y_n \mid \theta \sim \mathcal{N}(\theta, 1)$ with $H_0: \theta = 0$ and $H_1: \theta \sim \mathcal{N}(0, \tau^2)$. Then $\bar y \mid H_0 \sim \mathcal{N}(0, 1/n)$ and $\bar y \mid H_1 \sim \mathcal{N}(0, 1/n + \tau^2)$ (by the convolution of likelihood and prior). The Bayes factor is

$\mathrm{BF}_{10} \;=\; \frac{\phi(\bar y;\, 0,\, 1/n + \tau^2)}{\phi(\bar y;\, 0,\, 1/n)} \;=\; \frac{1}{\sqrt{1 + n\tau^2}}\, \exp\!\left(\frac{\bar y^2}{2}\left[\frac{1}{1/n} - \frac{1}{1/n + \tau^2}\right]\right).$

This is the Lindley setup — the featured case of §27.5. For $\bar y = 0.3, n = 100, \tau = 1$: $z = \bar y \sqrt n = 3$, $r = n\tau^2 = 100$, and $\mathrm{BF}_{10} \approx 8.57$ (moderate evidence for $\mathcal{M}_1$; lindleyBayesFactor(3, 100, 1) in bayes.ts). Same $z$ with $\tau = 100$: $\mathrm{BF}_{10} \approx 0.09$ (strong evidence for $\mathcal{M}_0$). The paradox lives between these two lines.

The Bayes factor is prior-free over model probabilities but prior-dependent over parameter distributions. The $\pi_k$ in $m_k(\mathbf{y}) = \int L_k(\theta) \pi_k(\theta)\,d\theta$ is a hard input: a diffuse $\pi_k$ under $\mathcal{M}_1$ diluting probability mass over a large parameter volume yields a smaller $m_1$ than a concentrated $\pi_k$ matching the data. This is not a nuisance — it’s the Occam factor (§27.6 Thm 4) acting as a model-complexity penalty. But it does mean you must think about your within-model priors when you use Bayes factors; vague “non-informative” priors are actively informative here.

The scale above is Jeffreys’ original proposal and remains common in applied work, but nothing in Bayesian probability requires you to use it — $\mathrm{BF}_{10} = 5$ has a specific probabilistic meaning (multiplies the prior odds by 5) regardless of whether a convention calls it “substantial.” Treat the labels as rough communication tools; report the numerical $\mathrm{BF}$ and let the reader apply their own threshold. Kass & Raftery 1995 §3.2 proposes a slightly coarser alternative scale (thresholds at 3, 20, 150); the distinction rarely matters in practice.

Under $H_0: \theta = \theta_0$ and $H_1: \theta \neq \theta_0$ with prior $\pi_1$ on the alternative,

$\mathrm{BF}_{10} \;=\; \int \frac{L(\theta)}{L(\theta_0)}\, \pi_1(\theta)\, d\theta \;=\; \mathbb{E}_{\pi_1}\!\left[\Lambda(\mathbf{y}; \theta, \theta_0)\right],$

where $\Lambda$ is the Topic 18 likelihood ratio. The BF is thus the prior-averaged LRT — Neyman–Pearson picks a specific alternative, Bayes averages over the prior. This is why §17.10 Rem 15 and §18 Rem 25 point here: Bayes factors are the evidence-level analog of the LRT decision rule.

Let $Y_1, \ldots, Y_n \overset{\mathrm{iid}}{\sim} \mathcal{N}(\theta, \sigma^2)$ with $\sigma^2$ known. Consider $H_0: \theta = 0$ against $H_1: \theta \neq 0$ with prior $\theta \mid H_1 \sim \mathcal{N}(0, \tau^2)$. Let $z = \bar Y \sqrt{n}/\sigma$ be the observed z-statistic. Then

$\mathrm{BF}_{10}(z) \;=\; \frac{1}{\sqrt{1 + r}}\, \exp\!\left(\frac{z^2}{2} \cdot \frac{r}{1+r}\right), \qquad r := \frac{n\tau^2}{\sigma^2}.$

In particular, for any fixed $z$ and $n$, $\mathrm{BF}_{10}(z) \to 0$ as $\tau \to \infty$ — the frequentist rejects $H_0$ (a function of $z$ alone) while the Bayesian overwhelmingly favors $H_0$.

Three values of $\tau$, same $z = 3$ (p-value $\approx 0.0027$, strong frequentist rejection):

• $\tau = 1$: $r = 100$. $\mathrm{BF}_{10} = (1/\sqrt{101}) \exp(4.5 \cdot 100/101) \approx 0.0995 \cdot 86.1 \approx 8.57$.
• $\tau = 10$: $r = 10{,}000$. $\mathrm{BF}_{10} = (1/\sqrt{10001}) \exp(4.5 \cdot 10000/10001) \approx 0.01 \cdot 90.0 \approx 0.90$.
• $\tau = 100$: $r = 10^6$. $\mathrm{BF}_{10} \approx 0.001 \cdot 90.0 \approx 0.0900$.

Same data, same likelihood, same model — the prior-scale parameter alone flips the conclusion from “moderate evidence for $\mathcal{M}_1$” to “strong evidence for $\mathcal{M}_0$.” This is not a numerical accident; it is the diffuse-prior penalty in action. The LindleyParadoxExplorer component below traces the full $\mathrm{BF}_{10}$ curve across $\tau$ at adjustable $(z, n)$.

Lindley’s paradox is specific to the setting where $H_0$ assigns a point mass to a single parameter value (here $\theta = 0$). If instead $H_0: |\theta| < \epsilon$ (a small interval) — an “interval null” — the paradox dissolves: both $m_0$ and $m_1$ scale similarly as $\tau$ grows, so their ratio stabilizes. Applied Bayesian practice often recommends interval nulls for this reason (cf. Berger & Sellke 1987). Topic 27 keeps the point-null framing because that’s what sharpens the contrast with frequentist hypothesis testing, but the pragmatic advice is: if your null is approximate, make the interval null explicit.

A related phenomenon: fix $\bar y$ (not $z$) and let $n \to \infty$. Then $z = \bar y \sqrt n \to \infty$ too, and one might expect $\mathrm{BF}_{10} \to \infty$ always. Not so — for large $\tau$, the $1/\sqrt{1+r}$ prefactor can dominate the exponential. Bartlett (1957, “A comment on D. V. Lindley’s statistical paradox”) notes this: the paradox arises because the prior-scale $\tau$ and sample-scale $\sqrt n$ interact non-trivially. Keeping $\tau$ fixed as $n$ grows is itself a modeling choice — often the wrong one.

The paradox doesn’t reveal a flaw in Bayesian inference; it reveals a demand for care in prior specification. A diffuse $\tau \gg 1$ under $H_1$ says: “before seeing data, I believed $\theta$ could plausibly be anywhere in $[-3\tau, 3\tau]$ — a range of hundreds of standard deviations.” Under that prior, observing $z = 3$ is modest evidence for $\mathcal{M}_1$ relative to the prior commitment, and $\mathcal{M}_0$‘s concentrated mass at 0 wins the comparison. The resolution: calibrate $\tau$ to what you actually believe about the effect size. For a standardized coefficient in a well-designed study, $\tau \sim 1$ is realistic. For a “diffuse” effect in an exploratory analysis, $\tau \sim 10$ might be appropriate — but don’t go to $\tau = 100$ unless your substantive knowledge genuinely says so. This is the advice of applied Bayesian workflow: prior elicitation matters precisely because Bayes factors take priors seriously.

Let $\mathcal{M}$ be a regular parametric model with $k$ parameters and prior $\pi$ continuous and positive at the posterior mode $\hat\theta$. Under the regularity conditions of Topic 14 §14.7 (existence of $\hat\theta_{\mathrm{MLE}}$; positive-definite observed information $\mathbf{J}_n = n \mathbf{I}(\hat\theta) + O_p(\sqrt n)$),

$-2 \log m(\mathbf{y}) \;=\; -2 \log L(\hat\theta_{\mathrm{MLE}}) \,+\, k \log n \,+\, O_p(1).$

The right-hand side is BIC; the $O_p(1)$ remainder absorbs the Occam factor (prior volume, Hessian determinant, Gaussian-integral constants) that doesn’t scale with $n$.

For $(n, y, \alpha, \beta) = (20, 12, 1, 1)$, the unnormalized log-posterior is $\log L + \log \pi = 12 \log\theta + 8 \log(1-\theta)$ (ignoring the binomial coefficient). The mode is $\hat\theta = 12/20 = 0.6$, and the second derivative at the mode is $-12/\hat\theta^2 - 8/(1-\hat\theta)^2 = -12/0.36 - 8/0.16 = -33.33 - 50 = -83.33$. So $H = 83.33$ and the Laplace estimate (after adding the binomial normalization) is

$\log m \approx \log \binom{20}{12} + 12 \log 0.6 + 8 \log 0.4 + \tfrac{1}{2}\log(2\pi) - \tfrac{1}{2}\log 83.33 \approx -3.06,$

close to the closed-form $-3.0445$ (Ex 1). The Laplace approximation is accurate to $O(1/n)$ under regularity — impressive given its closed-form simplicity — but deteriorates for small $n$ or multimodal posteriors. marginalLikelihoodLaplace in bayes.ts handles the general multivariate case via Cholesky log-determinant.

Let $q_1, q_2$ be unnormalized densities on $\Theta$ with normalizing constants $Z_i = \int q_i(\theta)\, d\theta$ and normalized densities $p_i = q_i/Z_i$. For any bridge function $h: \Theta \to \mathbb{R}$ such that both expectations below exist and are nonzero,

$\frac{Z_1}{Z_2} \;=\; \frac{\mathbb{E}_{p_2}[\,q_1(\theta)\, h(\theta)\,]}{\mathbb{E}_{p_1}[\,q_2(\theta)\, h(\theta)\,]}.$

Run bridge sampling on $(n, y, \alpha, \beta) = (20, 12, 1, 1)$ with 8000 posterior draws from Beta(13, 9) and 8000 proposal draws from a deliberately-offset Beta(12.7, 9.1). The iterative fixed point converges in 3 iterations to $\log \hat m \approx -3.0446$, matching the closed-form $-3.04452$ (Ex 1) to $\sim 10^{-4}$. Over 30 independent MC seeds the standard deviation is $\approx 7 \times 10^{-4}$ — two orders of magnitude tighter than naive importance sampling on the same sample budget. The BridgeSamplingConvergence component below traces the iteration on three preset problems.

Let $p_\beta(\theta) \propto L(\theta)^\beta \pi(\theta)$ be the power posterior at inverse temperature $\beta \in [0, 1]$. Then

$\log m(\mathbf{y}) \;=\; \int_0^1 \mathbb{E}_{p_\beta}\!\left[\log L(\theta)\right]\, d\beta.$

The estimator replaces the integral with trapezoidal quadrature over a discrete $\beta$-grid $\\{0 = \beta_0 < \beta_1 < \cdots < \beta_K = 1\\}$, and each expectation with an MCMC average over draws from $p_{\beta_j}$.

Same inputs as Ex 5. Discretize $\beta$ on $\\{0, 0.1, 0.2, \ldots, 1\\}$ (11 points), sample 2000 draws per $\beta$ via MCMC targeting $L(\theta)^\beta \pi(\theta)$, compute $\widehat{\mathbb{E}_{p_\beta}[\log L]}$ as the sample mean at each temperature, and trapezoidally integrate. Result: $\log \hat m \approx -3.05$, within $\sim 5 \times 10^{-3}$ of the closed form — slightly looser than bridge at the same total sample budget because path sampling spends compute on off-target temperatures. Where path sampling wins: the integrand $\mathbb{E}_{p_\beta}[\log L]$ is smooth in $\beta$, so even a coarse grid gives a decent answer; where it loses: choosing the $\beta$-schedule is a tuning decision with no universal rule (equidistant is fine for quick-and-dirty, geometric $\beta_j \propto 2^{j/K}$ is better for skewed posteriors).

Define the prior-mass function $X(\lambda) := \int_{\\{\theta : L(\theta) > \lambda\\}} \pi(\theta)\,d\theta$ — the prior volume at likelihood level $\lambda$ or higher. $X$ is monotone decreasing in $\lambda$ with $X(0) = 1$ and $X(\infty) = 0$. Then

$m(\mathbf{y}) \;=\; \int_0^\infty X(\lambda)\, d\lambda \;=\; \int_0^1 \lambda(X)\, dX,$

where $\lambda(X) = X^{-1}$ is the likelihood at prior-mass level $X$.

Skilling’s estimator simulates iid live-point draws from the prior, removes the lowest-likelihood point, replaces it with a constrained-prior draw at higher likelihood, and records the removed point’s likelihood. After $K$ iterations with $N$ live points, the dead-point likelihoods $L_i$ are associated with expected prior-mass shrinkage $X_i = \exp(-i/N)$, and

$\hat m(\mathbf{y}) \;\approx\; \sum_{i=1}^K L_i\, (X_{i-1} - X_i) \,+\, \text{live-point remainder}.$

The four estimators above all assume a continuous prior $\pi$. For sparsity-inducing priors — horseshoe (Carvalho–Polson–Scott 2010), spike-and-slab (Mitchell–Beauchamp 1988), Bayesian lasso — the prior has a sharp peak or mixture structure that Laplace handles poorly and bridge sampling requires careful proposal tuning for. Path sampling and nested sampling cope better but need fine $\beta$-grids or dense live-point arrangements near the origin. §27.10 Rem 26 points to formalml for the full development.

Newton & Raftery’s (1994) harmonic-mean estimator $\hat m = 1 / \mathbb{E}_{p(\theta \mid \mathbf{y})}[1/L(\theta)]$ converts posterior draws into a marginal-likelihood estimate with no proposal — just $1/L(\theta_s)$ averaged over posterior draws $\theta_s$. Elegant and seductive, but almost always has infinite variance: the integrand $1/L(\theta)$ has a heavy right tail under the posterior (since the posterior is proportional to $L\pi$, not $1/L$). Radford Neal famously called it “the worst Monte Carlo method ever proposed.” harmonicMeanEstimate in bayes.ts is kept for pedagogical contrast — use it to demonstrate the pathology, never in production.

• Bridge sampling is usually the default — stable, fast, handles non-Gaussian posteriors well, works for any dimension where MCMC works.
• Path sampling is the right choice when you already have power-posterior draws (e.g., from parallel tempering) — the marginal-likelihood estimate is a cheap byproduct.
• Nested sampling wins for strongly multimodal posteriors or when the marginal likelihood is the primary target (astrophysics model comparison is where nested sampling originated).
• Laplace is the starting point: a 30-second sanity check you should run before investing in any MCMC-based estimator. If Laplace and bridge disagree by orders of magnitude, one of them is wrong and you need to figure out which.

An even sharper variant of Lindley’s paradox, due to Stone (1982), Dawid, Stone & Zidek (1973): with improper (non-integrable) priors, the marginal likelihood is undefined, and any Bayes factor between models with improper priors is arbitrary up to a constant. The consequence: never use improper priors for model comparison. Proper priors with sensible scale (Rem 12) are the robust default. Topic 27 doesn’t belabor this case because the brief prescribed scope ends at proper priors; formalml’s “Objective Bayes” topic develops the full story.

Four downstream constructions consume $\hat m(\mathbf{y})$: (i) posterior model probabilities via §27.2 Thm 1, (ii) Bayes factors via §27.4 Thm 2, (iii) BMA weights via §27.7, (iv) evidence-level hypothesis testing as in §27.5’s Lindley paradox. All four treat $\hat m$ as a deterministic input, but MC variance in $\hat m$ propagates to each — bridge’s $\sim 10^{-3}$ log-scale SE on Beta-Binomial becomes $\sim 10^{-3}$ SE in the log-BF, which at BF $= 8.57$ translates to $\pm 0.01$ in BF space. Good enough for applied work; worth propagating via bootstrap for high-stakes comparisons.

Under the log-score loss $\mathcal{L}(\tilde y, q) = -\log q(\tilde y)$, the optimal predictive distribution — the one minimizing expected loss under the joint posterior over $(\mathcal{M}_k, \theta_k)$ — is the BMA mixture

$q^*(\tilde y) \;=\; \sum_{k=1}^M w_k^{\mathrm{BMA}}\, p(\tilde y \mid \mathcal{M}_k, \mathbf{y}), \qquad w_k^{\mathrm{BMA}} = p(\mathcal{M}_k \mid \mathbf{y}).$

Any non-BMA predictive — hard model selection, stacking, equal-weight averaging — incurs strictly greater expected log-loss against the posterior-truth. Proof: standard Bayesian-decision-theory argument; see HOE1999 §2.

Synthetic data: $y_i = \sin(2\pi x_i / 10) + 0.3 x_i + \varepsilon_i$, $\varepsilon_i \sim \mathcal{N}(0, 0.25^2)$, $x_i \sim U(0, 10)$, $n = 20$. Fit three polynomial models — degree 1 (linear), degree 2 (quadratic), degree 3 (cubic) — and compare via BIC-approximated marginal likelihoods. For this specific sample, BIC strongly prefers degree 3 (posterior model probability $\approx 1.00$, weights $(0, 0, 1)$), so BMA collapses to hard selection. The BMAPredictiveComparison component lets you flip to uniform weights $(1/3, 1/3, 1/3)$ or single-model weights and watch the predictive-fan change — uniform BMA produces a noticeably wider predictive band in the mid-range where the three models disagree most. The lesson: BMA’s predictive uncertainty reflects model-set uncertainty; hard selection under-reports it.

Thm 8 says BMA is the log-loss optimal predictive — it does not say BMA’s posterior model probabilities are calibrated measures of “which model is true.” A prior over models times a marginal likelihood gives a posterior over models; calling one of those probabilities “the probability that $\mathcal{M}_1$ is true” requires M-closed. In applied work, treat BMA weights as predictive weights (sum to 1, used for averaging) rather than truth probabilities (calibrated measures of model correctness). The distinction matters when you report “the probability that the effect is real is 0.95” — that claim is BMA’s weight, not a frequentist coverage guarantee.

Applied BMA over $M \gg 10$ candidate models is expensive and often fruitless — most models have negligible posterior probability. Madigan & Raftery (1994) propose Occam’s window: restrict the candidate set to models whose posterior probability is within a factor $O$ (typically 20 or 100) of the best model. This loses strict decision-theoretic optimality but is tractable. Most applied BMA software implements Occam’s window by default; bmaPredictive in bayes.ts takes pre-computed weights and assumes the pruning (if any) has already happened.

Stacking (Breiman 1996, and Yao, Vehtari, Simpson & Gelman 2018 for the Bayesian version) directly optimizes predictive weights on held-out data rather than deriving them from marginal likelihoods. Under M-closed, stacking and BMA converge asymptotically; under M-open, stacking typically wins (it’s optimizing the thing you care about — predictive performance — rather than assuming the true model is enumerated). The BMAPredictiveComparison component sticks to BMA because Topic 27 is about model comparison under Bayesian principles; formalml’s stacking topic develops the M-open alternative.

Back to the canonical tiny problem: $n = 5$ Normal observations, $\sigma = 1$, posterior $\mu \mid \mathbf{y} \sim \mathcal{N}(\bar y, \sigma^2/n) = \mathcal{N}(0.04, 0.2)$ with 2000 posterior draws. Compute the deviance $D(\mu) = -2 \sum \log \phi(y_i \mid \mu, 1)$ at each draw, take the posterior mean $\bar D$, compare to $D(\bar\mu)$ at the posterior-mean $\bar\mu = 0.04$. The effective parameter count $p_{\mathrm{DIC}} = \bar D - D(\bar\mu) \approx 1$ (as expected for a one-parameter model under weak prior regularity), and $\mathrm{DIC} = D(\bar\mu) + 2 p_{\mathrm{DIC}} \approx 12.06$. The dic function in bayes.ts computes this in one line per input.

Fit a simple linear regression $y \mid \beta \sim \mathcal{N}(X\beta, \sigma^2)$ with flat prior and compute both WAIC and PSIS-LOO on the same posterior draws. For a correctly-specified well-conditioned problem both land within 0.1% of each other; for a poorly-specified or outlier-heavy regression, PSIS-LOO can be $\sim 1$ larger (deviance scale) with high Pareto-$\hat k$ diagnostics on the outliers. The waic and psisLoo functions in bayes.ts both accept the same pointwise log-likelihood matrix and return compatible outputs — a direct drop-in comparison with no re-sampling required.

PSIS-LOO’s reliability hinges on the Pareto shape parameter $\hat k$ per observation. Vehtari’s rule of thumb: $\hat k < 0.5$ means PSIS is reliable; $0.5 \le \hat k < 0.7$ is borderline; $\hat k \ge 0.7$ flags the observation as an outlier whose LOO estimate is unreliable — usually a sign of model misspecification on that point. The psisLoo function returns both the aggregate LOO value and a per-observation $\hat k$ vector plus a nProblematic count. Operationally: always check the $\hat k$ distribution before reporting LOO values.

DIC has well-documented pathologies: $p_{\mathrm{DIC}}$ can be negative (spurious), doesn’t handle mixture models or hierarchical models well, and depends on the choice of posterior summary ($\bar\mu$ vs posterior median, say). Spiegelhalter himself later (Spiegelhalter et al. 2014) recommended treating DIC as a legacy tool — useful for historical comparison with BUGS-era literature, not a default for new analyses. Modern Bayesian practice defaults to PSIS-LOO or WAIC.

Watanabe (2010) proved $\mathrm{WAIC} = \mathrm{LOO} + O_p(1/n^2)$ under singular learning theory — an asymptotic equivalence that holds under substantially weaker regularity than DIC/AIC/BIC (no requirement of Fisher-information regularity). The proof requires machinery outside Topic 27’s scope; the result is stated-and-cited here. Practically: WAIC and LOO agree to 3–4 significant figures on most well-behaved problems.