Bayesian Inference and Parameter Estimation

Every measurement is an inverse problem. Nature runs the forward map — pick a parameter θ, and the apparatus produces data x with a probability the physics fixes. The experimenter holds the data and wants the parameter, so the map must be run backwards. Bayesian inference is the exact rule for that inversion, and its claim is radical in its simplicity: treat the unknown θ as a random variable, not because nature gambles but because your knowledge is incomplete, and let a single equation carry belief from before the experiment to after it. The machinery has three moving parts. A prior π(θ) encodes what you knew before looking. A likelihood 𝓛(x∣θ) — the forward model read as a function of θ — grades how well each candidate explains the data. Their product, renormalised, is the posterior p(θ∣x). The renormalising constant, the evidence Z, does double duty: it makes the posterior a genuine distribution, and its value is the figure of merit for comparing whole models. Everything an experiment reports flows from the posterior: a best estimate is a posterior average, an error bar a posterior width, a prediction integrates the parameters out, a model verdict is a ratio of evidences. One line of algebra, an ocean of consequence — from Laplace weighing Saturn in 1810 to LIGO weighing black holes today, the same theorem does the learning.

Key takeaways
  • Bayesian inference treats an unknown parameter θ as a random variable and updates a prior belief π(θ) into a posterior p(θ∣x) using observed data x through Bayes' theorem.
  • The engine is p(θ∣x) = 𝓛(x∣θ)π(θ)/Z, where 𝓛 is the likelihood and Z = ∫𝓛(x∣θ)π(θ)dθ is the marginal likelihood, or evidence.
  • Conjugate priors (Beta–Binomial, Gamma–Poisson, Normal–Normal) give closed-form posteriors and expose the prior as pseudo-counts of imagined prior data.
  • Point estimates are posterior functionals: the mean minimises squared loss, the median absolute loss, the mode (MAP) 0–1 loss; credible intervals read straight off the posterior.
  • As data accumulate the posterior concentrates into a Gaussian about the maximum-likelihood estimate with covariance the inverse Fisher information (Bernstein–von Mises).
  • Model comparison uses the Bayes factor Z₁/Z₂, and the evidence's integral over parameter space builds in an automatic Occam penalty for unused complexity.
  • Key idea: Bayes' theorem is the rule for learning — posterior ∝ likelihood × prior, with the proportionality fixed by integrating over all parameter values; the parameter is random because knowledge is incomplete, and the posterior is the complete summary of that knowledge after the data.
The core insight

Bayes' theorem is the rule for learning: posterior ∝ likelihood × prior, with the proportionality fixed by integrating over all parameter values. The parameter is treated as random not because nature is random but because our knowledge is incomplete, and the posterior is the complete summary of that knowledge after the data.

Where this lives today

Bayesian parameter estimation is the reporting language of precision physics. In cosmology, the Planck collaboration's headline numbers — the matter density Ω_m, the Hubble constant H₀, the six parameters of ΛCDM — are posterior means with credible intervals, extracted from the microwave-background power spectrum by Markov-chain Monte Carlo with foreground and calibration nuisances marginalised away; the "Hubble tension" is the statement that the CMB posterior and the local distance-ladder posterior overlap at less than the 5σ level. In gravitational-wave astronomy, every LIGO–Virgo–KAGRA event is announced with posteriors on masses, spins, distance, and sky location from samplers comparing strain against waveform templates through a Gaussian noise likelihood, and the detection significance is a Bayes factor of signal against noise. Combining particle-physics measurements with correlated systematics is the multivariate Normal–Normal update. Beyond the frontier the same machinery drives nested sampling, variational inference, and the probabilistic-programming toolkits now standard in data analysis — all numerical answers to Laplace's one integral, the evidence Z.

The physical picture

Picture the parameter axis laid flat, and over it a landscape of plausibility — the prior — high where you expected θ to lie. The data arrive as a second shape, the likelihood, tall over the values that would have produced what you saw and near zero elsewhere. Bayesian learning is nothing more exotic than multiplying these two shapes point by point and rescaling to unit area: where both are high the posterior is high, where either vetoes it collapses. The first figure shows the update in one frame — prior centred near a half, a likelihood that saw three successes in ten trials and votes for 0.30, and a posterior at a data-weighted compromise.

