Bayesian Model Comparison

In machine learning, we constantly face a fundamental question: Which model should we use? Given a dataset, we might consider a linear model, a polynomial model, a neural network, or countless other possibilities. Each model makes different assumptions about the underlying data-generating process, and each has different capacity to fit complex patterns.

The classical approach to model selection—choosing based on training error—is deeply flawed. A model with more parameters can always achieve lower training error, but this often leads to overfitting. Cross-validation offers some remedy, but it requires holding out data and doesn't provide a principled probabilistic framework.

Bayesian inference offers an elegant solution: the marginal likelihood. Also called the model evidence, the marginal likelihood quantifies how well a model explains the observed data while automatically penalizing model complexity. It's the foundation upon which all Bayesian model comparison rests.

To understand marginal likelihood, we must first revisit the structure of Bayesian inference. Given a model $\mathcal{M}$ with parameters $\boldsymbol{\theta}$ and observed data $\mathcal{D}$, Bayes' theorem gives us the posterior distribution over parameters:

$$p(\boldsymbol{\theta} | \mathcal{D}, \mathcal{M}) = \frac{p(\mathcal{D} | \boldsymbol{\theta}, \mathcal{M}) \cdot p(\boldsymbol{\theta} | \mathcal{M})}{p(\mathcal{D} | \mathcal{M})}$$

The denominator $p(\mathcal{D} | \mathcal{M})$ is the marginal likelihood. It's called 'marginal' because it's obtained by marginalizing (integrating) over all possible parameter values:

$$p(\mathcal{D} | \mathcal{M}) = \int p(\mathcal{D} | \boldsymbol{\theta}, \mathcal{M}) \cdot p(\boldsymbol{\theta} | \mathcal{M}) , d\boldsymbol{\theta}$$

Let's parse this integral carefully:

• $p(\mathcal{D} | \boldsymbol{\theta}, \mathcal{M})$: The likelihood function. Given specific parameter values $\boldsymbol{\theta}$, how probable is the observed data?

• $p(\boldsymbol{\theta} | \mathcal{M})$: The prior distribution. Before seeing data, how plausible are different parameter values under our model assumptions?

• The integral: We average the likelihood over all possible parameter values, weighted by their prior probability.

The key insight: The marginal likelihood doesn't evaluate the model at any single point in parameter space. Instead, it asks: 'Considering all possible parameter values this model could have, weighted by their prior probability, how well does this model explain the data on average?'

The marginal likelihood has a beautiful geometric interpretation that illuminates why it naturally penalizes complexity.

Imagine the parameter space as a landscape. The height at each point $\boldsymbol{\theta}$ represents the likelihood $p(\mathcal{D} | \boldsymbol{\theta})$—how well those particular parameters explain the data. The prior $p(\boldsymbol{\theta})$ acts as a 'sampling weight' over this landscape.

The marginal likelihood is the volume under this weighted likelihood surface. It's the integral of the likelihood, weighted by the prior, over the entire parameter space.

Now consider two models:

Model A: Simple model (few parameters)

• Parameter space is low-dimensional
• Prior is concentrated (each parameter region gets more prior mass)
• If the model fits well, the likelihood peak is high
• Volume = limited extent × high peak × concentrated prior = potentially high

Model B: Complex model (many parameters)

• Parameter space is high-dimensional
• Prior is spread thin (prior mass distributed over vast space)
• Maximum likelihood may be higher (more flexibility)
• But volume = vast extent × high peak × diffuse prior = often lower

An analogy from probability: Consider a Gaussian distribution. A narrow Gaussian (small variance) concentrates its probability mass in a small region, achieving high density there. A wide Gaussian (large variance) spreads its mass thin, achieving low density everywhere. The marginal likelihood follows similar logic—a simpler model that 'commits' to a smaller region of data space can achieve higher evidence if its predictions are correct.

The key principle: Marginal likelihood rewards models that make specific, accurate predictions. A complex model that could explain almost any data pattern gets penalized because it hasn't 'committed' to this particular pattern—its prior probability is spread across all possibilities.

The complexity-penalizing behavior of marginal likelihood can be made explicit through the concept of the Occam factor.

For models where the posterior is approximately Gaussian (which is common due to Bernstein-von Mises theorem), we can decompose the log marginal likelihood as:

$$\log p(\mathcal{D} | \mathcal{M}) \approx \log p(\mathcal{D} | \hat{\boldsymbol{\theta}}) + \log \text{Occam Factor}$$

where $\hat{\boldsymbol{\theta}}$ is the maximum a posteriori (MAP) estimate.

