46. Conjugate priors

The last lesson's Beta prior produced a Beta posterior. That wasn't luck — it's conjugacy, and it's what makes Bayesian updating possible with algebra instead of computation.

The definition

A prior is conjugate to a likelihood if the posterior belongs to the same family as the prior.

\text{Beta prior} \times \text{Binomial likelihood} \;\longrightarrow\; \text{Beta posterior}

The practical payoff: updating becomes arithmetic on the parameters. No integration, no normalising constant to compute, no numerical method. And because the family is closed, you can update again and again as data arrives, each time just adjusting numbers.

Beta–Binomial

Prior: p \sim \text{Beta}(\alpha, \beta) Data: k successes in n trials Posterior: p \mid \text{data} \sim \text{Beta}(\alpha + k,\; \beta + n - k)

The derivation is one line of matching exponents:

\underbrace{p^{\alpha-1}(1-p)^{\beta-1}}_{\text{prior}} \times \underbrace{p^k(1-p)^{n-k}}_{\text{likelihood}} = p^{\alpha+k-1}(1-p)^{\beta+n-k-1}

Interpreting the parameters: \alpha acts like "prior successes" and \beta like "prior failures". A \text{Beta}(3,7) prior carries the same weight as having already seen 3 successes and 7 failures — 10 pseudo-observations. This makes priors easy to elicit and easy to explain.

E[p \mid \text{data}] = \frac{\alpha + k}{\alpha + \beta + n}

Rewrite that as a weighted average and the shrinkage becomes explicit:

E[p \mid \text{data}] = \underbrace{\frac{n}{\alpha+\beta+n}}_{w}\cdot\underbrace{\frac{k}{n}}_{\text{MLE}} + \underbrace{\frac{\alpha+\beta}{\alpha+\beta+n}}_{1-w}\cdot\underbrace{\frac{\alpha}{\alpha+\beta}}_{\text{prior mean}}

The posterior mean is a weighted average of the data and the prior, with weights given by their respective sample sizes. As n \to \infty, w \to 1 and the prior vanishes — the mechanism from the last lesson, now visible in a formula.

Normal–Normal (known variance)

Prior: \mu \sim N(\mu_0, \tau_0^2) Data: n observations with known \sigma^2, sample mean \bar x Posterior: \mu \mid \text{data} \sim N(\mu_n, \tau_n^2) where

\frac{1}{\tau_n^2} = \frac{1}{\tau_0^2} + \frac{n}{\sigma^2}, \qquad \mu_n = \tau_n^2\left(\frac{\mu_0}{\tau_0^2} + \frac{n\bar x}{\sigma^2}\right)

This is cleanest in terms of precision — the reciprocal of variance:

\text{posterior precision} = \text{prior precision} + \text{data precision}

Precisions add. Information accumulates additively, and the posterior mean is a precision-weighted average of the prior mean and the sample mean. Whichever source is more precise gets more say.

Note the posterior variance is always smaller than either input — combining two sources of information can only reduce uncertainty.

The standard table

Likelihood	Parameter	Conjugate prior	Posterior
Binomial	p	Beta(\alpha, \beta)	Beta(\alpha + k,\ \beta + n - k)
Poisson	\lambda	Gamma(\alpha, \beta)	Gamma(\alpha + \sum x_i,\ \beta + n)
Normal (known \sigma^2)	\mu	Normal	Normal
Normal (known \mu)	\sigma^2	Inverse-Gamma	Inverse-Gamma
Exponential	\lambda	Gamma(\alpha, \beta)	Gamma(\alpha + n,\ \beta + \sum x_i)
Multinomial	\mathbf p	Dirichlet	Dirichlet

The pattern behind the table: every exponential-family likelihood has a conjugate prior, obtained by treating the likelihood's own functional form as a density in the parameter. That's why the same distributions from §1 keep reappearing here.

The catch

Conjugacy is a mathematical convenience, not a statement about your beliefs.

Your actual prior may not be conjugate. If you genuinely think p is either near 0.1 or near 0.9 but not in between, no Beta distribution says that — Betas are unimodal (or U-shaped, but symmetric about that). Forcing a conjugate prior because the algebra is nicer is letting convenience distort the analysis.

Most real models have no conjugate prior at all. Logistic regression, hierarchical models, anything with several interacting parameters — none of them update in closed form.

The modern answer is computation. MCMC (§11 gives you the sampling foundations) and variational methods handle arbitrary priors and likelihoods. Conjugacy is still valuable — for intuition, for fast inner loops, for components of larger models — but it is no longer a constraint on what you can fit.

Worked example

A/B test. Version A: 30 conversions from 100 visitors. Version B: 40 from 100. Is B better?

Use a \text{Beta}(1,1) (uniform) prior for each rate.

p_A \mid \text{data} \sim \text{Beta}(31, 71), \qquad p_B \mid \text{data} \sim \text{Beta}(41, 61)

E[p_A] = \frac{31}{102} = 0.304, \qquad E[p_B] = \frac{41}{102} = 0.402

Now the question a frequentist test can't answer directly:

P(p_B > p_A \mid \text{data}) \approx 0.93

A 93% probability that B is genuinely better. That's a direct statement about the parameters, and it's exactly what a decision-maker wants — not "we reject H_0 at \alpha = 0.05".

Better still, you can compute the expected loss of choosing wrongly: E[\max(p_A - p_B, 0)] \approx 0.004, or 0.4 percentage points. Even in the 7% of cases where A is really better, it's better by very little. That combination — 93% confident, tiny downside — is a much more actionable summary than a p-value.

Note this is a sequentially valid analysis in a way p-values are not: you can look at the posterior after every visitor without inflating any error rate, because there's no repeated-testing correction to violate.

Doing it in Python

Sequential updating, one batch at a time:

from scipy.stats import beta

a, b = 2, 2          # prior: Beta(2,2)
print(f"{'batch':>18} {'posterior':>18} {'mean':>8} {'95% interval':>22}")
print(f"{'(prior)':>18} {f'Beta({a}, {b})':>18} {a/(a+b):>8.4f} "
      f"{str(tuple(beta.ppf([0.025, 0.975], a, b).round(3))):>22}")

for heads, n in [(7, 10), (12, 20), (35, 50), (150, 200)]:
    a, b = a + heads, b + (n - heads)      # conjugacy: just add
    lo, hi = beta.ppf([0.025, 0.975], a, b)
    print(f"{f'{heads}/{n}':>18} {f'Beta({a}, {b})':>18} {a/(a+b):>8.4f} "
          f"{f'({lo:.3f}, {hi:.3f})':>22}")

print("\nEach update is two additions. The interval tightens as evidence accumulates.")

The posterior mean as a weighted average, made explicit:

alpha, beta_ = 5, 15          # prior mean 0.25, worth 20 pseudo-observations
prior_mean = alpha / (alpha + beta_)

print(f"prior mean {prior_mean:.4f} (weight {alpha + beta_} pseudo-observations)\n")
print(f"{'n':>8} {'MLE':>8} {'weight on data':>16} {'posterior mean':>16}")
for k, n in [(8, 10), (40, 50), (400, 500), (4000, 5000)]:
    w = n / (alpha + beta_ + n)
    post = (alpha + k) / (alpha + beta_ + n)
    check = w * (k / n) + (1 - w) * prior_mean
    print(f"{n:>8} {k/n:>8.4f} {w:>16.4f} {post:>16.4f}")
    assert abs(post - check) < 1e-12
print("\nThe weighted-average identity holds exactly at every n.")

Normal–Normal, where precisions add:

import numpy as np

mu0, tau0 = 100.0, 15.0        # prior: mu ~ N(100, 15^2)
sigma = 10.0                   # known observation sd

print(f"{'n':>6} {'xbar':>8} {'post mean':>12} {'post sd':>10} "
      f"{'prior prec':>12} {'data prec':>11}")
for n, xbar in [(1, 130.0), (5, 130.0), (25, 130.0), (100, 130.0)]:
    prior_prec = 1 / tau0**2
    data_prec = n / sigma**2
    post_var = 1 / (prior_prec + data_prec)
    post_mean = post_var * (mu0 / tau0**2 + n * xbar / sigma**2)
    print(f"{n:>6} {xbar:>8.1f} {post_mean:>12.4f} {np.sqrt(post_var):>10.4f} "
          f"{prior_prec:>12.6f} {data_prec:>11.6f}")

print("\nPosterior mean slides from the prior (100) toward the data (130)")
print("as the data's precision overtakes the prior's.")

The A/B test, answered the way a decision-maker asks it:

import numpy as np
from scipy.stats import beta

rng = np.random.default_rng(0)

# Beta(1,1) prior + data -> Beta(1+conversions, 1+failures)
A = beta(1 + 30, 1 + 70)
B = beta(1 + 40, 1 + 60)

draws = 400_000
pa, pb = A.rvs(draws, random_state=rng), B.rvs(draws, random_state=rng)

print(f"E[p_A] = {A.mean():.4f}   E[p_B] = {B.mean():.4f}")
print(f"P(B > A) = {(pb > pa).mean():.4f}")
print(f"expected uplift  E[p_B - p_A] = {(pb - pa).mean():+.4f}")
print(f"expected loss if we pick B    = {np.maximum(pa - pb, 0).mean():.5f}")
print(f"95% interval for the difference: "
      f"({np.quantile(pb - pa, 0.025):+.4f}, {np.quantile(pb - pa, 0.975):+.4f})")

And what conjugacy can't express:

import numpy as np
from scipy.stats import beta
from scipy.integrate import trapezoid

grid = np.linspace(0.001, 0.999, 2000)

# A genuinely bimodal belief: the coin is either heavily biased one way or the other
bimodal = beta.pdf(grid, 8, 2) + beta.pdf(grid, 2, 8)
bimodal /= trapezoid(bimodal, grid)

# The best-fitting Beta (matching mean and variance) cannot reproduce it
m = trapezoid(grid * bimodal, grid)
v = trapezoid((grid - m) ** 2 * bimodal, grid)
common = m * (1 - m) / v - 1
fit = beta.pdf(grid, m * common, (1 - m) * common)