Figure
Bayesian updating of a rate θ: an informative prior Beta(3,3) centred at 0.50 is multiplied by the likelihood of 3 successes in 10 trials (peaked at 0.30), giving the posterior Beta(6,10) with mode 0.357 and mean 0.375 — pulled toward the data because 10 observations outweigh only 4 prior pseudo-counts
posterior ∝ likelihood × prior — the update as a pointwise reweighting posterior Beta(6,10) mode 0.357 · mean 0.375 likelihood 3 of 10 → 0.30 prior Beta(3,3) centred at 0.50 0 0.30 0.50 1 rate θ probability density

The evidence Z is the quieter hero of the same picture: the area under the product curve before renormalising — the likelihood averaged over the prior. A model that smears its prior over values it never needs pays for it, since most of that mass lands where the likelihood is small and Z is diluted — Occam's razor made quantitative, and why the integral that normalises the posterior also arbitrates between theories. The second figure shows the update in its most classical guise: three Gaussians in which combining a prior measurement with a fresh one adds their precisions and sharpens the estimate.

Figure
Normal–Normal combination of two Gaussian measurements of one quantity: a prior 100.0 ± 2.0 and a measurement 103.0 ± 1.0 combine by adding precisions (1/2² + 1/1² = 1.25) to give the posterior 102.4 ± 0.894, narrower than either input and pulled toward the sharper measurement
Precisions add: two measurements sharpen into one posterior posterior 102.4 ± 0.894 prior 100.0 ± 2.0 measurement 103.0 ± 1.0 96 100 104 108 quantity θ probability density

Two threads run through both pictures. The temporal one: yesterday's posterior is today's prior, so the answer is blind to whether data came in one batch or one point at a time. The asymptotic one: pile up enough data and the likelihood overwhelms any sensible prior, the posterior collapses to a sharp Gaussian at the MLE, and the Bayesian answer merges with the frequentist one.

The core equations

The whole framework is one theorem, read as a statement about parameters. The forward model gives the probability of data given θ; Bayes' theorem turns it into the probability of θ given the data, with Z the normalising integral.

p(θ∣x)=L(x∣θ) π(θ)Z,Z=∫L(x∣θ) π(θ) dθ(1)p(\theta\mid x)=\frac{\mathcal{L}(x\mid\theta)\,\pi(\theta)}{Z},\qquad Z=\int \mathcal{L}(x\mid\theta)\,\pi(\theta)\,\mathrm{d}\theta \qquad (1)

where: θ — unknown parameter(s), a random variable; x — observed data, fixed once collected; π(θ) — prior, belief before the data; 𝓛(x∣θ) — likelihood, the probability (density) of the data given θ; p(θ∣x) — posterior, belief after the data; Z — marginal likelihood, or evidence. Reading: multiply prior by likelihood, divide by the constant Z, and the result is everything you now know about θ. For independent, identically distributed observations the likelihood factorises, making data additive on a log scale.

L(x1:n∣θ)=∏i=1np(xi∣θ)(2)\mathcal{L}(x_{1:n}\mid\theta)=\prod_{i=1}^{n} p(x_i\mid\theta) \qquad (2)

where: x₁:ₙ — a sample of n independent observations; p(xᵢ∣θ) — the sampling density of a single datum. Reading: the log-likelihood is a sum, so information accumulates linearly in n. When prior and likelihood belong to matched families the algebra closes and the posterior is the prior's shape with shifted parameters. The Beta–Binomial pair is the canonical case: a Beta prior over a probability θ meets binomial count data and returns a Beta posterior.

θ∼Beta(α,β),    x∼Binomial(n,θ)  ⇒  θ∣x∼Beta(α+x,  β+n−x)(3)\theta\sim\mathrm{Beta}(\alpha,\beta),\;\; x\sim\mathrm{Binomial}(n,\theta)\;\Rightarrow\; \theta\mid x\sim\mathrm{Beta}(\alpha+x,\;\beta+n-x) \qquad (3)

where: α,β — prior pseudo-counts of successes and failures; x — observed successes in n trials. Reading: updating is arithmetic — add successes to α and failures to β, exposing the prior as (α+β) trials of imagined data. The Gamma–Poisson pair does the same for a counting rate — a detector logging events over a fixed exposure.