The Occam factor captures the ratio of posterior volume to prior volume in parameter space:

$$\text{Occam Factor} = \frac{\text{Posterior Volume}}{\text{Prior Volume}} = \frac{\Delta\boldsymbol{\theta}{\text{posterior}}}{\Delta\boldsymbol{\theta}{\text{prior}}}$$

Since data constrains parameters, the posterior is always more concentrated than the prior. Thus the Occam factor is always less than 1, contributing a negative term to the log marginal likelihood.

For a d-dimensional Gaussian posterior approximation, the Occam factor can be computed explicitly:

$$\log \text{Occam Factor} = \log \frac{(2\pi)^{d/2} |\boldsymbol{\Sigma}_{\text{post}}|^{1/2}}{\int p(\boldsymbol{\theta}) d\boldsymbol{\theta}}$$

where $\boldsymbol{\Sigma}_{\text{post}}$ is the posterior covariance matrix.

For a model with isotropic priors of width $\sigma_{\text{prior}}$ and isotropic posterior of width $\sigma_{\text{post}}$:

$$\log \text{Occam Factor} \approx d \log \frac{\sigma_{\text{post}}}{\sigma_{\text{prior}}}$$

This is negative (a penalty) because $\sigma_{\text{post}} < \sigma_{\text{prior}}$. The penalty grows linearly with the number of parameters $d$, but is modulated by how much each parameter is constrained by the data.

For certain model classes, the marginal likelihood can be computed exactly in closed form. These cases typically involve conjugate prior relationships, where the prior and likelihood combine to produce a posterior in the same family as the prior.

Example 1: Bayesian Linear Regression

Consider linear regression with Gaussian noise: $$y = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}, \quad \boldsymbol{\epsilon} \sim \mathcal{N}(\mathbf{0}, \sigma^2\mathbf{I})$$

With a Gaussian prior on weights: $$\boldsymbol{\beta} \sim \mathcal{N}(\mathbf{0}, \tau^2\mathbf{I})$$

The marginal likelihood (assuming known $\sigma^2, \tau^2$) is:

$$p(\mathbf{y} | \mathbf{X}, \sigma^2, \tau^2) = \mathcal{N}(\mathbf{y} | \mathbf{0}, \tau^2\mathbf{X}\mathbf{X}^\top + \sigma^2\mathbf{I})$$

This is a multivariate Gaussian with mean zero and covariance $\mathbf{K} = \tau^2\mathbf{X}\mathbf{X}^\top + \sigma^2\mathbf{I}$.

Example 2: Beta-Bernoulli Model

For coin flipping with a Beta prior on the probability $\theta$:

• Prior: $\theta \sim \text{Beta}(\alpha, \beta)$
• Likelihood: $k$ successes in $n$ trials: $p(k | \theta) = \binom{n}{k}\theta^k(1-\theta)^{n-k}$

The marginal likelihood is:

$$p(k | \alpha, \beta) = \binom{n}{k} \cdot \frac{B(\alpha + k, \beta + n - k)}{B(\alpha, \beta)}$$

where $B(\cdot, \cdot)$ is the Beta function.

This closed-form result arises from the conjugacy of the Beta prior with the Binomial likelihood.

When closed-form computation is impossible, we need approximation methods. Each has distinct strengths and is appropriate for different situations.

5.1 Laplace Approximation

The Laplace approximation treats the posterior as approximately Gaussian, centered at its mode (MAP estimate). This gives:

$$\log p(\mathcal{D} | \mathcal{M}) \approx \log p(\mathcal{D} | \hat{\boldsymbol{\theta}}) + \log p(\hat{\boldsymbol{\theta}}) + \frac{d}{2}\log(2\pi) - \frac{1}{2}\log|\mathbf{H}|$$

where $\hat{\boldsymbol{\theta}}$ is the MAP estimate and $\mathbf{H}$ is the Hessian of the negative log posterior at the mode.

Intuition: We approximate the integrand with a Gaussian bump centered at the MAP, then compute the volume under this Gaussian exactly.

5.2 Bayesian Information Criterion (BIC)

BIC is a rough approximation to the log marginal likelihood:

$$\text{BIC} = \log p(\mathcal{D} | \hat{\boldsymbol{\theta}}) - \frac{d}{2}\log n$$

where $d$ is the number of parameters and $n$ is the number of data points.

Derivation sketch: BIC emerges from the Laplace approximation when we assume:

1. The prior is flat (uniform) over a reasonable region
2. The data provides $O(n)$ information, so posterior variance scales as $O(1/n)$

