Marginal Likelihood

Marginal Likelihood (Evidence)

Given a model with prior $p(\theta \mid \mathbf{\varphi })$ and likelihood $p(y\mid \theta )$, the marginal likelihood (also called evidence, prior predictive density, or model evidence) is the prior-weighted average of the likelihood,
$p(y\mid \mathbf{\varphi })=\int p(y\mid \theta )p(\theta \mid \mathbf{\varphi })d\theta .$
It is the normalizing constant in Bayes’ rule: $p(\theta \mid y,\mathbf{\varphi })=p(y\mid \theta )p(\theta \mid \mathbf{\varphi })/p(y\mid \mathbf{\varphi })$. It is the probability of the observed data under the model before any $\theta$ is committed to. Three roles:

1. Hyperparameter estimation (empirical Bayes). Maximizing $p(y\mid \mathbf{\varphi })$ over $\mathbf{\varphi }$ tunes the prior to the data. The resulting Type-II MLE $\hat{\mathbf{\varphi }}$ is used as a plug-in for unit-level inference. This is the standard way to choose Dirichlet concentrations, Beta-Binomial pseudo-counts, Gaussian-prior variances, or Gaussian-process kernel lengthscales in practice.
2. Model comparison. The ratio of evidences for two models is the Bayes factor ${\mathrm{BF}}_{10}=p(y\mid M_1)/p(y\mid M_0)$, which Bayesians use to quantify how strongly the data favor one model over another.
3. Posterior normalization. It is the denominator of Bayes’ rule. In conjugate models it is closed form; in non-conjugate models it is computed by Laplace approximation, bridge sampling, or nested sampling.

Closed-form evidence in conjugate families. When $p(\theta )$ and $p(y\mid \theta )$ form a conjugate pair, the evidence has a clean ratio-of-normalizing-constants form. For the Beta-Binomial: $p(y\mid \alpha ,\beta )=\left(\genfrac{}{}{0ex}{}{n}{y}\right)B(\alpha +y,\beta +n-y)/B(\alpha ,\beta )$. For the Gamma-Poisson: $p(y\mid \alpha ,\beta )={\beta }^{\alpha }\Gamma (\alpha +\sum y_i)/[\Gamma (\alpha )(\beta +n{)}^{\alpha +\sum y_i}\prod y_i!]$. For the Dirichlet-Multinomial: $p(\mathbf{y}\mid \mathbf{\alpha })=\left(\genfrac{}{}{0ex}{}{n}{\mathbf{y}}\right)B(\mathbf{\alpha }+\mathbf{y})/B(\mathbf{\alpha })$.

Laplace approximation. When closed form is unavailable,
$p(y\mid \mathbf{\varphi })\approx (2\pi {)}^{d/2}|-{\nabla }^2\log p(\hat{\theta },y\mid \mathbf{\varphi }){|}^{-1/2}p(\hat{\theta }\mid \mathbf{\varphi })p(y\mid \hat{\theta }),$
where $\hat{\theta }$ is the MAP and the Hessian is evaluated there. See 3.9 Laplace Approximation.

Asymptotic approximation — BIC. For large $n$ with a proper prior,
$-2\log p(y\mid \mathbf{\varphi })\approx -2\log p(y\mid \hat{\theta })+d\log n+O(1),$
which is the Bayesian Information Criterion. BIC approximates $-2\log$ evidence to within an $O(1)$ error, giving a fast model-comparison surrogate that trades off fit against dimension.

Prerequisites

• 3.1 Bayesian Framework — Bayes’ rule, prior, posterior, evidence.
• 3.3 Bayes Factors — ratio of evidences used for model comparison.
• 3.9 Laplace Approximation — standard tool for approximating evidence when closed form is unavailable.

Intuition

The marginal likelihood measures how well a model predicts the observed data before any parameters are fit. It is the average likelihood across all parameter values the prior considers plausible. A model that concentrates its prior on values producing the observed data receives high evidence; a model that spreads its prior over regions of parameter space far from the data’s implications receives low evidence, because most of the prior mass is “wasted.”

This produces an automatic Occam-razor effect that no frequentist model-selection criterion gets for free. A flexible model with many parameters can fit almost any dataset, but its prior must spread mass across the whole parameter space — so the prior density at any particular $\theta$ is small. A simpler model has less flexibility but concentrates its prior predictions, so prior density near the truth is larger. The marginal likelihood weighs fit against how sharply the prior had to be tuned to produce that fit, and rewards the model that used its flexibility efficiently. This is the formal content of the 3.3 Bayes Factors “built-in Occam” property.