θ∼Gamma(a,b),    x∼Poisson(θ)  ⇒  θ∣x∼Gamma(a+x,  b+1)(4)\theta\sim\mathrm{Gamma}(a,b),\;\; x\sim\mathrm{Poisson}(\theta)\;\Rightarrow\;\theta\mid x\sim\mathrm{Gamma}(a+x,\;b+1) \qquad (4)

where: a,b — prior shape and rate of the Gamma; x — observed Poisson count with mean θ. Reading: the count adds to the shape, one unit of exposure to the rate — again a pure shift. The Normal–Normal pair is the workhorse of metrology: a Gaussian prior meets Gaussian measurements of known noise, and the posterior mean is a precision-weighted average while the posterior precision is the plain sum of the inputs.

θpost=σ−2 xˉ+σ0−2 μ0n σ−2+σ0−2,σpost−2=σ0−2+n σ−2(5)\theta_{\mathrm{post}}=\frac{\sigma^{-2}\,\bar x+\sigma_0^{-2}\,\mu_0}{n\,\sigma^{-2}+\sigma_0^{-2}},\qquad \sigma_{\mathrm{post}}^{-2}=\sigma_0^{-2}+n\,\sigma^{-2} \qquad (5)

where: μ₀,σ₀² — prior mean and variance of θ; x̄ — sample mean of n observations of known noise σ²; σ²post — posterior variance. Reading: precision (inverse variance) is the currency of information and it simply adds, so every measurement sharpens the estimate. Turning a posterior into a single number is a decision, and the best number depends on the loss you pay for being wrong — the Bayes estimator minimises the posterior expected loss.

θ^Bayes=arg⁡min⁡θ^ E[L(θ,θ^)∣x](6)\hat\theta_{\mathrm{Bayes}}=\arg\min_{\hat\theta}\,\mathbb{E}\big[L(\theta,\hat\theta)\mid x\big] \qquad (6)

where: L(θ,θ̂) — loss incurred by estimate θ̂ when the truth is θ; 𝔼[·∣x] — expectation against the posterior. Reading: squared-error loss picks the posterior mean, absolute-error loss the median, and 0–1 loss the mode (the MAP estimate) — one posterior, three summaries. To predict a future datum rather than estimate a parameter, integrate the parameter out against the posterior, folding parameter uncertainty into the forecast rather than betting on one best-fit value.

p(x~∣x)=∫p(x~∣θ) p(θ∣x) dθ(7)p(\tilde x\mid x)=\int p(\tilde x\mid\theta)\,p(\theta\mid x)\,\mathrm{d}\theta \qquad (7)

where: x̃ — a future, unobserved datum; p(x̃∣θ) — its sampling density. Reading: averaging the forward model over the posterior makes the prediction broader than any point estimate would. Comparing two whole models uses the ratio of their evidences — the Bayes factor — weighing how well each predicted the data after averaging over its own parameters.

B12=Z1Z2=∫L(x∣θ1,M1) π(θ1∣M1) dθ1∫L(x∣θ2,M2) π(θ2∣M2) dθ2(8)B_{12}=\frac{Z_1}{Z_2}=\frac{\int \mathcal{L}(x\mid\theta_1,M_1)\,\pi(\theta_1\mid M_1)\,\mathrm{d}\theta_1}{\int \mathcal{L}(x\mid\theta_2,M_2)\,\pi(\theta_2\mid M_2)\,\mathrm{d}\theta_2} \qquad (8)

where: M₁,M₂ — competing models, each with its own parameters and prior; Zᵢ — the evidence of model i; B₁₂ — the factor by which the data shift the odds toward M₁. Reading: because each Zᵢ averages over parameter space, a model wasting prior mass on values the data reject is penalised automatically. Finally, the curvature of the log-likelihood at its peak sets how tightly the data pin the parameter. The Fisher information is that expected curvature; its inverse is the smallest achievable variance.

Ijk(θ)=− E ⁣[∂2ln⁡p(x∣θ)∂θj ∂θk](9)\mathcal{I}_{jk}(\theta)=-\,\mathbb{E}\!\left[\frac{\partial^2 \ln p(x\mid\theta)}{\partial\theta_j\,\partial\theta_k}\right] \qquad (9)

where: 𝓘ⱼₖ — the (j,k) element of the Fisher information matrix; θⱼ — the j-th parameter component. Reading: a sharply peaked log-likelihood is a large 𝓘 and a small error bar — why (9) sets the asymptotic posterior width below.

