Conjugacy and the Exponential Family

§6.5 showed that the Gaussian is closed under everything you want to do to it. This section asks the obvious follow-up: why, and which other distributions share that?

The book sets out three requirements for a distribution to be usable in machine learning:

1. A closure property when you apply the rules of probability — “applying a particular operation returns an object of the same type.”
2. The parameter count does not grow. “As we collect more data, we do not need more parameters to describe the distribution.”
3. Parameter estimation behaves nicely.

The answer to “which distributions satisfy all three” turns out to be a single family, and the answer is not a coincidence: Pitman, Darmois and Koopman proved in 1935–36 that the exponential family is the only family with finite-dimensional sufficient statistics under repeated independent sampling.

What you’ll learn

• Three more named distributions: Bernoulli (Example 6.8), Binomial (6.9), Beta (6.10), and the Beta’s four qualitative regimes.
• Definition 6.13: a conjugate prior, and why conjugacy turns Bayesian updating into arithmetic on parameters — measured against numerical Bayes.
• Examples 6.11 and 6.12: Beta–Binomial and Beta–Bernoulli conjugacy, and what $\alpha$ and $\beta$ actually count.
• Table 6.2’s catalogue of conjugate pairs.
• §6.6.2: sufficient statistics and the Fisher–Neyman theorem (6.14), demonstrated with two different datasets that no likelihood can tell apart.
• §6.6.3: the exponential family, Equation 6.107 — natural parameters, sufficient statistics, and the log-partition function whose derivative is the mean.
• Examples 6.13 and 6.14: the Gaussian and the Bernoulli written in exponential-family form — and why the sigmoid appears.
• Equation 6.120: every exponential family has a conjugate prior, and Example 6.15 derives the Beta rather than guessing it.

Intuition: what makes a distribution convenient

Bayes’ theorem multiplies a prior by a likelihood. Do that with two arbitrary densities and you get a function with no name, which you can only handle numerically — on a grid, or by sampling.

But sometimes the product lands back in the same family as the prior, with different parameters. Then updating is not integration; it is adding numbers to the parameters. Observe a head: add one to $\alpha$. Observe a tail: add one to $\beta$. That is the whole computation, and it stays two numbers whether you have seen ten observations or ten million.

The second idea makes the first one possible. A sufficient statistic is a summary of the data that the likelihood cannot see past. For a Gaussian, the likelihood only ever touches $n$, $\sum x_i$ and $\sum x_i^2$ — so two completely different datasets sharing those three numbers are indistinguishable to any inference about $\mu$ and $\sigma$. That is why a fixed number of parameters can absorb an unbounded amount of data: the data compresses to a fixed-size summary first.

The exponential family is the class of distributions built to have exactly that property, and Equation 6.107 wears it on its face: the parameters and the data meet only through an inner product $\langle\boldsymbol{\theta},\phi(\mathbf{x})\rangle$. The data enters only as $\phi(\mathbf{x})$, so $\phi$ is the sufficient statistic, by construction.