When the prior has hyperparameters $\mathbf{\varphi }$, viewing $p(y\mid \mathbf{\varphi })$ as a function of $\mathbf{\varphi }$ yields a standard likelihood function for $\mathbf{\varphi }$ alone — the “marginal” in marginal likelihood is over $\theta$. Maximizing this likelihood picks the prior that best explains the observed data, which is the empirical-Bayes recipe. In conjugate models the maximization is analytic; in non-conjugate models it reduces to numerical optimization of a (usually log-concave) function.

Properties

(Prior Predictive) Evidence is literally the prior-averaged likelihood.

$p(y\mid \mathbf{\varphi })=\int p(y\mid \theta )p(\theta \mid \mathbf{\varphi })d\theta ={\mathbb{E}}_{\theta \sim p(\cdot \mid \mathbf{\varphi })}\left[p(y\mid \theta )\right].$
This is the probability of $y$ before the data is used, with uncertainty about $\theta$ integrated out per the prior. It is not a function of $\theta$ — that has been integrated away. It is a function of $\mathbf{\varphi }$ (and the observed $y$). The term “prior predictive density” is sometimes preferred for this reason.

(Closed Form in Conjugate Families) Each conjugate pair has an evidence expression in terms of ratios of normalizing constants.

For the Beta-Binomial:
$p(y\mid \alpha ,\beta )=\left(\genfrac{}{}{0ex}{}{n}{y}\right)\frac{B(\alpha +y,\beta +n-y)}{B(\alpha ,\beta )},B(\alpha ,\beta )=\frac{\Gamma (\alpha )\Gamma (\beta )}{\Gamma (\alpha +\beta )}.$
For the Gamma-Poisson with i.i.d. counts and unit exposure:
$p(\mathbf{y}\mid \alpha ,\beta )=\frac{{\beta }^{\alpha }}{\Gamma (\alpha )}\cdot \frac{\Gamma (\alpha +\sum y_i)}{(\beta +n{)}^{\alpha +\sum y_i}}\cdot {\prod }_i\frac{1}{y_i!}.$
For the Normal-Normal with known ${\sigma }^2$:
$p(\mathbf{y}\mid {\mu }_0,{\tau }_0^2)=(2\pi {\sigma }^2{)}^{-n/2}\sqrt{\frac{{\tau }_n^2}{{\tau }_0^2}}\exp \left[-\frac{1}{2{\sigma }^2}{\sum }_i{y}_i^2+\frac{{\mu }_n^2}{2{\tau }_n^2}-\frac{{\mu }_0^2}{2{\tau }_0^2}\right],$
with ${\tau }_n^2$ and ${\mu }_n$ the posterior parameters. These closed forms allow Type-II MLE over $\mathbf{\varphi }$ by simple root-finding.

(Empirical-Bayes Estimator) Maximizing evidence over $\mathbf{\varphi }$ gives the Type-II MLE.

$\hat{\mathbf{\varphi }}=\arg {\max }_{\mathbf{\varphi }}p(y\mid \mathbf{\varphi })=\arg {\max }_{\mathbf{\varphi }}\int p(y\mid \theta )p(\theta \mid \mathbf{\varphi })d\theta .$
In the Beta-Binomial hierarchy, the first-order condition involves the digamma function $\psi (x)=d/dx\log \Gamma (x)$: setting $\partial \log p/\partial \alpha =0$ and $\partial \log p/\partial \beta =0$ gives two digamma equations solved numerically by Newton’s method. The Type-II MLE is then used as a plug-in prior; see 3.6 Hierarchical Models.

(Built-in Occam Penalty via BIC) Evidence penalizes model dimension logarithmically via the Laplace approximation.

For a proper prior on $\theta \in {\mathbb{R}}^d$ with log-likelihood smooth and $n$ growing, the Laplace approximation gives
$\log p(y\mid \mathbf{\varphi })\approx \log p(y\mid \hat{\theta })-\frac{d}{2}\log n+O(1).$
Hence $-2\log p(y\mid \mathbf{\varphi })\approx -2\log p(y\mid \hat{\theta })+d\log n+O(1)$, recovering the BIC. Evidence thus performs model selection that rewards fit and penalizes dimension — an asymptotic surrogate for the full Bayes factor that needs only the MLE and the dimension of $\theta$. See 3.9 Laplace Approximation and 3.3 Bayes Factors.

(Log-Evidence Decomposes as Prequential Prediction) Total evidence equals the sum of one-step-ahead predictive log-likelihoods.

By the chain rule,
$\log p(y_1,\dots ,y_n\mid \mathbf{\varphi })={\sum }_{i=1}^n\log p(y_i\mid y_{1:i-1},\mathbf{\varphi }),$
where each term is the posterior-predictive likelihood conditional on all earlier data. This prequential or leave-future-out decomposition is a Bayesian version of cross-validation: the evidence is the cumulative predictive accuracy along the data stream. Two models with comparable overall evidence but different prequential profiles differ in when they predict well — useful diagnostic information.