Derivation APosterior as likelihood times prior and the role of the evidence

The goal here is to show that Bayes' theorem for parameters is not a new axiom but a two-line consequence of the definition of conditional probability, and to pin down exactly what the evidence Z is and why it must be an integral over the whole parameter space. We factorise the joint density both ways, equate them, and normalise to force Z into existence.

Step 1. Write the joint density and factorise it both ways: p(θ,x) = p(θ∣x)p(x) = p(x∣θ)π(θ), the first conditioning on the data, the second on the parameter. Why: the entire theorem is the statement that these two decompositions of one joint density must agree — no extra physics, only bookkeeping of conditional probability.

Step 2. Equate the two factorisations and solve for the posterior:

p(θ∣x)=p(x∣θ) π(θ)p(x)=L(x∣θ) π(θ)p(x)p(\theta\mid x)=\frac{p(x\mid\theta)\,\pi(\theta)}{p(x)}=\frac{\mathcal{L}(x\mid\theta)\,\pi(\theta)}{p(x)}

Here p(x∣θ), read as a function of θ with x fixed, is precisely the likelihood 𝓛(x∣θ). Why: the same conditional density is a probability model for the data when θ is fixed, but a likelihood when the data are fixed — this change of viewpoint is the whole trick of inference.

Step 3. The denominator p(x) does not depend on θ, so demand that the posterior integrate to one:

∫p(θ∣x) dθ=1  ⇒  p(x)=∫L(x∣θ) π(θ) dθ≡Z\int p(\theta\mid x)\,\mathrm{d}\theta=1 \;\Rightarrow\; p(x)=\int \mathcal{L}(x\mid\theta)\,\pi(\theta)\,\mathrm{d}\theta \equiv Z

Why: normalisation forces Z to be the prior-weighted integral of the likelihood — the average probability of the data under the model, what a Bayes factor compares.

Step 4. Read off the two faces of Z. For a single model it is a constant in θ that cancels in every posterior expectation, so it can be ignored while estimating parameters; across models it differs, and B₁₂ = Z₁/Z₂ carries the verdict. Why: parameter estimation can use an unnormalised posterior, yet model comparison cannot — the very integral you were allowed to drop is the model's score. *Verified: for the Beta–Binomial pair Z = C(n,x) B(α+x, β+n−x)/B(α,β), a finite closed form, so Z is a genuine number and the posterior a proper Beta.*

Go deeper · Derivation B · 1 minNormal–Normal update and the addition of precisions
Derivation BNormal–Normal update and the addition of precisions

The aim is to derive equation (5) from scratch by completing the square in the log-posterior, and in doing so to expose the single most useful fact in all of measurement combination: that precisions, not variances, are what add. Because a quadratic log-density is exactly a Gaussian, no approximation is involved for this conjugate pair.

Step 1. Add the log-likelihood of n observations to the log-prior, discarding θ-independent constants:

ln⁡p(θ∣x)=−12σ2∑i=1n(xi−θ)2−12σ02(θ−μ0)2+const\ln p(\theta\mid x)=-\frac{1}{2\sigma^2}\sum_{i=1}^{n}(x_i-\theta)^2-\frac{1}{2\sigma_0^2}(\theta-\mu_0)^2+\text{const}

Why: on the log scale the update is addition, so the posterior is the sum of two parabolas — itself a parabola, which guarantees a Gaussian posterior.

Step 2. Expand and collect powers of θ:

ln⁡p(θ∣x)=−12(nσ−2+σ0−2)θ2+(σ−2∑ixi+σ0−2μ0)θ+const\ln p(\theta\mid x)=-\tfrac12\big(n\sigma^{-2}+\sigma_0^{-2}\big)\theta^2+\big(\sigma^{-2}\textstyle\sum_i x_i+\sigma_0^{-2}\mu_0\big)\theta+\text{const}

Why: only the quadratic and linear coefficients survive, and for a Gaussian they are the whole story — the quadratic coefficient is the precision, their ratio is the mean.

Step 3. Read the precision off as twice the quadratic coefficient and the mean as the linear-to-quadratic ratio, using Σᵢxᵢ = nx̄:

σpost−2=σ0−2+nσ−2,θpost=σ−2 nxˉ+σ0−2μ0nσ−2+σ0−2\sigma_{\mathrm{post}}^{-2}=\sigma_0^{-2}+n\sigma^{-2},\qquad \theta_{\mathrm{post}}=\frac{\sigma^{-2}\,n\bar x+\sigma_0^{-2}\mu_0}{n\sigma^{-2}+\sigma_0^{-2}}

Why: the posterior precision is the plain sum of prior and n data precisions — information adds — and the mean weights each source by its sharpness, so a tight measurement dominates a loose prior.

Step 4. For two independent measurements x₁ ± σ₁ and x₂ ± σ₂ under a flat prior the estimate is (σ₁⁻²x₁ + σ₂⁻²x₂)/(σ₁⁻² + σ₂⁻²) with variance (σ₁⁻² + σ₂⁻²)⁻¹. Why: this is the inverse-variance weighted average of metrology, and the combined error bar is always smaller than either input. *Verified: with μ₀ = 100.0, σ₀ = 2.0 and one datum 103.0 ± 1.0, precision 1/4 + 1 = 1.25 gives σ_post = 0.894 and θ_post = 102.4 — matching the second figure.*

Go deeper · Derivation C · 1 minAsymptotic Gaussian posterior and the Laplace evidence
Derivation CAsymptotic Gaussian posterior and the Laplace evidence

The purpose here is to explain why, with enough data, the posterior forgets the prior and becomes a Gaussian centred on the maximum-likelihood estimate with covariance the inverse Fisher information — the Bernstein–von Mises theorem — and to extract as a bonus the Laplace approximation to the evidence, whose determinant factor is the Occam penalty.

Step 1. Let θ★ maximise g(θ) = ln[𝓛(x∣θ)π(θ)] and Taylor-expand to second order:

g(θ)≈g(θ⋆)+12(θ−θ⋆) ⁣⊤∇2g(θ⋆)(θ−θ⋆)g(\theta)\approx g(\theta^\star)+\tfrac12(\theta-\theta^\star)^{\!\top}\nabla^2 g(\theta^\star)(\theta-\theta^\star)

The linear term vanishes at the stationary point θ★. Why: near any smooth peak a function looks parabolic, and a parabolic log-density is exactly a Gaussian, so the posterior is locally Gaussian whatever the model.

Step 2. Exponentiate to identify the local Gaussian with precision matrix H = −∇²g(θ★):

p(θ∣x)≈N ⁣(θ; θ⋆, H−1),Hjk=−∂2g∂θj∂θk∣θ⋆p(\theta\mid x)\approx\mathcal{N}\!\big(\theta;\,\theta^\star,\,H^{-1}\big),\qquad H_{jk}=-\frac{\partial^2 g}{\partial\theta_j\partial\theta_k}\bigg|_{\theta^\star}

Why: the curvature of the log-posterior is its precision — a sharp peak (large H) is a narrow posterior (small H⁻¹): more curvature, smaller error bar.

Step 3. Integrate the Gaussian to approximate the evidence:

Z≈L(x∣θ⋆) π(θ⋆) (2π)d/2 ∣H∣−1/2Z\approx \mathcal{L}(x\mid\theta^\star)\,\pi(\theta^\star)\,(2\pi)^{d/2}\,|H|^{-1/2}

with d the number of parameters. Why: the best-fit likelihood is multiplied by a volume factor |H|⁻¹/², the width of the posterior peak — this shrinking volume is the Occam factor penalising extra parameters or wider priors.

Step 4. Take n independent observations so that H = −∇²(ln𝓛) − ∇²(lnπ) → n𝓘(θ★) as n → ∞, since the log-likelihood curvature grows like n while the fixed prior curvature does not:

p(θ∣x)  → n→∞   N ⁣(θ; θ^MLE, [ n I(θ^) ]−1)p(\theta\mid x)\;\xrightarrow{\ n\to\infty\ }\;\mathcal{N}\!\big(\theta;\,\hat\theta_{\mathrm{MLE}},\,[\,n\,\mathcal{I}(\hat\theta)\,]^{-1}\big)

Why: the prior becomes negligible against n data points, θ★ slides to the MLE, and the covariance is the inverse Fisher information — Bayesian and frequentist intervals converge, the inferential echo of the central limit theorem. *Verified: for a Gaussian likelihood with a Gaussian prior the quadratic expansion is exact, and Step 4 reproduces the precision-addition result σ_post⁻² = σ₀⁻² + nσ⁻² of the previous derivation, confirming the approximation becomes an identity in the conjugate case.*