diagram Diagram mermaid

  flowchart TD
  W["three desiderata:
closure, fixed parameter count, nice estimation"]
  W --> EF["the exponential family, Eq 6.107
p(x|theta) = h(x) exp(theta . phi(x) - A(theta))"]
  EF --> SS["phi(x) is the sufficient statistic
by construction"]
  SS --> FN["Fisher-Neyman, Thm 6.14
p(x|theta) = h(x) g_theta(phi(x))"]
  FN --> FIX["data compresses to a fixed-size summary
so parameters need not grow"]
  EF --> A["A(theta), the log-partition function"]
  A --> MOM["dA/dtheta = E[phi(x)]
differentiating it gives moments"]
  EF --> CONJ["Eq 6.120: every member has a conjugate prior"]
  CONJ --> ARITH["Bayes becomes parameter arithmetic
Def 6.13"]
  ARITH --> EX["Beta-Binomial, Eq 6.104d
alpha + heads, beta + tails"]
  EF --> G["Gaussian, Ex 6.13"]
  EF --> B["Bernoulli, Ex 6.14"]
  B --> SIG["theta = log(mu/(1-mu))
inverse is the SIGMOID"]

The math

Three named distributions

Example 6.8 — Bernoulli. A single binary $x \in \{0,1\}$, governed by $\mu \in [0,1]$:

$$p(x\mid\mu) = \mu^x(1-\mu)^{1-x}, \qquad \mathbb{E}[x] = \mu, \qquad \mathbb{V}[x] = \mu(1-\mu) \tag{6.92--6.94}$$

The exponent trick

The book flags $\mu^x(1-\mu)^{1-x}$ as “a trick that is often used in machine learning textbooks” — putting a Boolean into an exponent to select between two cases without an if. At $x=1$ it is $\mu$; at $x=0$ it is $1-\mu$. The same device appears in the Multinomial and in cross-entropy loss.

Example 6.9 — Binomial. The number of ones in $N$ Bernoulli draws:

$$p(m\mid N,\mu) = \binom{N}{m}\mu^m(1-\mu)^{N-m}, \qquad \mathbb{E}[m] = N\mu, \qquad \mathbb{V}[m] = N\mu(1-\mu) \tag{6.95--6.97}$$

Example 6.10 — Beta. A distribution over a continuous $\mu \in [0,1]$ — so, a distribution over the parameter of a Bernoulli:

$$p(\mu\mid\alpha,\beta) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\mu^{\alpha-1}(1-\mu)^{\beta-1} \tag{6.98}$$ $$\mathbb{E}[\mu] = \frac{\alpha}{\alpha+\beta}, \qquad \mathbb{V}[\mu] = \frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)} \tag{6.99}$$

with the Gamma function $\Gamma(t) := \int_0^\infty x^{t-1}e^{-x}\,\mathrm{d}x$ and its recursion $\Gamma(t+1) = t\,\Gamma(t)$ (Equations 6.100, 6.101) supplying the normaliser.

$\alpha$ pushes mass toward $1$ and $\beta$ toward $0$. Four regimes, all verified below:

Parameters	Shape
$\alpha=\beta=1$	uniform on $[0,1]$
$\alpha,\beta<1$	bimodal, with spikes at $0$ and $1$
$\alpha,\beta>1$	unimodal
$\alpha=\beta>1$	unimodal, symmetric, mode at $\tfrac12$

§6.6.1: conjugate priors

Definition 6.13. A prior is conjugate for the likelihood function if the posterior is of the same form as the prior.

The payoff: “we can algebraically calculate our posterior distribution by updating the parameters of the prior distribution.”

Example 6.11 — Beta–Binomial. With a Binomial likelihood and a Beta prior, observe $h$ heads in $N$ flips:

$$p(\mu\mid x=h,N,\alpha,\beta) \propto p(x\mid N,\mu)\,p(\mu\mid\alpha,\beta) \tag{6.104a}$$ $$\propto \mu^h(1-\mu)^{N-h}\cdot\mu^{\alpha-1}(1-\mu)^{\beta-1} \tag{6.104b}$$ $$= \mu^{h+\alpha-1}(1-\mu)^{(N-h)+\beta-1} \tag{6.104c}$$ $$\propto \operatorname{Beta}(h+\alpha,\ N-h+\beta) \tag{6.104d}$$

Look at what happened: the exponents added. No integral was computed, and the normaliser never had to be worked out — it is whatever the Beta’s normaliser is at the new parameters.

Example 6.12 — Beta–Bernoulli. One observation at a time:

$$p(\theta\mid x,\alpha,\beta) \propto \theta^{x}(1-\theta)^{1-x}\cdot\theta^{\alpha-1}(1-\theta)^{\beta-1} = \theta^{\alpha+x-1}(1-\theta)^{\beta+(1-x)-1} \tag{6.105}$$ $$\propto p\bigl(\theta \mid \alpha+x,\ \beta+(1-x)\bigr) \tag{6.105d}$$

So $\alpha$ counts the ones and $\beta$ counts the zeros. The prior parameters are pseudo-counts, which is exactly why $\operatorname{Beta}(1,1)$ — “one of each” — is the uniform prior.

Table 6.2, the catalogue:

Likelihood	Conjugate prior	Posterior
Bernoulli	Beta	Beta
Binomial	Beta	Beta
Gaussian	Gaussian / inverse Gamma	Gaussian / inverse Gamma
Gaussian	Gaussian / inverse Wishart	Gaussian / inverse Wishart
Multinomial	Dirichlet	Dirichlet

The Gaussian appears twice because the univariate and multivariate cases differ: the inverse Gamma is conjugate for a scalar variance, the inverse Wishart for a covariance matrix. (Equivalently, the book’s margin note: the Gamma is conjugate for the precision, and the Wishart for the precision matrix.)

§6.6.2: sufficient statistics

Fisher’s idea: some statistics “contain all available information that can be inferred from data” about the parameter.

Theorem 6.14 (Fisher–Neyman). $\phi(\mathbf{x})$ are sufficient for $\boldsymbol{\theta}$ if and only if the density factorises as

$$p(\mathbf{x}\mid\boldsymbol{\theta}) = h(\mathbf{x})\,g_{\boldsymbol{\theta}}\bigl(\phi(\mathbf{x})\bigr) \tag{6.106}$$

where $h$ does not depend on $\boldsymbol{\theta}$ and $g_{\boldsymbol{\theta}}$ captures all the dependence on $\boldsymbol{\theta}$ through $\phi(\mathbf{x})$ alone.

The consequence is sharp: if two datasets have the same $\phi$, the likelihood is the same function of $\boldsymbol{\theta}$ for both, so no inference about $\boldsymbol{\theta}$ can distinguish them. Measured below on a Gaussian with $\phi = (n, \sum x_i, \sum x_i^2)$: two datasets that differ by up to $0.26$ per order statistic give log-likelihoods agreeing to $2.3\times10^{-13}$ across the parameter grid.

And then the question that motivates the next subsection: as we see more data, do we need more parameters? “It turns out that the answer is yes in general” — that is non-parametric statistics. The class of distributions with finite-dimensional sufficient statistics is exactly the exponential family.

§6.6.3: the exponential family

An exponential family is a family parameterised by $\boldsymbol{\theta}\in\mathbb{R}^D$ of the form

$$p(\mathbf{x}\mid\boldsymbol{\theta}) = h(\mathbf{x})\exp\bigl(\langle\boldsymbol{\theta},\phi(\mathbf{x})\rangle - A(\boldsymbol{\theta})\bigr) \tag{6.107}$$

where $\phi(\mathbf{x})$ is the vector of sufficient statistics. This is a particular expression of $g_{\boldsymbol{\theta}}(\phi(\mathbf{x}))$ from Theorem 6.14 — the sufficiency is built in.

Ignoring $h$ and $A$ gives the intuition:

$$p(\mathbf{x}\mid\boldsymbol{\theta}) \propto \exp\bigl(\boldsymbol{\theta}^\top\phi(\mathbf{x})\bigr) \tag{6.108}$$

and in this parameterisation $\boldsymbol{\theta}$ are the natural parameters. $A(\boldsymbol{\theta})$ is the normaliser and is called the log-partition function.

A(theta) is not bookkeeping

The log-partition function generates the moments: $\dfrac{\partial A}{\partial\boldsymbol{\theta}} = \mathbb{E}[\phi(\mathbf{x})]$. Measured for the Bernoulli, whose $A(\theta) = \log(1+e^\theta)$: the numerical derivative matches $\sigma(\theta)=\mu$ to about $4\times10^{-11}$ at every $\theta$ tried. So the normalising constant you were tempted to ignore is the thing that hands you the mean — and its second derivative gives the covariance.

Example 6.13 — the Gaussian. Take $\phi(x) = \begin{bmatrix}x\\x^2\end{bmatrix}$. Then $p(x\mid\boldsymbol{\theta}) \propto \exp(\theta_1 x + \theta_2 x^2)$, and setting

$$\boldsymbol{\theta} = \begin{bmatrix}\dfrac{\mu}{\sigma^2} & -\dfrac{1}{2\sigma^2}\end{bmatrix}^\top \tag{6.110}$$

gives

$$p(x\mid\boldsymbol{\theta}) \propto \exp\!\left(\frac{\mu x}{\sigma^2}-\frac{x^2}{2\sigma^2}\right) \propto \exp\!\left(-\frac{1}{2\sigma^2}(x-\mu)^2\right) \tag{6.111}$$

Verified below: the ratio of the two expressions is constant to $1.8\times10^{-15}$. And the map inverts — $\mu = -\theta_1/(2\theta_2)$, $\sigma^2 = -1/(2\theta_2)$.

Example 6.14 — the Bernoulli. Rewrite $\mu^x(1-\mu)^{1-x}$ by taking logs and regrouping:

$$p(x\mid\mu) = \exp\!\left[x\log\frac{\mu}{1-\mu} + \log(1-\mu)\right] \tag{6.113d}$$

which is Equation 6.107 with

$$h(x) = 1, \qquad \theta = \log\frac{\mu}{1-\mu}, \qquad \phi(x) = x, \qquad A(\theta) = -\log(1-\mu) = \log(1+\exp\theta) \tag{6.114--6.117}$$

and the relationship inverts to

$$\mu = \frac{1}{1+\exp(-\theta)} \tag{6.118}$$

This is where the sigmoid comes from

The book’s Remark: the map between the Bernoulli’s natural parameter and $\mu$ is the sigmoid. Note the domains — $\mu\in(0,1)$ but $\theta\in\mathbb{R}$ — so the sigmoid is exactly the function that squeezes an unconstrained real number into a probability. That is why logistic regression puts a linear predictor through a sigmoid: the linear predictor is the natural parameter of a Bernoulli. The same function reappears as a neural-network activation.

Equation 6.120: every member has a conjugate prior. For a likelihood in exponential-family form,

$$p(\boldsymbol{\theta}\mid\boldsymbol{\gamma}) = h_c(\boldsymbol{\theta})\exp\!\left(\left\langle\begin{bmatrix}\gamma_1\\\gamma_2\end{bmatrix},\begin{bmatrix}\boldsymbol{\theta}\\-A(\boldsymbol{\theta})\end{bmatrix}\right\rangle - A_c(\boldsymbol{\gamma})\right) \tag{6.120}$$

with $\boldsymbol{\gamma}$ of dimension $\dim(\boldsymbol{\theta})+1$. So conjugate priors are not lucky finds — they are produced by a formula.

Example 6.15 applies it to the Bernoulli. Substituting $\theta = \log\frac{\mu}{1-\mu}$, $A(\theta) = -\log(1-\mu)$ and $\boldsymbol{\gamma} = [\alpha,\ \beta+\alpha]^\top$, the inner product collapses:

$$\alpha\log\frac{\mu}{1-\mu} + (\beta+\alpha)\log(1-\mu) = \alpha\log\mu + \beta\log(1-\mu)$$

so the exponential factor alone is $\mu^\alpha(1-\mu)^\beta$ — verified to $1.0\times10^{-17}$. Equation 6.124 wants $\mu^{\alpha-1}(1-\mu)^{\beta-1}$, so $h_c$ must supply both $-1$ exponents, giving $h_c(\mu) = \dfrac{1}{\mu(1-\mu)}$ and

$$p(\mu\mid\alpha,\beta) \propto \mu^{\alpha-1}(1-\mu)^{\beta-1} \tag{6.124}$$

which is the Beta. The Beta was not guessed — it is what Equation 6.120 produces for a Bernoulli likelihood.

Worked example by hand

Example 6.11 with numbers

Prior $\operatorname{Beta}(2,5)$ — mean $2/7 = 0.2857$, so mildly pessimistic about $\mu$. Observe $h=26$ heads in $N=40$ flips.

By Equation 6.104d, no calculus:

$$\alpha' = h+\alpha = 26+2 = 28, \qquad \beta' = N-h+\beta = 14+5 = 19$$

so the posterior is $\operatorname{Beta}(28,19)$ with mean $28/47 = 0.5957446809$.

By brute force, to check: put a $20\,001$-point grid on $[0,1]$, multiply the Binomial likelihood by the Beta prior, normalise numerically. Measured against the closed form:

	value
worst absolute density gap	$5.4\times10^{-14}$
relative to the peak	$9.8\times10^{-15}$
posterior mean, numerical	$0.5957446809$
posterior mean, closed form	$0.5957446809$

The same distribution, to every printed digit. Conjugacy is not an approximation that happens to be close — it is an algebraic identity.

Note also where the posterior mean sits: $0.5957$, between the prior mean $0.2857$ and the raw data proportion $26/40 = 0.65$. With $\alpha+\beta=7$ pseudo-counts against $40$ real observations, the data dominates but has not erased the prior.

The parameter count, measured

Start from $\operatorname{Beta}(1,1)$ and stream observations from a coin with $\mu = 0.62$:

observations	$\alpha$	$\beta$	numbers stored	posterior mean	posterior sd
$0$	$1$	$1$	2	$0.500000$	$0.288675$
$10$	$5$	$7$	2	$0.416667$	$0.136735$
$1\,000$	$616$	$386$	2	$0.614770$	$0.015366$
$100\,000$	$61\,789$	$38\,213$	2	$0.617878$	$0.001537$

Two numbers, at every scale. That is desideratum 2 satisfied exactly. A non-conjugate prior would need a grid, and at $200$ points per axis a grid costs $200^D$: $40\,000$ at $D=2$, $8$ million at $D=3$, $1.0\times10^{23}$ at $D=10$.

Sufficiency, made concrete

Build two datasets of $400$ points each with the same $n$, $\sum x_i$ and $\sum x_i^2$ — one drawn from $\mathcal{N}(2,1.5^2)$, the other from $\mathcal{N}(-3,4^2)$ then rescaled to match.

	dataset 1	dataset 2	gap
$\sum x_i$	$820.48115864$	$820.48115864$	$0.0\text{e}{+}00$
$\sum x_i^2$	$2603.52559396$	$2603.52559396$	$0.0\text{e}{+}00$
worst sorted difference	—	—	$0.260077$

They are genuinely different data. Now evaluate the Gaussian log-likelihood at several $(\mu,\sigma)$:

$\mu$	$\sigma$	dataset 1	dataset 2	gap
$0.0$	$1.5$	$-1108.32269963$	$-1108.32269963$	$0.0\text{e}{+}00$
$1.0$	$1.5$	$-832.55329579$	$-832.55329579$	$0.0\text{e}{+}00$
$2.0$	$1.5$	$-734.56166973$	$-734.56166973$	$1.1\text{e}{-}13$
$3.0$	$1.5$	$-814.34782144$	$-814.34782144$	$0.0\text{e}{+}00$

Worst gap over twelve parameter settings: $2.3\times10^{-13}$. The likelihood cannot see the difference. Everything inferable about $(\mu,\sigma)$ lives in those two sums, which is exactly what Theorem 6.14 asserts.

See it move

sketch Conjugacy is parameter arithmetic p5.js Pause Restart

Flip coins and watch the posterior update. Every observation moves alpha or beta by exactly one, and the whole posterior is those two numbers — the curve is drawn from them, never stored. The counter shows how many numbers a grid posterior would need for the same job.

sketch The natural parameter, and the sigmoid p5.js Pause Restart

Drag mu and watch the Bernoulli's natural parameter theta. The map is unbounded in theta and squeezed in mu, and its inverse is exactly the sigmoid — which is why a linear predictor pushed through a sigmoid is a Bernoulli likelihood. The lower panel shows the log-partition function and the fact that its slope is mu.

From scratch

conjugacy_and_exponential_family.py

"""Section 6.6 — conjugacy, sufficient statistics and the exponential family."""
 
import math
 
import numpy as np
 
np.set_printoptions(precision=6, suppress=True, linewidth=150)
 
print("########## three_named_distributions")
rng = np.random.default_rng(6)
# Example 6.8, Eq 6.92 to 6.94.
mu = 0.3
x = (rng.random(4_000_000) < mu).astype(float)
print(f"  Bernoulli(mu = {mu})")
print(f"    Eq 6.93  E[x] = mu          = {mu:.6f}   measured {x.mean():.6f}")
print(f"    Eq 6.94  V[x] = mu(1 - mu)  = {mu*(1-mu):.6f}   measured {x.var():.6f}")
 
# Example 6.9, Eq 6.95 to 6.97.
N_, mu2 = 15, 0.4
m = rng.binomial(N_, mu2, 2_000_000)
print(f"  Binomial(N = {N_}, mu = {mu2})")
print(f"    Eq 6.96  E[m] = N mu           = {N_*mu2:.6f}   measured {m.mean():.6f}")
print(f"    Eq 6.97  V[m] = N mu(1 - mu)   = {N_*mu2*(1-mu2):.6f}   measured {m.var():.6f}")
# Eq 6.95 against a direct count.
counts = np.bincount(m, minlength=N_ + 1) / m.size
pmf = np.array([math.comb(N_, k) * mu2**k * (1-mu2)**(N_-k) for k in range(N_+1)])
print(f"    Eq 6.95 pmf worst error over all {N_+1} outcomes: "
      f"{np.abs(counts - pmf).max():.6f}   sums to {pmf.sum():.10f}")
 
# Example 6.10, Eq 6.98 to 6.99.
al, be = 4.0, 10.0
b = rng.beta(al, be, 2_000_000)
e_th = al / (al + be)
v_th = al * be / ((al + be) ** 2 * (al + be + 1))
print(f"  Beta(alpha = {al}, beta = {be})")
print(f"    Eq 6.99  E[mu] = a/(a+b)   = {e_th:.6f}   measured {b.mean():.6f}")
print(f"    Eq 6.99  V[mu]             = {v_th:.6f}   measured {b.var():.6f}")
# Eq 6.100 / 6.101: the Gamma recursion, which normalises the Beta.
print(f"    Eq 6.101  Gamma(t+1) = t Gamma(t):  "
      f"Gamma(5.5) = {math.gamma(5.5):.6f}  vs 4.5*Gamma(4.5) = {4.5*math.gamma(4.5):.6f}")
 
print()
print("########## beta_special_cases")
# The book's four cases, checked by counting modes of the density.
grid = np.linspace(1e-6, 1 - 1e-6, 20001)
 
 
def beta_pdf(a, bb):
    logc = math.lgamma(a + bb) - math.lgamma(a) - math.lgamma(bb)
    return np.exp(logc + (a - 1) * np.log(grid) + (bb - 1) * np.log1p(-grid))
 
 
def modes_of(v):
    inner = [i for i in range(1, len(v) - 1) if v[i] > v[i-1] and v[i] > v[i+1]]
    edges = []
    if v[0] > v[1]:
        edges.append("0")
    if v[-1] > v[-2]:
        edges.append("1")
    return len(inner), edges
 
 
print(f"  {'alpha, beta':>14}  {'interior modes':>15}  {'spikes at':>10}  verdict")
for a, bb, label in ((1.0, 1.0, "uniform"), (0.5, 0.5, "bimodal, spikes at 0 and 1"),
                     (2.0, 0.3, "mass toward 1"), (4.0, 10.0, "unimodal"),
                     (5.0, 5.0, "unimodal, symmetric, mode at 1/2"),
                     (5.0, 1.0, "mass toward 1")):
    v = beta_pdf(a, bb)
    ni, ed = modes_of(v)
    print(f"  {f'{a}, {bb}':>14}  {ni:>15}  {','.join(ed) or '-':>10}  {label}")
v11 = beta_pdf(1.0, 1.0)
print(f"  alpha = beta = 1 really is uniform: density range "
      f"[{v11.min():.10f}, {v11.max():.10f}]")
v55 = beta_pdf(5.0, 5.0)
print(f"  alpha = beta = 5 mode at {grid[int(np.argmax(v55))]:.6f}  (theory 0.5)")
 
print()
print("########## beta_binomial_conjugacy")
# Example 6.11, Eq 6.102 to 6.104d, checked against numerical Bayes.
a0, b0, N_obs, h = 2.0, 5.0, 40, 26
print(f"  prior Beta({a0}, {b0}), observe h = {h} heads in N = {N_obs} flips")
a1, b1 = h + a0, N_obs - h + b0
print(f"  Eq 6.104d  posterior = Beta(h + alpha, N - h + beta) = Beta({a1}, {b1})")
# Now do it the hard way: grid, multiply, normalise.
lik = np.array([math.comb(N_obs, h)]) * grid**h * (1 - grid)**(N_obs - h)
prior = beta_pdf(a0, b0)
post_num = lik * prior
post_num /= np.trapezoid(post_num, grid)
post_closed = beta_pdf(a1, b1)
print(f"  numerical posterior vs the closed form:")
print(f"    worst absolute density gap  {np.abs(post_num - post_closed).max():.3e}")
print(f"    relative to the peak        "
      f"{np.abs(post_num - post_closed).max()/post_closed.max():.3e}")
print(f"    posterior mean, numerical  {np.trapezoid(grid*post_num, grid):.10f}")
print(f"    posterior mean, closed     {a1/(a1+b1):.10f}")
print("  conjugacy is not an approximation: it is the same distribution, reached by")
print("  adding the counts to the prior's parameters.")
 
print()
print("########## example_6_12_bernoulli")
# Eq 6.105a to 6.105d: one observation moves alpha or beta by exactly one.
print(f"  Beta-Bernoulli: prior Beta({a0}, {b0}), a single observation x")
for xv in (1, 0):
    print(f"    x = {xv}  ->  Beta({a0 + xv}, {b0 + (1 - xv)})   "
          f"(Eq 6.105d: alpha + x, beta + 1 - x)")
print("  so alpha counts the ones and beta counts the zeros. The prior parameters")
print("  are PSEUDO-COUNTS, which is why Beta(1,1) is 'one of each', i.e. uniform.")
 
print()
print("########## the_parameter_count_does_not_grow")
# The book's desideratum 2, made concrete.
rng2 = np.random.default_rng(21)
truth = 0.62
data = (rng2.random(100_000) < truth).astype(int)
a, bb = 1.0, 1.0
print(f"  sequential updating from Beta(1,1), truth = {truth}")
print(f"  {'observations':>13}  {'alpha':>10}  {'beta':>10}  {'numbers stored':>15}"
      f"  {'post. mean':>11}  {'post. sd':>9}")
seen = 0
for target in (0, 1, 10, 100, 1_000, 10_000, 100_000):
    while seen < target:
        a += data[seen]
        bb += 1 - data[seen]
        seen += 1
    m_ = a / (a + bb)
    s_ = math.sqrt(a * bb / ((a + bb) ** 2 * (a + bb + 1)))
    print(f"  {seen:>13,}  {a:>10.0f}  {bb:>10.0f}  {2:>15}  {m_:>11.6f}  {s_:>9.6f}")
print("  two numbers, forever. A non-conjugate prior would need a grid, and to keep")
print("  the same resolution in D dimensions that grid costs:")
for D in (1, 2, 3, 5, 10):
    print(f"    D = {D:>2}: {200**D:,} grid points at 200 per axis")
 
print()
print("########## sufficient_statistics")
# Theorem 6.14. For a Gaussian, (sum x, sum x^2) is sufficient: two datasets
# sharing them have IDENTICAL likelihood at every parameter value.
rng3 = np.random.default_rng(4)
d1 = rng3.normal(2.0, 1.5, 400)
# Build a completely different dataset with the same n, sum and sum of squares.
d2 = rng3.normal(-3.0, 4.0, 400)
d2 = (d2 - d2.mean()) / d2.std() * d1.std() + d1.mean()
print(f"  two datasets of {d1.size} points each")
print(f"    sum x   : {d1.sum():.8f}  and  {d2.sum():.8f}   gap {abs(d1.sum()-d2.sum()):.1e}")
print(f"    sum x^2 : {np.sum(d1**2):.8f}  and  {np.sum(d2**2):.8f}   "
      f"gap {abs(np.sum(d1**2)-np.sum(d2**2)):.1e}")
print(f"    but the data differs: worst |sorted difference| "
      f"{np.abs(np.sort(d1)-np.sort(d2)).max():.6f}")
 
 
def gauss_loglik(d, m_, s_):
    return float(np.sum(-0.5 * np.log(2 * np.pi * s_**2) - (d - m_)**2 / (2 * s_**2)))
 
 
print(f"  {'mu':>6}  {'sigma':>6}  {'loglik dataset 1':>18}  {'loglik dataset 2':>18}  {'gap':>9}")
worst = 0.0
for m_ in (0.0, 1.0, 2.0, 3.0):
    for s_ in (1.0, 1.5, 3.0):
        l1, l2 = gauss_loglik(d1, m_, s_), gauss_loglik(d2, m_, s_)
        worst = max(worst, abs(l1 - l2))
        if s_ == 1.5:
            print(f"  {m_:>6.1f}  {s_:>6.1f}  {l1:>18.8f}  {l2:>18.8f}  {abs(l1-l2):>9.1e}")
print(f"  worst gap over all 12 parameter settings tried: {worst:.1e}")
print("  the likelihood cannot tell these datasets apart. That is what Theorem")
print("  6.14 means: everything inferable about (mu, sigma) is in those two sums.")
 
print()
print("########## gaussian_as_exponential_family")
# Example 6.13, Eq 6.109 to 6.111.
mu_g, sd_g = 1.3, 0.8
th = np.array([mu_g / sd_g**2, -1.0 / (2 * sd_g**2)])        # Eq 6.110
print(f"  N({mu_g}, {sd_g**2:.4f})  ->  natural parameters theta = {th}")
xs = np.linspace(-3, 6, 9)
phi = np.stack([xs, xs**2], axis=1)                           # sufficient statistics
unnorm = np.exp(phi @ th)                                     # Eq 6.109
true = np.exp(-0.5 * (xs - mu_g)**2 / sd_g**2)                # Eq 6.111, unnormalised
ratio = unnorm / true
print(f"  Eq 6.111: exp(theta . phi(x)) / exp(-(x-mu)^2/2sigma^2) should be constant")
print(f"    ratios: {ratio}")
print(f"    spread of the ratio: {ratio.max()/ratio.min() - 1:.2e}   (constant, as claimed)")
print(f"  and inverting: mu = -theta1/(2 theta2) = {-th[0]/(2*th[1]):.6f}")
print(f"                 sigma^2 = -1/(2 theta2) = {-1/(2*th[1]):.6f}")
 
print()
print("########## bernoulli_as_exponential_family")
# Example 6.14, Eq 6.112 to 6.118.
print(f"  {'mu':>6}  {'theta = log(mu/(1-mu))':>23}  {'A(theta) = log(1+e^th)':>23}"
      f"  {'-log(1-mu)':>12}  {'sigmoid(theta)':>15}")
for m_ in (0.1, 0.3, 0.5, 0.7, 0.9):
    t_ = math.log(m_ / (1 - m_))                              # Eq 6.115
    A_ = math.log1p(math.exp(t_))                             # Eq 6.117 right side
    A2 = -math.log(1 - m_)                                    # Eq 6.117 left side
    inv = 1 / (1 + math.exp(-t_))                             # Eq 6.118
    print(f"  {m_:>6.1f}  {t_:>23.10f}  {A_:>23.10f}  {A2:>12.10f}  {inv:>15.10f}")
print("  Eq 6.117's two forms agree, and Eq 6.118 inverts the map: the link between")
print("  the natural parameter and mu IS the sigmoid, which is why logistic")
print("  regression and a Bernoulli likelihood are the same statement.")
print()
print("  and the log-partition function generates the moments -- dA/dtheta = mu:")
print(f"  {'theta':>8}  {'dA/dtheta (numeric)':>20}  {'sigmoid(theta)':>15}  {'gap':>9}")
for t_ in (-2.0, -0.5, 0.0, 0.5, 2.0):
    hh = 1e-6
    dA = (math.log1p(math.exp(t_ + hh)) - math.log1p(math.exp(t_ - hh))) / (2 * hh)
    sg = 1 / (1 + math.exp(-t_))
    print(f"  {t_:>8.1f}  {dA:>20.10f}  {sg:>15.10f}  {abs(dA-sg):>9.1e}")
print("  so A(theta) is not bookkeeping: differentiating it hands you E[phi(x)].")
 
print()
print("########## example_6_15_derive_the_beta")
# Eq 6.119 to 6.124: the canonical conjugate of the Bernoulli IS the Beta.
print("  Eq 6.120's canonical conjugate for the Bernoulli, with gamma = [alpha, beta+alpha],")
print("  simplifies through Eq 6.123 to Eq 6.124:")
print("    p(mu | alpha, beta) proportional to mu^(alpha-1) (1-mu)^(beta-1)")
aa, bbb = 3.0, 5.0
# The inner product of Eq 6.120, with theta = log(mu/(1-mu)), A(theta) =
# -log(1-mu) and gamma = [alpha, beta+alpha]:
#     alpha log(mu/(1-mu)) + (beta+alpha) log(1-mu)
#   = alpha log mu - alpha log(1-mu) + (beta+alpha) log(1-mu)
#   = alpha log mu + beta log(1-mu)
# so the exponential factor alone is mu^alpha (1-mu)^beta.
core = np.exp(aa * np.log(grid / (1 - grid)) + (bbb + aa) * np.log1p(-grid))
simplified = grid ** aa * (1 - grid) ** bbb
print(f"  the exponential factor of Eq 6.122, at alpha = {aa}, beta = {bbb}:")
print(f"    it equals mu^alpha (1-mu)^beta to "
      f"{np.abs(core - simplified).max():.2e}")
print()
print("  Eq 6.124 wants mu^(alpha-1) (1-mu)^(beta-1), so h_c has to supply BOTH")
print("  minus-one exponents. Only one candidate does:")
rhs = grid ** (aa - 1) * (1 - grid) ** (bbb - 1)                     # Eq 6.124
for label, hc in (("mu/(1-mu)", grid / (1 - grid)),
                  ("1/(mu(1-mu))", 1.0 / (grid * (1 - grid)))):
    r = (hc * core) / rhs
    flat = r.max() / r.min() - 1
    verdict = "CONSTANT -> this is h_c" if flat < 1e-9 else "not constant"
    print(f"    h_c = {label:<14} ratio spread {flat:.2e}   {verdict}")
print("  (the PDF renders h_c as a two-line fraction that extracts ambiguously;")
print("   the algebra settles it.)")
norm_rhs = rhs / np.trapezoid(rhs, grid)
print(f"  and normalising Eq 6.124 gives Beta({aa}, {bbb}) to "
      f"{np.abs(norm_rhs - beta_pdf(aa, bbb)).max():.2e}")
print("  so the Beta was not guessed: it is what Eq 6.120 produces for a Bernoulli.")

output

########## three_named_distributions
  Bernoulli(mu = 0.3)
    Eq 6.93  E[x] = mu          = 0.300000   measured 0.299841
    Eq 6.94  V[x] = mu(1 - mu)  = 0.210000   measured 0.209936
  Binomial(N = 15, mu = 0.4)
    Eq 6.96  E[m] = N mu           = 6.000000   measured 5.997410
    Eq 6.97  V[m] = N mu(1 - mu)   = 3.600000   measured 3.596868
    Eq 6.95 pmf worst error over all 16 outcomes: 0.000235   sums to 1.0000000000
  Beta(alpha = 4.0, beta = 10.0)
    Eq 6.99  E[mu] = a/(a+b)   = 0.285714   measured 0.285665
    Eq 6.99  V[mu]             = 0.013605   measured 0.013591
    Eq 6.101  Gamma(t+1) = t Gamma(t):  Gamma(5.5) = 52.342778  vs 4.5*Gamma(4.5) = 52.342778
 
########## beta_special_cases
     alpha, beta   interior modes   spikes at  verdict
        1.0, 1.0                0           -  uniform
        0.5, 0.5                0         0,1  bimodal, spikes at 0 and 1
        2.0, 0.3                0           1  mass toward 1
       4.0, 10.0                1           -  unimodal
        5.0, 5.0                1           -  unimodal, symmetric, mode at 1/2
        5.0, 1.0                0           1  mass toward 1
  alpha = beta = 1 really is uniform: density range [1.0000000000, 1.0000000000]
  alpha = beta = 5 mode at 0.500000  (theory 0.5)
 
########## beta_binomial_conjugacy
  prior Beta(2.0, 5.0), observe h = 26 heads in N = 40 flips
  Eq 6.104d  posterior = Beta(h + alpha, N - h + beta) = Beta(28.0, 19.0)
  numerical posterior vs the closed form:
    worst absolute density gap  5.418e-14
    relative to the peak        9.759e-15
    posterior mean, numerical  0.5957446809
    posterior mean, closed     0.5957446809
  conjugacy is not an approximation: it is the same distribution, reached by
  adding the counts to the prior's parameters.
 
########## example_6_12_bernoulli
  Beta-Bernoulli: prior Beta(2.0, 5.0), a single observation x
    x = 1  ->  Beta(3.0, 5.0)   (Eq 6.105d: alpha + x, beta + 1 - x)
    x = 0  ->  Beta(2.0, 6.0)   (Eq 6.105d: alpha + x, beta + 1 - x)
  so alpha counts the ones and beta counts the zeros. The prior parameters
  are PSEUDO-COUNTS, which is why Beta(1,1) is 'one of each', i.e. uniform.
 
########## the_parameter_count_does_not_grow
  sequential updating from Beta(1,1), truth = 0.62
   observations       alpha        beta   numbers stored   post. mean   post. sd
              0           1           1                2     0.500000   0.288675
              1           1           2                2     0.333333   0.235702
             10           5           7                2     0.416667   0.136735
            100          57          45                2     0.558824   0.048924
          1,000         616         386                2     0.614770   0.015366
         10,000        6157        3845                2     0.615577   0.004864
        100,000       61789       38213                2     0.617878   0.001537
  two numbers, forever. A non-conjugate prior would need a grid, and to keep
  the same resolution in D dimensions that grid costs:
    D =  1: 200 grid points at 200 per axis
    D =  2: 40,000 grid points at 200 per axis
    D =  3: 8,000,000 grid points at 200 per axis
    D =  5: 320,000,000,000 grid points at 200 per axis
    D = 10: 102,400,000,000,000,000,000,000 grid points at 200 per axis
 
########## sufficient_statistics
  two datasets of 400 points each
    sum x   : 820.48115864  and  820.48115864   gap 0.0e+00
    sum x^2 : 2603.52559396  and  2603.52559396   gap 0.0e+00
    but the data differs: worst |sorted difference| 0.260077
      mu   sigma    loglik dataset 1    loglik dataset 2        gap
     0.0     1.5      -1108.32269963      -1108.32269963    0.0e+00
     1.0     1.5       -832.55329579       -832.55329579    0.0e+00
     2.0     1.5       -734.56166973       -734.56166973    1.1e-13
     3.0     1.5       -814.34782144       -814.34782144    0.0e+00
  worst gap over all 12 parameter settings tried: 2.3e-13
  the likelihood cannot tell these datasets apart. That is what Theorem
  6.14 means: everything inferable about (mu, sigma) is in those two sums.
 
########## gaussian_as_exponential_family
  N(1.3, 0.6400)  ->  natural parameters theta = [ 2.03125 -0.78125]
  Eq 6.111: exp(theta . phi(x)) / exp(-(x-mu)^2/2sigma^2) should be constant
    ratios: [3.744591 3.744591 3.744591 3.744591 3.744591 3.744591 3.744591 3.744591 3.744591]
    spread of the ratio: 1.78e-15   (constant, as claimed)
  and inverting: mu = -theta1/(2 theta2) = 1.300000
                 sigma^2 = -1/(2 theta2) = 0.640000
 
########## bernoulli_as_exponential_family
      mu   theta = log(mu/(1-mu))   A(theta) = log(1+e^th)    -log(1-mu)   sigmoid(theta)
     0.1            -2.1972245773             0.1053605157  0.1053605157     0.1000000000
     0.3            -0.8472978604             0.3566749439  0.3566749439     0.3000000000
     0.5             0.0000000000             0.6931471806  0.6931471806     0.5000000000
     0.7             0.8472978604             1.2039728043  1.2039728043     0.7000000000
     0.9             2.1972245773             2.3025850930  2.3025850930     0.9000000000
  Eq 6.117's two forms agree, and Eq 6.118 inverts the map: the link between
  the natural parameter and mu IS the sigmoid, which is why logistic
  regression and a Bernoulli likelihood are the same statement.
 
  and the log-partition function generates the moments -- dA/dtheta = mu:
     theta   dA/dtheta (numeric)   sigmoid(theta)        gap
      -2.0          0.1192029220     0.1192029220    2.8e-12
      -0.5          0.3775406688     0.3775406688    1.3e-11
       0.0          0.5000000000     0.5000000000    4.1e-11
       0.5          0.6224593312     0.6224593312    4.4e-11
       2.0          0.8807970779     0.8807970780    3.8e-11
  so A(theta) is not bookkeeping: differentiating it hands you E[phi(x)].
 
########## example_6_15_derive_the_beta
  Eq 6.120's canonical conjugate for the Bernoulli, with gamma = [alpha, beta+alpha],
  simplifies through Eq 6.123 to Eq 6.124:
    p(mu | alpha, beta) proportional to mu^(alpha-1) (1-mu)^(beta-1)
  the exponential factor of Eq 6.122, at alpha = 3.0, beta = 5.0:
    it equals mu^alpha (1-mu)^beta to 1.04e-17
 
  Eq 6.124 wants mu^(alpha-1) (1-mu)^(beta-1), so h_c has to supply BOTH
  minus-one exponents. Only one candidate does:
    h_c = mu/(1-mu)      ratio spread 1.00e+12   not constant
    h_c = 1/(mu(1-mu))   ratio spread 1.71e-14   CONSTANT -> this is h_c
  (the PDF renders h_c as a two-line fraction that extracts ambiguously;
   the algebra settles it.)
  and normalising Eq 6.124 gives Beta(3.0, 5.0) to 1.02e-13
  so the Beta was not guessed: it is what Eq 6.120 produces for a Bernoulli.

Five things worth stopping on.

Conjugacy is exact. The closed-form $\operatorname{Beta}(28,19)$ matches grid-based numerical Bayes to a relative $9.8\times10^{-15}$, and both posterior means read $0.5957446809$.

The Beta’s four regimes are real. $\alpha=\beta=1$ gives a density flat to ten decimal places; $\alpha=\beta=0.5$ has zero interior modes and spikes at both ends; $\alpha=\beta=5$ has one mode, measured at exactly $0.500000$.

Sufficiency is total, not approximate. Two datasets differing by up to $0.26$ per order statistic produce log-likelihoods agreeing to $2.3\times10^{-13}$ across twelve parameter settings.

The exponential-family rewrite is an identity. The Gaussian’s ratio check is constant to $1.8\times10^{-15}$; the Bernoulli’s two forms of $A(\theta)$ agree exactly and the sigmoid inverts the link exactly.

$A(\theta)$ generates the mean. $\mathrm{d}A/\mathrm{d}\theta$ measured numerically matches $\sigma(\theta)$ to about $4\times10^{-11}$ at every $\theta$.

And one thing I had to work out rather than read off: the book’s $h_c$ in Equation 6.122 renders in the PDF as a two-line fraction that extracts ambiguously. The algebra settles it — the exponential factor simplifies to $\mu^\alpha(1-\mu)^\beta$, so $h_c$ must supply both $-1$ exponents. Testing both candidates: $\mu/(1-\mu)$ gives a ratio spread of $1.0\times10^{12}$ (not constant), and $1/(\mu(1-\mu))$ gives $1.7\times10^{-14}$ — constant, so that is the one.

On real data

figure The book's Figures 6.10 and 6.11 matplotlib

Left: the Binomial at the three values of mu the book plots, with dotted lines at N mu = 1.5, 6.0 and 11.25 — Equation 6.96. Right: the Beta at the parameter pairs the book lists, covering all four of its regimes. alpha equals beta equals one is flat to ten decimal places; both below one spikes at zero and one; both above one is unimodal; and equal parameters above one is symmetric with its mode measured at exactly 0.5.

figure What conjugacy actually buys matplotlib

Left: every posterior is a Beta, and the update is arithmetic — add the heads to alpha, the tails to beta. No integral is computed at any point. Right: the second desideratum made concrete. The conjugate posterior stays at dim(gamma) numbers however much data arrives; a grid posterior at the same resolution costs 200 to the power D, which is 40 thousand in two dimensions and 1.0e+23 in ten.

figure Theorem 6.14: two datasets, one likelihood matplotlib

The two histograms are different data — up to 0.26 apart per order statistic. They were constructed to share n, the sum and the sum of squares, which for a Gaussian are the sufficient statistics. The right panel overlays their log-likelihood surfaces over mu and sigma: solid and dashed coincide everywhere, worst gap 2.3e-13. No inference about mu or sigma can distinguish these datasets, which is exactly what sufficiency means — and why a fixed parameter count can absorb unbounded data.

figure Example 6.14: where the sigmoid comes from matplotlib

Left: the Bernoulli's natural parameter theta = log(mu/(1-mu)) runs over the whole real line while mu is confined to the unit interval. Right: the inverse is the sigmoid, which is precisely the function that turns an unconstrained linear predictor into a probability — logistic regression is a Bernoulli likelihood in natural parameters. Below: the log-partition function A(theta) = log(1 + e^theta), whose derivative is the mean. The measured derivative matches the sigmoid to about 4e-11, so the normalising constant is what hands you the moments.

Reading the plot

From the named-distributions figure. The right panel’s U-shaped red curve is the one worth noticing. A Beta with both parameters below $1$ is bimodal with spikes at the endpoints — a “prior over a probability” that believes the coin is either almost always heads or almost always tails, and is nearly certain it is not fair. That is a legitimate and occasionally useful belief, and it is what $\operatorname{Beta}(0.5, 0.5)$ encodes.

From the conjugacy figure. The right panel is the whole argument for conjugacy. The green line does not move. Everything else in Bayesian computation — MCMC, variational inference, all of Chapter 11 — exists because most models are not on that green line.

From the sufficiency figure. Look at the two histograms, then at the overlaid contours. Your eye can trivially distinguish the datasets; the likelihood cannot. Anything that only ever consults the likelihood — maximum likelihood, Bayesian posteriors, likelihood-ratio tests — is blind to every difference between them.

From the link figure. The lower panel makes a point that is easy to skip: the tangent’s slope is $\mu$. So if you have $A(\boldsymbol{\theta})$ in closed form, you have every moment by differentiation, with no integration — which is the computational reason exponential families are the family of choice.

From the book

Mathematics for Machine Learning, §6.6. The section opens with the three desiderata (closure, non-growing parameter count, well-behaved estimation), then gives Example 6.8 (Bernoulli, Eq 6.92–6.94, with the exponent-trick Remark), Example 6.9 (Binomial, 6.95–6.97, Figure 6.10) and Example 6.10 (Beta, 6.98–6.101, Figure 6.11, with the four special cases). §6.6.1 gives Definition 6.13, Example 6.11 (Beta–Binomial, 6.102–6.104d), Table 6.2 and Example 6.12 (Beta–Bernoulli, 6.105). §6.6.2 introduces sufficient statistics and the Fisher–Neyman theorem (6.14), Equation 6.106, and notes that the exponential family is the class with finite-dimensional sufficient statistics. §6.6.3 gives the three levels of abstraction, the Pitman–Darmois–Koopman anecdote, Equation 6.107, natural parameters (6.108), the log-partition function, Example 6.13 (Gaussian, 6.109–6.111), Example 6.14 (Bernoulli, 6.112–6.118, with the sigmoid Remark), the general conjugate prior Equation 6.120, and Example 6.15 deriving the Beta (6.121–6.124).

Pitfalls

Conjugacy is a property of the PAIR, not of the prior

“Beta is a conjugate prior” is incomplete — Beta is conjugate for a Bernoulli or Binomial likelihood. Put a Beta prior on the mean of a Gaussian and nothing cancels. Table 6.2 is a table of pairs for that reason, and the Gaussian appears twice because the univariate and multivariate cases need different partners (inverse Gamma against inverse Wishart).

Sufficiency means the likelihood is blind, not that the data is

Two datasets sharing $(n, \sum x_i, \sum x_i^2)$ give identical Gaussian log-likelihoods — measured to $2.3\times10^{-13}$ — while visibly differing. Everything downstream of the likelihood inherits that blindness. If you care about skew, multimodality or outliers, a Gaussian model cannot see them no matter how much data you supply, because they are not in the sufficient statistics.

A conjugate prior is a modelling choice, not a neutral one

$\alpha$ and $\beta$ are pseudo-counts, so choosing $\operatorname{Beta}(2,5)$ asserts you have already seen seven observations, two of them heads. That is a real claim. Measured on the worked example: with $40$ real observations the posterior mean is $0.5957$, pulled down from the data’s own $0.65$ toward the prior mean $0.2857$. Convenient priors are not uninformative priors.

Do not confuse the natural parameter with the usual one

For a Bernoulli, $\mu \in (0,1)$ but $\theta \in \mathbb{R}$, and they are related by the sigmoid, not by equality. Optimising in $\theta$ is unconstrained; optimising in $\mu$ needs a box constraint. This is exactly why logistic regression parameterises the linear predictor and not the probability, and why a model that predicts $\mu$ directly needs a squashing function bolted on.

The log-partition function is not a constant to be dropped

$A(\boldsymbol{\theta})$ depends on $\boldsymbol{\theta}$, so it does not vanish when you differentiate a log-likelihood with respect to the parameters — and its derivative is $\mathbb{E}[\phi(\mathbf{x})]$. Dropping it is how gradients of exponential-family likelihoods get silently wrong; keeping it is how maximum likelihood turns into “match the empirical sufficient statistics to their expectations”.

Compare

	Bernoulli	Binomial	Beta
Support	$\{0,1\}$	$\{0,\ldots,N\}$	$[0,1]$, continuous
Parameters	$\mu$	$N$, $\mu$	$\alpha>0$, $\beta>0$
Mean	$\mu$	$N\mu$	$\alpha/(\alpha+\beta)$
Variance	$\mu(1-\mu)$	$N\mu(1-\mu)$	Eq 6.99
Role here	likelihood	likelihood	prior on $\mu$

	conjugate posterior	grid posterior
Update	add to the parameters	multiply and renormalise
Numbers stored	$\dim(\boldsymbol{\gamma})$, fixed	$200^D$
Exact?	yes, algebraically	to grid resolution
At $D=10$	a handful	$1.0\times10^{23}$
Available for	exponential families	anything

Piece of Equation 6.107	Name	What it does
$\phi(\mathbf{x})$	sufficient statistics	the only route from data to parameters
$\boldsymbol{\theta}$	natural parameters	unconstrained, unlike the usual parameters
$A(\boldsymbol{\theta})$	log-partition function	normalises; its derivative is $\mathbb{E}[\phi]$
$h(\mathbf{x})$	base measure	parameter-free; absorbable into $\phi$

pch.quizTag Check your understanding

1. What does conjugacy buy you, concretely?

   A A more accurate posterior B The posterior stays in the prior's family, so updating is arithmetic on parameters rather than an integral — and the number of stored parameters never grows. Measured: still two numbers after 100000 observations, against 200^D for a grid C A faster sampler D An uninformative prior
   pch.quizShowAnswer

   B — The posterior stays in the prior's family, so updating is arithmetic on parameters rather than an integral — and the number of stored parameters never grows. Measured: still two numbers after 100000 observations, against 200^D for a grid — And it is exact, not approximate: the closed-form Beta(28,19) matched grid-based numerical Bayes to a relative 9.8e-15. Conjugacy is an algebraic identity.

2. Two datasets share the same n, sum and sum of squares but look completely different. What can a Gaussian likelihood say about them?

   A It will prefer whichever is closer to Gaussian B Nothing — the log-likelihoods are identical to 2.3e-13 at every parameter value, because those three numbers are the sufficient statistics. Theorem 6.14 says all the parameter information is in them C It will detect the difference at large sample sizes D The difference shows up in the variance estimate
   pch.quizShowAnswer

   B — Nothing — the log-likelihoods are identical to 2.3e-13 at every parameter value, because those three numbers are the sufficient statistics. Theorem 6.14 says all the parameter information is in them — Everything downstream of the likelihood inherits that blindness. If skew or multimodality matters to you, a Gaussian model cannot see it at any sample size — the information is not in phi(x).

3. In Example 6.14 the Bernoulli's natural parameter is theta = log(mu/(1-mu)). Why does that matter?

   A It makes the density easier to normalise B Because its inverse is the SIGMOID: mu is confined to (0,1) while theta ranges over all of R, so the sigmoid is what turns an unconstrained linear predictor into a probability. Logistic regression is a Bernoulli likelihood in natural parameters C Because it removes the need for a prior D Because theta equals the mean
   pch.quizShowAnswer

   B — Because its inverse is the SIGMOID: mu is confined to (0,1) while theta ranges over all of R, so the sigmoid is what turns an unconstrained linear predictor into a probability. Logistic regression is a Bernoulli likelihood in natural parameters — The same function reappears as a neural-network activation. Optimising in theta is unconstrained; optimising in mu would need a box constraint — which is precisely why models parameterise the predictor rather than the probability.

4. What is the log-partition function A(theta) for?

   A It is a constant that can be ignored B It normalises the density AND generates the moments — dA/dtheta = E[phi(x)], measured against the sigmoid to about 4e-11 for the Bernoulli. It depends on theta, so it does not vanish when you differentiate C It is the prior D It measures the sufficient statistic
   pch.quizShowAnswer

   B — It normalises the density AND generates the moments — dA/dtheta = E[phi(x)], measured against the sigmoid to about 4e-11 for the Bernoulli. It depends on theta, so it does not vanish when you differentiate — Dropping it is how gradients of exponential-family likelihoods go silently wrong. Keeping it turns maximum likelihood into 'match the empirical sufficient statistics to their expectations'.

5. Why is the Beta the conjugate prior for a Bernoulli — luck, or something else?

   A Historical convention B Equation 6.120 produces it. Substituting the Bernoulli's natural parameter and log-partition function collapses the inner product to mu^alpha (1-mu)^beta — verified to 1e-17 — and Example 6.15 reads off Equation 6.124, which is the Beta C Because both are defined on [0,1] D Because the Beta has two parameters
   pch.quizShowAnswer

   B — Equation 6.120 produces it. Substituting the Bernoulli's natural parameter and log-partition function collapses the inner product to mu^alpha (1-mu)^beta — verified to 1e-17 — and Example 6.15 reads off Equation 6.124, which is the Beta — Every exponential family member has a conjugate prior given by a formula, so Table 6.2 is derivable rather than a list to memorise. Example 6.12 assumed the Beta and checked it; Example 6.15 derives it.

🧪 Try It Yourself

Exercise 1 – Beta–Binomial conjugacy

pch.exerciseTag pch.exerciseTitle

Exercise 2 – Check it against brute force

pch.exerciseTag pch.exerciseTitle

Exercise 3 – Sufficiency

pch.exerciseTag pch.exerciseTitle

Exercise 4 – The Bernoulli in exponential-family form

pch.exerciseTag pch.exerciseTitle

Exercise 5 – The log-partition function gives the mean

pch.exerciseTag pch.exerciseTitle

Recall card

• Three desiderata: a closure property under the rules of probability, a parameter count that does not grow with the data, and well-behaved estimation. The exponential family is the class that satisfies all three.
• Pitman, Darmois and Koopman (1935-36) proved the exponential families are the ONLY families with finite-dimensional sufficient statistics under repeated independent sampling. The convenience is a theorem, not a coincidence.
• Def 6.13: a prior is conjugate for a LIKELIHOOD if the posterior stays in the prior’s family. Conjugacy is a property of the pair, which is why Table 6.2 lists pairs.
• Eq 6.104d: Beta-Binomial gives Beta(h + alpha, N - h + beta). The exponents simply add; no integral is computed and the normaliser never needs working out.
• alpha and beta are PSEUDO-COUNTS — alpha counts ones, beta counts zeros — which is why Beta(1,1) (“one of each”) is the uniform prior, and why a convenient prior is still an informative one.
• Conjugacy is exact. Closed-form Beta(28,19) matched grid-based numerical Bayes to a relative 9.8e-15, with both posterior means at 0.5957446809.
• The parameter count really does stay fixed. Measured: two numbers after 100000 observations, Beta(61789, 38213). A grid at the same resolution costs 200^D — 1.0e+23 at D = 10.
• The Beta has four regimes: uniform at alpha = beta = 1; bimodal with endpoint spikes when both are below 1; unimodal when both exceed 1; symmetric with mode 1/2 when they are equal and above 1.
• Theorem 6.14 (Fisher-Neyman): phi(x) is sufficient iff p(x|theta) = h(x) g_theta(phi(x)). All the parameter dependence flows through phi.
• Sufficiency is total. Two datasets sharing n, sum x and sum x-squared gave Gaussian log-likelihoods agreeing to 2.3e-13 while differing by 0.26 per order statistic. Anything downstream of the likelihood is equally blind.
• Eq 6.107: p(x|theta) = h(x) exp(theta . phi(x) - A(theta)). The data reaches the parameters ONLY through the inner product, so phi is sufficient by construction.
• theta are the NATURAL parameters and are unconstrained, unlike the usual ones. For a Bernoulli, theta = log(mu/(1-mu)) ranges over all of R while mu is confined to (0,1).
• The inverse of that link is the SIGMOID (Eq 6.118), which is why logistic regression is a Bernoulli likelihood in natural parameters, and why the same function turns up as an activation.
• A(theta) is not a constant to drop: dA/dtheta = E[phi(x)]. Measured against the sigmoid to about 4e-11. Its second derivative gives the covariance.
• Eq 6.120 generates a conjugate prior for EVERY exponential family member. Example 6.15 uses it to derive the Beta from the Bernoulli rather than guessing it — the exponential factor collapses to mu^alpha (1-mu)^beta, verified to 1e-17.

Next: Change of Variables and the Inverse Transform — what happens to a density when you transform the random variable, and the Jacobian that Chapter 5 built for exactly this moment.