(Computation in Non-Conjugate Models) Without conjugacy, evidence is computed by Laplace, thermodynamic integration, bridge sampling, or nested sampling.

• Laplace: $p(y)\approx (2\pi {)}^{d/2}|H{|}^{-1/2}p(\hat{\theta })p(y\mid \hat{\theta })$; fast and accurate when posterior is approximately Gaussian. See 3.9 Laplace Approximation.
• Bridge sampling / thermodynamic integration: interpolate between prior and posterior and integrate; robust for multimodal posteriors.
• Nested sampling (Skilling 2006): evaluate $p(y\mid \theta )$ at decreasing prior-mass levels; standard in astrophysical inference.
• Harmonic-mean estimator: $\hat{p}(y)=1/{\mathbb{E}}_{\theta \mid y}[1/p(y\mid \theta )]$; infamous for infinite variance and should not be used without stabilization.

(Evidence Depends on the Prior, Not Just the Likelihood) Unlike the posterior for a well-identified parameter, evidence is sensitive to prior choice.

As the prior width ${\tau }_0\to \infty$, evidence $p(y\mid \mathbf{\varphi })\to 0$ at rate ${\tau }_0^{-d}$ (the prior spreads mass ever more thinly). Two models with identical likelihoods but different prior widths therefore have different evidences, and the Bayes factor inherits this sensitivity — the source of the Jeffreys-Lindley paradox in 3.3 Bayes Factors. Weakly informative priors, not improper ones, are mandatory for Bayes-factor analyses.

(Marginal Likelihood Under Improper Priors Is Undefined) Model comparison via Bayes factors requires proper priors.

If $p(\theta )\propto 1$ is improper, the evidence $\int p(y\mid \theta )d\theta$ is finite only if the likelihood is integrable in $\theta$, and the result has no absolute scale. The Bayes factor ${\mathrm{BF}}_{10}$ then inherits an arbitrary scaling constant. Posterior inference within a single model is immune to this (the normalizing constant cancels), but cross-model comparison is not. See the 3.3 Bayes Factors Jeffreys-Lindley discussion for what to do when only improper priors feel principled.

Examples

Example 1: Beta-Binomial Evidence and Hyperparameter Estimation

A single binomial observation $y\mid \theta \sim \text{Binomial}(n,\theta )$ with a $\text{Beta}(\alpha ,\beta )$ prior has evidence
$p(y\mid \alpha ,\beta )=\left(\genfrac{}{}{0ex}{}{n}{y}\right)\frac{B(\alpha +y,\beta +n-y)}{B(\alpha ,\beta )}.$
For $n=10$, $y=3$, and the non-informative uniform prior $\alpha =\beta =1$:
$p(y\mid 1,1)=\left(\genfrac{}{}{0ex}{}{10}{3}\right)\frac{B(4,8)}{B(1,1)}=120\cdot \frac{\Gamma (4)\Gamma (8)/\Gamma (12)}{1}=120\cdot \frac{6\cdot 5040}{39916800}=120\cdot \frac{1}{1320}=\frac{1}{11}.$
This is exactly what the uniform prior on $\theta$ predicts: if $\theta$ is uniform, then $y$ is uniform on $\{0,1,\dots ,10\}$, and $P(y=3)=1/11$.

Empirical Bayes over a hierarchy: With $J$ groups $y_j\mid n_j,{\theta }_j$, ${\theta }_j\sim \text{Beta}(\alpha ,\beta )$, the hierarchical evidence is $p(\mathbf{y}\mid \alpha ,\beta )={\prod }_j\text{BetaBin}(y_j;n_j,\alpha ,\beta )$. Its log-derivatives involve digamma functions:
$\frac{\partial \log p}{\partial \alpha }={\sum }_j[\psi (\alpha +y_j)-\psi (\alpha +\beta +n_j)-\psi (\alpha )+\psi (\alpha +\beta )],$
and similarly for $\beta$. Setting both to zero and solving numerically gives $\hat{\alpha },\hat{\beta }$; Newton iteration converges in a handful of steps.

Example 2: Evidence for Model Comparison

Compare $M_0:y\sim N(0,1)$ (point null) against $M_1:y\sim N(\mu ,1)$ with $\mu \sim N(0,{\tau }^2)$ under $n$ i.i.d. observations with sample mean $\bar{y}$. The evidences are
$p(y\mid M_0)=(2\pi {)}^{-n/2}\exp \left(-\frac{1}{2}{\sum }_i{y}_i^2\right),$
$p(y\mid M_1)=(2\pi {)}^{-n/2}(1+n{\tau }^2{)}^{-1/2}\exp \left(-\frac{1}{2}{\sum }_i{y}_i^2+\frac{n^2{\bar{y}}^2{\tau }^2}{2(1+n{\tau }^2)}\right).$
Their ratio is the Bayes factor
${\mathrm{BF}}_{10}=(1+n{\tau }^2{)}^{-1/2}\exp \left(\frac{n^2{\bar{y}}^2{\tau }^2}{2(1+n{\tau }^2)}\right),$
which exhibits the Jeffreys-Lindley paradox: for fixed $z=\sqrt{n}\bar{y}$ (frequentist significance), ${\mathrm{BF}}_{10}\to 0$ as $\tau \to \infty$, i.e. a vague prior on $\mu$ guarantees evidence against any non-null hypothesis regardless of the data’s frequentist strength.