Limiting cases

The general rule p(θ∣x) ∝ 𝓛(x∣θ)π(θ) collapses to familiar textbook answers in the regimes below.

CaseResultMeaning
Large data / likelihood dominancep(θ∣x) → 𝒩(θ; θ̂_MLE, [n𝓘(θ̂)]⁻¹) as n → ∞Bernstein–von Mises: the prior washes out and the posterior is asymptotically Gaussian about the MLE with covariance the inverse of n times the Fisher information.
Flat (uniform) priorp(θ∣x) ∝ 𝓛(x∣θ), so θ̂_MAP = θ̂_MLEWith π constant over the support the posterior is the normalised likelihood, and the MAP estimate coincides with maximum likelihood.
Conjugate priorπ(θ) ∈ 𝓕 ⇒ p(θ∣x) ∈ 𝓕For exponential-family likelihoods a conjugate prior keeps the posterior in the same family; the update is an arithmetic shift of parameters.
Single Gaussian datum, Gaussian priorσ_post⁻² = σ₀⁻² + σ⁻²Precisions add: combining a prior of precision σ₀⁻² with a measurement of precision σ⁻² always sharpens the estimate.
No data (n = 0)p(θ∣x) = π(θ)With nothing observed the posterior equals the prior — the sanity check every update must pass.

Worked examples

1 · Detector efficiency from a Beta–Binomial update (fully worked)

A detector fires on 3 of 10 injected test pulses. Starting from a uniform prior Beta(1,1) on the efficiency θ, find the posterior, its point estimates, and a 95% credible interval, and compare with the naïve estimate 3/10.

Step 1. The uniform prior is Beta(1,1), so α = β = 1. With x = 3 successes in n = 10, conjugacy (equation 3) gives the posterior Beta(1+3, 1+7) = Beta(4, 8) directly — the uniform prior contributed one pseudo-success and one pseudo-failure.

Step 2. From the standard Beta formulas the mean is α'/(α'+β') = 4/12 = 0.333, the MAP mode (α'−1)/(α'+β'−2) = 3/10 = 0.300, and the standard deviation √(32/1872) = 0.131. The mean sits above the naïve 3/10 because the uniform prior nudges toward 1/2 — Laplace's rule of succession, mean = (x+1)/(n+2) = 4/12.

Step 3. The 95% central credible interval is the 2.5% and 97.5% quantiles of Beta(4,8): [0.109, 0.610]. It is wide and asymmetric, honestly reflecting thin data, and it never touches 0 or 1 — so a run of all-successes could never fool it into false certainty, unlike the raw ratio.

Step 4. The mean (0.333), median (≈ 0.32), and mode (0.300) differ because Beta(4,8) is right-skewed: quote the mean under squared-error loss, the mode for the most probable efficiency. Posterior Beta(4,8): mean 0.333, mode 0.300, sd 0.131, 95% CI [0.109, 0.610] — the answer is the whole distribution; these are just summaries.

2 · Combining two measurements with the Normal–Normal model (fully worked)

Two independent experiments report the same mass as 100.0 ± 2.0 and 103.0 ± 1.0 (arbitrary units, Gaussian errors). Under a flat prior, what is the combined estimate and its uncertainty?

Step 1. Treat the first measurement as the prior, μ₀ = 100.0, σ₀ = 2.0, and the second as the datum, 103.0 ± 1.0 (n = 1). Work in precisions: prior 1/σ₀² = 0.25, data 1/σ² = 1.00.

Step 2. Add the precisions (equation 5): the posterior precision is 0.25 + 1.00 = 1.25, so the variance is 1/1.25 = 0.80 and the combined uncertainty √0.80 = 0.894 — smaller than either 2.0 or 1.0, as adding information must make it.

Step 3. The mean is the precision-weighted average: (0.25·100.0 + 1.00·103.0)/1.25 = 128.0/1.25 = 102.4, landing nearer the sharper measurement, which carries four times the weight. Combined result: 102.4 ± 0.894 — the inverse-variance average of metrology, identical to the second figure.

3 · An Occam penalty in a Bayes factor