Interpretation: The $\frac{d}{2}\log n$ term is the complexity penalty. It grows with the number of parameters but is modulated by sample size—more data justifies more parameters.

5.3 Monte Carlo Integration

For arbitrary posteriors, we can estimate the marginal likelihood using Monte Carlo sampling:

Simple Monte Carlo (often poor): $$p(\mathcal{D}) \approx \frac{1}{S}\sum_{s=1}^{S} p(\mathcal{D} | \boldsymbol{\theta}^{(s)}), \quad \boldsymbol{\theta}^{(s)} \sim p(\boldsymbol{\theta})$$

This samples from the prior and averages likelihoods. It's unbiased but has extremely high variance—prior samples rarely land near high-likelihood regions.

Importance Sampling (better): $$p(\mathcal{D}) \approx \frac{1}{S}\sum_{s=1}^{S} \frac{p(\mathcal{D} | \boldsymbol{\theta}^{(s)}) p(\boldsymbol{\theta}^{(s)})}{q(\boldsymbol{\theta}^{(s)})}, \quad \boldsymbol{\theta}^{(s)} \sim q(\boldsymbol{\theta})$$

We sample from a proposal distribution $q$ concentrated near high-posterior regions.

Bridge Sampling, Path Sampling, Thermodynamic Integration: Advanced methods that construct a 'bridge' between prior and posterior, often via a sequence of intermediate distributions.

A critical but often overlooked aspect of marginal likelihood is its sensitivity to prior specification. Unlike posterior inference, where reasonable priors often lead to similar conclusions when data is abundant, the marginal likelihood can be highly sensitive to prior choices.

Why Priors Matter More for Model Comparison

Consider two priors for a parameter $\theta$:

• Prior A: $\theta \sim \mathcal{N}(0, 1)$
• Prior B: $\theta \sim \mathcal{N}(0, 100)$

If the true parameter is $\theta = 0.5$ and data strongly supports this:

• For posterior inference: Both priors lead to posteriors concentrated near 0.5. The prior gets 'washed out' by data.

• For marginal likelihood: Prior B spreads probability mass over a huge range. Since most of that range has low likelihood, the marginal likelihood is lower. Prior A 'predicted' the true parameter region better, so it gets higher evidence.

The implication: When comparing models via marginal likelihood, you're also comparing how well each prior 'predicted' what the data would show. Vague, diffuse priors get penalized.

Best Practices for Prior Specification

1. Use informative priors when possible: If you have domain knowledge about reasonable parameter ranges, encode it. This isn't 'cheating'—it's using available information.

2. Conduct sensitivity analysis: Compute marginal likelihoods under several reasonable prior specifications. If conclusions are robust, you can be more confident.

3. Use proper priors: Improper priors (e.g., flat priors over infinite ranges) lead to undefined or zero marginal likelihoods. Always use proper probability distributions.

4. Consider objective Bayes approaches: Methods like reference priors or Jeffreys priors attempt to be 'least informative' while remaining proper, but have their own limitations.

5. Match prior scales across models: When comparing models, ensure that priors on shared parameters have compatible scales. A model isn't 'worse' just because its prior was poorly calibrated.

Having understood the theory, let's address practical challenges in using marginal likelihood for model selection.

Numerical Stability

Marginal likelihoods span enormous ranges. For typical datasets, log marginal likelihoods can be -1000 or lower. Always work in log space:

✗ BAD:  ml_ratio = ml_model1 / ml_model2  # Underflow/overflow
✓ GOOD: log_ml_diff = log_ml_model1 - log_ml_model2

Nested Model Comparison

When comparing nested models (e.g., linear vs. polynomial regression), marginal likelihood automatically handles the 'nesting' through the Occam factor. You don't need to account for it separately—it emerges from the integral.

Non-Nested Model Comparison

Marginal likelihood excels here. Unlike likelihood ratio tests (which require nested models), marginal likelihood compares any two models on equal footing. A Gaussian mixture vs. a kernel density estimator? Both get evaluated on how well they explain the data, weighted by their prior commitments.

Let's consolidate what we've learned about marginal likelihood:

Looking ahead: Marginal likelihood is rarely used in isolation. In the next page, we'll learn how to combine marginal likelihoods from different models to compute Bayes factors—the ratio of evidences that quantifies how much more one model is supported by the data compared to another.

Bayes factors provide a principled, calibrated scale for interpreting the strength of evidence, moving beyond mere p-values to a fully probabilistic comparison framework.