is_peak = (bimodal[1:-1] > bimodal[:-2]) & (bimodal[1:-1] > bimodal[2:])
print(f"bimodal prior: mean {m:.4f}, variance {v:.4f}")
print(f"peaks at p = {grid[1:-1][is_peak].round(3)}   ({is_peak.sum()} of them)")
print(f"best matching Beta({m*common:.2f}, {(1-m)*common:.2f}) has a single peak at "
      f"{grid[np.argmax(fit)]:.3f}")
print("\nNo Beta distribution is bimodal. If that's your actual belief,")
print("conjugacy cannot represent it -- and you should use numerical methods.")

Your turn

1. Prior \text{Gamma}(2, 1) for a Poisson rate; you observe counts 3, 5, 4. Find the posterior.

2. What does a \text{Beta}(1,1) prior represent, and what's the posterior after k successes in n trials?

3. Prior N(50, 10^2), and one observation of 70 with \sigma = 10. Find the posterior mean.

Solutions

1. For a Gamma(\alpha, \beta) prior with Poisson data, the posterior is Gamma(\alpha + \sum x_i,\ \beta + n).

Here \sum x_i = 3 + 5 + 4 = 12 and n = 3:

\lambda \mid \text{data} \sim \text{Gamma}(2 + 12,\; 1 + 3) = \text{Gamma}(14, 4)

E[\lambda \mid \text{data}] = \frac{14}{4} = 3.5

Compare with the MLE \bar x = 4. The posterior is pulled toward the prior mean of 2/1 = 2. Reading the parameters: \alpha counts prior events and \beta prior exposure, so Gamma(2,1) is worth one observation containing 2 events.

2. \text{Beta}(1,1) has density \propto p^0(1-p)^0 = 1 — it's the uniform distribution on [0,1], treating all values of p as equally plausible.

The posterior is \text{Beta}(1 + k,\; 1 + n - k), with mean

E[p \mid \text{data}] = \frac{k+1}{n+2}

This is Laplace's rule of succession. Its practical virtue shows up at the extremes: with 0 successes in 5 trials the MLE says p = 0 — claiming certainty that the event is impossible — while the posterior mean gives $1/7 \approx 0.143$, which is far more sensible from five observations.

3. Prior precision: 1/10^2 = 0.01. Data precision (one observation): 1/10^2 = 0.01.

\text{posterior precision} = 0.01 + 0.01 = 0.02 \implies \tau_n^2 = 50

\mu_n = 50\left(\frac{50}{100} + \frac{70}{100}\right) = 50(0.5 + 0.7) = 60

Exactly 60 — the midpoint of 50 and 70, because the prior and the single observation are equally precise. The posterior sd is \sqrt{50} \approx 7.07, smaller than either input's 10.

If the observation had been more precise (\sigma = 5, precision 0.04), the posterior mean would move to (0.5 + 2.8)/0.05 = 66 — closer to the data, in proportion to its precision.

Check yourself in code

Run the sequential Beta–Binomial update and verify the posterior mean equals the precision-style weighted average of the MLE and the prior mean.

Print exactly this:

posterior Beta(24, 20)
posterior mean 0.5455
weighted average 0.5455
identity holds: True

Start from a Beta(2,2) prior and apply two batches: 7 heads in 10, then 15 heads in 30. Round both means to 4 decimal places, and treat the identity as holding if they agree to within 10^{-12}.

prior_a, prior_b = 2, 2

a, b = prior_a, prior_b
for heads, n in [(7, 10), (15, 30)]:
    a, b = a + heads, b + (n - heads)

print(f"posterior Beta({a}, {b})")
print("posterior mean", round(a / (a + b), 4))

# Compute the same number as w * MLE + (1 - w) * prior_mean, where the total
# data is 22 heads in 40 and w = n / (prior_a + prior_b + n). Compare them.

prior_a, prior_b = 2, 2

a, b = prior_a, prior_b
for heads, n in [(7, 10), (15, 30)]:
    a, b = a + heads, b + (n - heads)

print(f"posterior Beta({a}, {b})")
posterior_mean = a / (a + b)
print("posterior mean", round(posterior_mean, 4))

total_heads, total_n = 22, 40
w = total_n / (prior_a + prior_b + total_n)
prior_mean = prior_a / (prior_a + prior_b)
weighted = w * (total_heads / total_n) + (1 - w) * prior_mean

print("weighted average", round(weighted, 4))
print("identity holds:", abs(posterior_mean - weighted) < 1e-12)

---

A conjugate prior keeps the posterior in the same family as the prior, turning Bayesian updating into arithmetic on parameters. Beta–Binomial, Gamma–Poisson and Normal–Normal are the ones worth memorising, and each makes the shrinkage explicit: the posterior mean is a weighted average of prior and data, weighted by their effective sample sizes. Conjugacy is a convenience — when your real belief or model doesn't fit it, use computation rather than distorting the prior.

Next: reading intervals off a posterior, and how they differ from confidence intervals.