Model 1 has one parameter; rival Model 2 fits equally well but carries a second parameter the data do not need. Each parameter is constrained to a posterior width Δθ_post = 0.2 inside a prior of width Δθ_prior = 10. Which model does the evidence favour, and by how much?

Step 1. For a flat prior the evidence factorises as Z ≈ 𝓛(θ★) × ∏ⱼ (Δθ_post/Δθ_prior). Each constrained parameter contributes an Occam factor 0.2/10 = 0.02, the fraction of prior space the data leave standing.

Step 2. Both models share the same best-fit likelihood 𝓛(θ★) by assumption, and both spend one Occam factor of 0.02 on the parameter they genuinely need. Model 2 spends a second factor of 0.02 on its superfluous parameter, so Z₂ = Z₁ × 0.02.

Step 3. The Bayes factor is B₁₂ = Z₁/Z₂ = 1/0.02 = 50: the data favour the simpler model by 50 to 1 even though the complex model fits equally well point-for-point, because the extra parameter's unused prior volume is charged against it. B₁₂ = 50, a quantitative Occam's razor with no penalty term inserted by hand.

Physics of the portrait

Bayesian inference carries two names — Thomas Bayes, who wrote down the special case, and Pierre-Simon Laplace, who built the general theory and made it the working language of physical science, first using it to weigh planets and to ask how sure we can be that the Sun will rise.

Origins · 1 minPierre-Simon Laplace and how this was discovered
PL
Pierre-Simon Laplace
1749–1827 · the founder of Bayesian probability as the calculus of scientific inference

Pierre-Simon Laplace did far more with inverse probability than the Reverend Thomas Bayes, whose 1763 essay (published posthumously by Richard Price) treated only a uniform prior on a binomial rate. Working independently and at first in ignorance of Bayes, Laplace stated the general principle in 1774: the probability of a cause given an effect is proportional to the probability of the effect given that cause. Over forty years he made it an engine of physical science — estimating the masses of Saturn and Jupiter from the planets' mutual perturbations and quoting what were in effect posterior error bars (his Saturn value stands to within one percent), deriving the rule of succession (x+1)/(n+2), and computing the celebrated odds that the Sun would rise tomorrow. His 1812 Théorie analytique des probabilités laid out the whole apparatus — priors, likelihoods, the asymptotic Gaussian posterior now named for Bernstein and von Mises, and least squares read as a Bayesian estimate under a Gaussian error law. The approach then fell from favour for a century under the frequentist critique of "subjective" priors, until Harold Jeffreys in the 1930s and later E. T. Jaynes restored it to physics. Every posterior error bar a modern cosmologist quotes is, in the end, Laplace's inverse probability on a bigger computer.

History
1763
— Bayes' essay on inverse probability is read to the Royal Society by Richard Price
1774
— Laplace states the general principle of probability of causes from effects
1810
— Laplace applies inverse probability with a Gaussian error law to estimate planetary masses
1812
— Laplace publishes Théorie analytique des probabilités, the first systematic treatise
1922
— Fisher formalises the likelihood and maximum-likelihood estimation
1939
— Jeffreys' Theory of Probability revives Bayesian methods for the physical sciences
2015
— LIGO applies full Bayesian parameter estimation to the first gravitational-wave detection

Going further

Texts. D. S. Sivia & J. Skilling, Data Analysis: A Bayesian Tutorial, 2nd ed. (Oxford, 2006) — the physicist's entry point · E. T. Jaynes, Probability Theory: The Logic of Science (Cambridge, 2003) · A. Gelman et al., Bayesian Data Analysis, 3rd ed. (CRC, 2013) · Stone & Goldbart, Mathematics for Physics (Cambridge, 2009) · Arfken, Weber & Harris, Mathematical Methods for Physicists, 7th ed., Ch. 23 · Riley, Hobson & Bence, Mathematical Methods for Physics and Engineering, 3rd ed., Ch. 30–31.

Provenance. The conjugate updates follow Sivia & Skilling and Gelman et al.; the decision-theoretic estimators, Laplace/Bernstein–von Mises asymptotics, and Occam-factor reading follow Jaynes and MacKay. A verifier should: (1) re-derive the Beta–Binomial, Gamma–Poisson, and Normal–Normal posteriors by multiplying densities; (2) confirm the Beta(4,8) summaries mean 1/3, mode 0.300, sd 0.131, 95% interval [0.109, 0.610]; (3) check the Normal–Normal combination 102.4 ± 0.894 from precisions 0.25 and 1.00; (4) verify the Laplace evidence and the Occam factor 0.02 giving B₁₂ = 50; (5) confirm the Bernstein–von Mises covariance [n𝓘]⁻¹ against definition (9). All quoted values are standard textbook results.

