45. Prior and posterior distributions
Everything in §4 and §5 treated \theta as a fixed unknown constant. That's why a confidence interval couldn't say "there's a 95% chance \mu is in here" — \mu isn't random, so it has no probability distribution.
Bayesian statistics makes one change, and everything follows from it:
Treat the parameter itself as a random variable, with a distribution representing your uncertainty about it.
Bayes' theorem, for parameters
The same theorem from §0, now with \theta in place of an event:
p(\theta \mid \text{data}) = \frac{p(\text{data} \mid \theta)\,p(\theta)}{p(\text{data})}
Each piece has a name:
| Term | Name | Meaning |
|---|---|---|
| p(\theta) | prior | belief about \theta before seeing data |
| p(\text{data} \mid \theta) | likelihood | how well each \theta explains the data |
| p(\theta \mid \text{data}) | posterior | belief after seeing data |
| p(\text{data}) | evidence | normalising constant |
The denominator doesn't depend on \theta — it's just whatever makes the posterior integrate to 1. So in practice you work with
\boxed{\;p(\theta \mid \text{data}) \;\propto\; p(\text{data} \mid \theta)\,p(\theta)\;}
Posterior ∝ likelihood × prior. That is the entire machinery of Bayesian inference. Everything else is computation.
What changes
The output is a whole distribution, not a point estimate plus a standard error. From it you can read off anything you want:
- Posterior mean E[\theta \mid \text{data}] — the usual point estimate.
- Posterior mode (MAP) — the most probable value.
- Credible interval — a range containing 95% of the posterior probability.
And crucially, you can now say:
P(0.3 < \theta < 0.5 \mid \text{data}) = 0.95
This is a genuine probability statement about \theta, which a frequentist confidence interval never was (§4). The price is that it depends on your prior.
Choosing a prior
Informative priors encode real knowledge — previous studies, physical constraints, expert judgement. If a coin came from a bank, a prior concentrated near p = 0.5 is honest, not a cheat.
Weakly informative priors rule out the absurd without pinning things down. Usually the sensible default.
Uninformative / flat priors try to "let the data speak". Uniform on [0,1] for a probability, for instance. But be careful: flat is not reparameterisation-invariant. A uniform prior on \theta is not uniform on \log\theta or on \theta^2, so "no information" isn't a well-defined notion. (Jeffreys priors, p(\theta) \propto \sqrt{I(\theta)} using the Fisher information from §4, are constructed precisely to be invariant.)
The honest position: the prior is an assumption, like the choice of likelihood, the choice of model, or the decision to trim outliers. It should be stated and its influence checked — that's what sensitivity analysis means. Frequentist analyses make assumptions too; they're just less explicit about them.
The data swamps the prior
This is the reassuring part. As n grows, the likelihood becomes sharply peaked while the prior stays fixed, so the posterior is increasingly dominated by the data.
Two analysts with quite different priors converge on nearly the same posterior given enough data. The prior matters most exactly when data is scarce — which is also when you most want to use whatever outside knowledge you have.
The one exception: a prior that assigns zero probability to some region keeps the posterior at zero there forever, no matter what the data says. Multiplying by zero is irreversible. Never rule anything out absolutely unless it's genuinely impossible.
Worked example: a coin
Prior: p \sim \text{Beta}(2, 2) — a gentle hump at 0.5, mildly favouring a fair coin without insisting on it.
Data: 7 heads in 10 flips.
Likelihood: p(\text{data} \mid p) \propto p^7(1-p)^3.
Posterior:
p(p \mid \text{data}) \propto \underbrace{p^7(1-p)^3}_{\text{likelihood}} \times \underbrace{p^{2-1}(1-p)^{2-1}}_{\text{prior}} = p^{8}(1-p)^{4}
which is the kernel of a \text{Beta}(9, 5).
The prior was Beta, the posterior is Beta — that's conjugacy, the subject of the next lesson, and it's why this worked out on paper.
Reading off the summaries:
E[p \mid \text{data}] = \frac{9}{9+5} = 0.643
Compare with the MLE of 7/10 = 0.7. The posterior mean is pulled toward the prior mean of 0.5 — that's shrinkage, from §6, appearing here as a natural consequence rather than an added penalty.
A 95% credible interval is roughly (0.39, 0.86), and it means exactly what it looks like: given the prior and data, there's a 95% probability p lies in that range.
With 700 heads in 1000 flips, the posterior becomes \text{Beta}(702, 302), with mean 0.699 — essentially the MLE. The same prior that moved the estimate by 0.06 at n = 10 moves it by 0.001 at n = 1000.
Doing it in Python
The update, and how the prior fades:
from scipy.stats import beta
prior_a, prior_b = 2, 2
print(f"{'data':>16} {'posterior':>18} {'mean':>8} {'MLE':>8}")
for heads, n in [(0, 0), (7, 10), (70, 100), (700, 1000)]:
a, b = prior_a + heads, prior_b + (n - heads)
mle = heads / n if n else float("nan")
print(f"{f'{heads}/{n}':>16} {f'Beta({a}, {b})':>18} "
f"{a/(a+b):>8.4f} {mle:>8.4f}")
print("\nThe posterior mean starts at the prior mean (0.5) and converges to the MLE.")
Computing the posterior numerically on a grid — the method that works for any prior, conjugate or not:
import numpy as np
from scipy.stats import beta
from scipy.integrate import trapezoid
grid = np.linspace(0, 1, 2001)
heads, n = 7, 10
prior = beta.pdf(grid, 2, 2)
likelihood = grid**heads * (1 - grid) ** (n - heads)
post = prior * likelihood
post /= trapezoid(post, grid) # normalise numerically
exact = beta.pdf(grid, 2 + heads, 2 + n - heads)
print("grid posterior matches Beta(9, 5):", bool(np.allclose(post, exact, atol=1e-6)))
mean = trapezoid(grid * post, grid)
print(f"posterior mean : {mean:.4f}")
print(f"posterior mode : {grid[np.argmax(post)]:.4f}")
print(f"MLE : {heads/n:.4f}")
cdf = np.cumsum(post) * (grid[1] - grid[0])
lo = grid[np.searchsorted(cdf, 0.025)]
hi = grid[np.searchsorted(cdf, 0.975)]
print(f"95% credible interval: ({lo:.4f}, {hi:.4f})")
Different priors, same data — and then the same priors with more data:
import numpy as np
from scipy.stats import beta
priors = {
"flat Beta(1,1)": (1, 1),
"gentle Beta(2,2)": (2, 2),
"strong fair Beta(50,50)": (50, 50),
"sceptical Beta(1,9)": (1, 9),
}
for label, n in [("n = 10 (7 heads)", (7, 10)), ("n = 1000 (700 heads)", (700, 1000))]:
heads, total = n
print(f"\n{label}")
for name, (a0, b0) in priors.items():
a, b = a0 + heads, b0 + total - heads
lo, hi = beta.ppf([0.025, 0.975], a, b)
print(f" {name:>24}: mean {a/(a+b):.4f} 95% CI ({lo:.3f}, {hi:.3f})")
print("\nAt n=10 the priors disagree noticeably. At n=1000 they have converged.")
And the one prior that never washes out:
import numpy as np
from scipy.integrate import trapezoid
grid = np.linspace(0, 1, 1001)
# A prior that flatly rules out p > 0.6
prior = np.where(grid <= 0.6, 1.0, 0.0)
for heads, n in [(7, 10), (700, 1000), (99_000, 100_000)]:
like = np.exp(heads * np.log(grid + 1e-300) + (n - heads) * np.log(1 - grid + 1e-300))
post = prior * like
post /= trapezoid(post, grid)
mean = trapezoid(grid * post, grid)
print(f"{heads}/{n}: posterior mean {mean:.4f}, "
f"P(p > 0.6) = {trapezoid(post[grid > 0.6], grid[grid > 0.6]):.4f}")
print("\nEven with 99% heads in 100,000 flips, the posterior still says p <= 0.6.")
print("Zero prior probability is permanent -- Bayes can never undo a multiplication by 0.")
Your turn
1. Prior \text{Beta}(1,1), data 3 heads in 5 flips. Find the posterior and its mean.
2. Why does the prior matter less as n grows?
3. What's wrong with assigning a prior of zero to a region you think is merely unlikely?
Solutions
1. \text{Beta}(1,1) is the uniform distribution on [0,1], so the posterior is proportional to the likelihood alone:
p(p \mid \text{data}) \propto p^3(1-p)^2 \times 1 = p^{4-1}(1-p)^{3-1}
That's a \text{Beta}(4, 3), and
E[p \mid \text{data}] = \frac{4}{4+3} = \frac{4}{7} \approx 0.571
Interesting: the MLE is 3/5 = 0.6, and the posterior mean is lower. Even a flat prior pulls toward 0.5, because it spreads weight uniformly across a region that includes values below the MLE. The general Beta-Binomial rule is
E[p \mid \text{data}] = \frac{a + k}{a + b + n}
which for \text{Beta}(1,1) becomes \frac{k+1}{n+2} — Laplace's rule of succession, the classic fix for the zero-count problem.
2. Because the likelihood concentrates while the prior doesn't.
The log-posterior is \log(\text{likelihood}) + \log(\text{prior}). The likelihood term is a sum over n observations, so it grows with n; the prior term is a single fixed quantity. Their ratio goes to infinity.
Geometrically: the likelihood's width shrinks like 1/\sqrt n (§3), so it becomes a spike. Multiplying a fixed smooth prior by a spike gives back essentially the spike, shifted negligibly.
This is why Bayesian and frequentist answers converge at large n — the Bernstein–von Mises theorem makes it precise: the posterior becomes approximately Normal, centred at the MLE, with variance 1/I_n(\theta), which is exactly the frequentist sampling distribution from §4.
3. Because the posterior can never recover from it:
p(\theta \mid \text{data}) \propto p(\text{data} \mid \theta) \times 0 = 0
No amount of evidence multiplies back up from zero. You have declared the region impossible, not merely unlikely, and Bayes' theorem will honour that declaration forever.
This is sometimes called Cromwell's rule, after Cromwell's plea: "think it possible that you may be mistaken."
The practical fix: use a prior that's small rather than zero in regions you doubt. A prior giving a region probability 10^{-6} still lets overwhelming evidence move the posterior there; a prior of exactly 0 does not.
The legitimate exception is genuine impossibility — a probability outside [0,1], a negative variance. Ruling those out is a fact about the parameter space, not a belief.
Check yourself in code
Compute the posterior for a coin with a Beta(2,2) prior and 7 heads in 10 flips, and verify the numerical grid posterior matches the exact conjugate answer.
Print exactly this:
posterior Beta(9, 5)
posterior mean 0.6429
MLE 0.7
grid matches exact: True
Round the posterior mean to 4 decimal places and the MLE to 4. Use a grid of
2001 points on [0,1] and scipy.integrate.trapezoid to normalise; treat the match as
holding if np.allclose passes with atol=1e-6.
import numpy as np
from scipy.stats import beta
prior_a, prior_b = 2, 2
heads, n = 7, 10
a, b = prior_a + heads, prior_b + (n - heads)
print(f"posterior Beta({a}, {b})")
print("posterior mean", round(a / (a + b), 4))
print("MLE", round(heads / n, 4))
# Build the posterior on a grid from prior * likelihood, normalise it,
# and check it matches beta.pdf(grid, a, b).
import numpy as np
from scipy.stats import beta
from scipy.integrate import trapezoid
prior_a, prior_b = 2, 2
heads, n = 7, 10
a, b = prior_a + heads, prior_b + (n - heads)
print(f"posterior Beta({a}, {b})")
print("posterior mean", round(a / (a + b), 4))
print("MLE", round(heads / n, 4))
grid = np.linspace(0, 1, 2001)
post = beta.pdf(grid, prior_a, prior_b) * grid**heads * (1 - grid) ** (n - heads)
post /= trapezoid(post, grid)
print("grid matches exact:", bool(np.allclose(post, beta.pdf(grid, a, b), atol=1e-6)))
Bayesian inference treats the parameter as random and updates a prior into a posterior: posterior ∝ likelihood × prior. The output is a full distribution, so probability statements about \theta become legitimate. The prior is an explicit assumption whose influence fades as data accumulates — except where it assigns exactly zero, which is permanent.
Next: the special priors that make this update work out in closed form.