Marginal Likelihood

The marginal likelihood is the probability of the data under a model, averaged over every parameter value the prior allows:

\[ P(D) = \int P(D \mid \theta)\,P(\theta)\,d\theta \]

It is the denominator of Bayes' theorem, and it has two names for its two jobs: a normalizing constant when you want a posterior, and the model evidence when you want to compare models.

The constant, and what it does

import numpy as np
from scipy.special import comb, beta as beta_fn
from scipy.stats import beta as beta_dist
from scipy.integrate import quad

n, k = 10, 7
print(f"  {k} heads in {n} tosses, with a uniform Beta(1,1) prior on p")
exact = comb(n, k) * beta_fn(k + 1, n - k + 1)
numeric, _ = quad(lambda p: comb(n, k) * p ** k * (1 - p) ** (n - k), 0, 1)
print(f"    P(D) = C({n},{k}) * B({k + 1}, {n - k + 1}) = {exact:.8f}")
print(f"    numerical integration agrees: {numeric:.8f}")
print(f"    with a uniform prior it is exactly 1/(n+1) = {1 / (n + 1):.8f} for every k")
"""The marginal likelihood averages the likelihood over the prior. Under a flat
   prior every possible count is equally likely in advance, hence 1/(n+1)."""

unnorm = lambda p: p ** k * (1 - p) ** (n - k)
area, _ = quad(unnorm, 0, 1)
print(f"    the unnormalized posterior integrates to {area:.8f}")
print(f"    dividing by it gives {quad(lambda p: unnorm(p) / area, 0, 1)[0]:.8f}")
print(f"    which is the Beta({k + 1}, {n - k + 1}) density: "
      f"{quad(lambda p: beta_dist(k + 1, n - k + 1).pdf(p), 0, 1)[0]:.8f}")

def evidence(a, b):
    return float(comb(n, k) * beta_fn(k + a, n - k + b) / beta_fn(a, b))

print(f"{'prior':>21}{'P(D)':>14}{'posterior mean':>18}")
for name, (a, b) in (("Beta(1, 1) flat", (1, 1)), ("Beta(2, 2) weak", (2, 2)),
                     ("Beta(20, 20) firm", (20, 20)), ("Beta(1, 20) opposed", (1, 20))):
    A, B = k + a, n - k + b
    print(f"{name:>21}{evidence(a, b):>14.8f}{A / (A + B):>18.6f}")
span = evidence(20, 20) / evidence(1, 20)
print(f"  a factor of {span:,.0f}. A prior you would tolerate for estimation can be")

fair = comb(n, k) * 0.5 ** n
for name, (a, b) in (("Beta(1, 1)", (1, 1)), ("Beta(20, 20)", (20, 20))):
    bf = evidence(a, b) / fair
    print(f"    evidence under {name:>13} / evidence under p = 0.5 exactly = {bf:.4f}")
print(f"  the same data reads as {evidence(1, 1) / fair:.4f} under one prior and"
      f" {evidence(20, 20) / fair:.4f} under another -")
#   7 heads in 10 tosses, with a uniform Beta(1,1) prior on p
#     P(D) = C(10,7) * B(8, 4) = 0.09090909
#     numerical integration agrees: 0.09090909
#     with a uniform prior it is exactly 1/(n+1) = 0.09090909 for every k
#     the unnormalized posterior integrates to 0.00075758
#     dividing by it gives 1.00000000
#     which is the Beta(8, 4) density: 1.00000000
#                 prior          P(D)    posterior mean
#       Beta(1, 1) flat    0.09090909          0.666667
#       Beta(2, 2) weak    0.11188811          0.642857
#     Beta(20, 20) firm    0.12327012          0.540000
#   Beta(1, 20) opposed    0.00005126          0.258065
#   a factor of 2,405. A prior you would tolerate for estimation can be
#     evidence under    Beta(1, 1) / evidence under p = 0.5 exactly = 0.7758
#     evidence under  Beta(20, 20) / evidence under p = 0.5 exactly = 1.0519
#   the same data reads as 0.7758 under one prior and 1.0519 under another -

For 7 heads in 10 tosses under a uniform prior, \(P(D) = 0.09090909\) — computed in closed form from the Beta function and confirmed by numerical integration. Under a flat prior it is exactly \(1/(n+1)\) for every possible count, which is the arithmetic version of saying the prior considered all counts equally likely in advance.

The second block shows the normalizing job. The unnormalized posterior integrates to 0.00075758; dividing by that gives a density integrating to 1, which is the Beta(8, 4) you would have written down directly.

This is why the constant is usually ignored. In conjugate problems it has a closed form, and in sampling methods it cancels: MCMC only ever needs ratios of the unnormalized posterior, which is precisely what lets it work when the integral is intractable.

As evidence, it is not the posterior

The third block is the part that catches people. Four priors give posterior means within 0.41 of each other — a difference of 41 percentage points, substantial for a probability — while the marginal likelihood spans a factor of 2,405, from 0.12327012 to 0.00005126.

The reason is structural. A posterior is dominated by the likelihood once there is enough data, so a vague prior washes out. The marginal likelihood is an average over the prior and never stops depending on it: spreading the prior across regions the data rules out lowers \(P(D)\) no matter how much data arrives. A prior you would happily use for estimation can be ruinous for evidence.

That is not a bug so much as the mechanism behind its main use. Because a model that spreads its prior over many possibilities pays for the ones the data rejects, the marginal likelihood penalizes unnecessary flexibility automatically — an Occam factor, with no complexity term written by hand. The same property makes it sensitive to priors chosen for convenience.

Bayes factors carry the same dependence

A ratio of two marginal likelihoods is a Bayes factor, and the last block prices the consequence. Against a model that fixes \(p = 0.5\), the flexible model scores 0.7758 under a Beta(1,1) prior and 1.0519 under a Beta(20,20) one — one below 1 and one above, so this comparison changes direction with the prior. The wider table also shows substantial estimation sensitivity.

Three practices follow. Report the prior alongside any Bayes factor, since the number is uninterpretable without it. Run the comparison under a materially different prior and report both. And avoid improper priors entirely for model comparison — they leave \(P(D)\) defined only up to an arbitrary constant, so the Bayes factor is arbitrary too, even where the posterior is perfectly well behaved.

Finally, the integral is hard. Outside conjugate cases it has to be estimated, and estimators of the marginal likelihood are substantially harder to get right than posterior samples — which is why cross-validation and information criteria are often used for model comparison instead.

References

• SciPy documentation, scipy.special.beta and scipy.integrate.quad. docs.scipy.org