Prerequisites. Conditional Probability and Bayes' Theorem · Random Variables and Probability Distributions · Joint, Marginal and Conditional Distributions · Expectation, Variance and Moments · Gaussian (Normal) Distribution. Leads to. Central Limit Theorem · Multivariate Gaussian Distribution · Moment-Generating and Characteristic Functions.

Ask the tutor: Ask me to run a full Beta–Binomial or Normal–Normal update on your own numbers, or to walk through why a Bayes factor rewards the simpler model.

What is the difference between a likelihood and a probability?

They are the same function read from opposite ends. The forward model p(x∣θ) is a probability distribution over data x when θ is fixed — its integral over x is one. Read as a function of θ with the data fixed, the identical expression is the likelihood 𝓛(x∣θ), which is not a distribution in θ: it need not integrate to one and has no units of probability. Bayesian inference turns that likelihood back into a genuine distribution over θ by multiplying by a prior and dividing by the evidence. Frequentist maximum likelihood stops at the peak of 𝓛; the Bayesian keeps the whole shape.

How do I choose a prior, and doesn't it make the answer subjective?

The prior encodes what you knew before the experiment; choosing it is a modelling decision — you state it, and readers can check how much it matters. Three routes dominate. If earlier data exist, use their posterior as your prior (just sequential updating). If you want the data to speak, use a weakly informative or "objective" prior — a Jeffreys prior or a broad Gaussian — and the large-data limit guarantees the answer stops depending on it. If physical knowledge constrains the parameter (a mass is positive), build that in. The safeguard is a sensitivity analysis: rerun with a different reasonable prior. If the conclusion moves, the data were too weak to settle the question — itself worth knowing.

What is a conjugate prior and why is it so convenient?

A prior is conjugate to a likelihood when the posterior stays in the same functional family, so the update is a change of parameters rather than an integral. Beta is conjugate to the Binomial, Gamma to the Poisson, Normal to a Normal likelihood of known variance. The convenience is twofold: a closed-form posterior, and a transparent reading — the prior acts like imagined data, so Beta(α,β) is "α+β prior trials with α successes" and real data just adds real counts. The limitation is rigidity; for realistic multi-parameter models you abandon conjugacy and sample the posterior with MCMC or nested sampling.

When do the posterior mean, median, and mode disagree?

They agree whenever the posterior is symmetric and unimodal — a Gaussian being the standard case. They part company when it is skewed, as Beta(4,8) is: mean 0.333, median ≈ 0.32, mode 0.300. Which to quote is set by the loss you face: squared-error rewards the mean, absolute-error the median, 0–1 loss the mode (MAP). Any single number throws information away; when the posterior is markedly skewed or multimodal, show the distribution or a credible interval, not a lone point.

Is maximum likelihood just a special case of Bayesian estimation?

Almost — it is the MAP estimate under a flat prior. If π(θ) is constant, the posterior is proportional to the likelihood, so its mode is exactly the maximum-likelihood estimate; and by Bernstein–von Mises the two agree asymptotically on the error bar too, since the posterior covariance tends to the inverse Fisher information that sets the frequentist variance. The differences are conceptual and finite-sample: Bayes reports a full distribution, can fold in prior information, and handles nuisance parameters by integration (marginalising) rather than maximisation (profiling), which matters when parameters are correlated.

What is the evidence Z actually good for if it cancels in parameter estimation?

Within a single model Z is a constant in θ, drops out of every posterior expectation, and is never needed to estimate parameters — which is why samplers like Metropolis–Hastings work with the unnormalised posterior. It matters when you compare models: Z is the prior-averaged probability a model assigned to the data, so Z₁/Z₂ is the Bayes factor, the amount the data shift the odds between two theories. Because Z integrates over each model's whole parameter space, it automatically penalises wasted prior mass on values the data reject — the built-in Occam's razor. Computing it is the hard part, and nested sampling was invented to estimate this one integral efficiently.

Before this

Where this leads