Example 3: BIC as a Fast Proxy for Log-Evidence

A Bayesian linear regression with $d$ predictors on $n$ observations has log-likelihood at the MLE $\ell (\hat{\beta })=-\frac{n}{2}\log (2\pi {\hat{\sigma }}^2)-\frac{n}{2}$. BIC is
$\mathrm{BIC}=-2\ell (\hat{\beta })+d\log n=n\log ({\hat{\sigma }}^2)+n\log (2\pi )+n+d\log n.$
For two nested models with $d_1>d_0$ predictors, the approximate log-evidence difference is
$\log {\mathrm{BF}}_{10}\approx -\frac{1}{2}({\mathrm{BIC}}_1-{\mathrm{BIC}}_0)=\ell ({\hat{\beta }}_1)-\ell ({\hat{\beta }}_0)-\frac{d_1-d_0}{2}\log n.$
The $\log n$ penalty for each additional parameter means that as $n$ grows, BIC (and asymptotic evidence) ever-more-strongly prefers the smaller model unless the larger model’s likelihood improvement also grows with $n$. This is stricter than AIC ($2$ per parameter, no $\log n$), and matches the evidence’s own Occam penalty.

Derivation

Goal: derive the Laplace approximation for evidence.

The posterior is proportional to $g(\theta ):=p(\theta \mid \mathbf{\varphi })p(y\mid \theta )$, and the evidence is $p(y\mid \mathbf{\varphi })=\int g(\theta )d\theta$. Write $h(\theta ):=\log g(\theta )$ and let $\hat{\theta }$ be its maximizer (the MAP). Taylor-expanding $h$ around $\hat{\theta }$ to second order:
$h(\theta )\approx h(\hat{\theta })+\frac{1}{2}(\theta -\hat{\theta }{)}^{\top }H(\hat{\theta })(\theta -\hat{\theta }),$
where $H(\hat{\theta })={\nabla }^{2}h(\hat{\theta })$ is the Hessian (negative definite at a maximum). Exponentiating and integrating:
$p(y\mid \mathbf{\varphi })\approx e^{h(\hat{\theta })}\int \exp \left[-\frac{1}{2}(\theta -\hat{\theta }{)}^{\top }(-H(\hat{\theta }))(\theta -\hat{\theta })\right]d\theta =e^{h(\hat{\theta })}(2\pi {)}^{d/2}|-H(\hat{\theta }){|}^{-1/2}.$
Since $e^{h(\hat{\theta })}=p(\hat{\theta }\mid \mathbf{\varphi })p(y\mid \hat{\theta })$,
$p(y\mid \mathbf{\varphi })\approx (2\pi {)}^{d/2}|-H(\hat{\theta }){|}^{-1/2}p(\hat{\theta }\mid \mathbf{\varphi })p(y\mid \hat{\theta }).$

Asymptotic BIC derivation. For $n$ i.i.d. observations, $H(\hat{\theta })=n\bar{H}$ where $\bar{H}$ is order 1. Taking logs:
$\log p(y\mid \mathbf{\varphi })\approx \log p(y\mid \hat{\theta })+\log p(\hat{\theta }\mid \mathbf{\varphi })+\frac{d}{2}\log (2\pi )-\frac{d}{2}\log n-\frac{1}{2}\log |\bar{H}|.$
All $\mathbf{\varphi }$- and $\theta$-independent terms are absorbed into a constant; the prior term $\log p(\hat{\theta }\mid \mathbf{\varphi })=O(1)$ does not grow with $n$. The dominant terms are $\log p(y\mid \hat{\theta })$ and $-\frac{d}{2}\log n$, giving BIC.

See Also

• 3.1 Bayesian Framework — evidence as the Bayes-rule normalizing constant.
• 3.3 Bayes Factors — evidence ratio for model comparison.
• 3.9 Laplace Approximation — standard approximation when evidence has no closed form.
• 3.5 Conjugate Families — conjugate pairs admit closed-form evidence.
• 3.6 Hierarchical Models — empirical Bayes plugs in the evidence-maximizing hyperparameters.
• 3.11 MCMC — full-Bayes alternatives to plugging in the Type-II MLE.