Sufficiency and the factorization theorem · Probability & Statistics

You flip a coin 1,000 times and want to estimate p. Do you need the full sequence of heads and tails, or is the count enough?

The count is enough — and "enough" has a precise meaning.

The definition

A statistic T(X) is sufficient for \theta if the conditional distribution of the data given T does not depend on \theta:

$P\big(X_1, \dots, X_n \mid T(X) = t, \;\theta\big) = P\big(X_1, \dots, X_n \mid T(X) = t\big)$

Read it as: once you know T, the rest of the data tells you nothing more about \theta. Every scrap of information the sample carries about the parameter has been compressed into T.

You could throw the raw data away, keep only T, and lose nothing — for the purpose of learning \theta.

The coin example, concretely

Observe 1,000 flips with 300 heads. Given that the count is 300, every specific arrangement of those 300 heads is equally likely — there are \binom{1000}{300} of them, each with conditional probability 1/\binom{1000}{300}.

$P(\text{specific sequence} \mid k = 300) = \frac{p^{300}(1-p)^{700}}{\binom{1000}{300}p^{300}(1-p)^{700}} = \frac{1}{\binom{1000}{300}}$

The p cancels completely. The particular pattern is uninformative about p; only the count matters. So T = \sum X_i is sufficient.

The factorization theorem

Verifying sufficiency from the definition means computing conditional distributions, which is painful. The Fisher–Neyman factorization theorem makes it mechanical:

$T(X) \text{ is sufficient for } \theta \iff f(x \mid \theta) = g\big(T(x), \theta\big)\, h(x)$

The likelihood must split into:

g, which involves \theta and touches the data only through T(x)
h, which may involve the data however it likes but must not contain \theta

That's it. Factor the likelihood and read off T.

Why it works: the h(x) part is a constant multiplier as far as \theta is concerned. It scales the likelihood but never changes its shape, so it cannot affect which \theta looks best. All the parameter-relevant structure lives in g(T(x), \theta).

Worked example 1: Bernoulli

$f(x \mid p) = \prod p^{x_i}(1-p)^{1-x_i} = p^{\sum x_i}(1-p)^{n - \sum x_i}$

Take T(x) = \sum x_i, then

$g(T, p) = p^{T}(1-p)^{n-T}, \qquad h(x) = 1$

T = \sum X_i is sufficient. The full sequence is irrelevant once you have the total.

Worked example 2: Normal with known \sigma^2

$f(x \mid \mu) \propto \exp\!\left(-\frac{1}{2\sigma^2}\sum(x_i - \mu)^2\right)$

Expand the square:

$\sum(x_i - \mu)^2 = \sum x_i^2 - 2\mu\sum x_i + n\mu^2$

$f(x \mid \mu) = \underbrace{\exp\!\left(\frac{\mu\sum x_i}{\sigma^2} - \frac{n\mu^2}{2\sigma^2}\right)}_{g(\sum x_i,\; \mu)} \times \underbrace{(2\pi\sigma^2)^{-n/2}\exp\!\left(-\frac{\sum x_i^2}{2\sigma^2}\right)}_{h(x)}$

\mu appears only alongside \sum x_i, and the \sum x_i^2 term has no \mu in it. T = \sum X_i (equivalently \bar X) is sufficient for \mu.

Any one-to-one function of a sufficient statistic is also sufficient — $\sum X_i$ and \bar X carry identical information.

Worked example 3: Normal with both unknown

Now \theta = (\mu, \sigma^2) and the \sum x_i^2 term can no longer be dumped into h, since h must be free of parameters:

$f(x \mid \mu, \sigma^2) = (2\pi\sigma^2)^{-n/2}\exp\!\left(-\frac{\sum x_i^2 - 2\mu\sum x_i + n\mu^2}{2\sigma^2}\right)$

Everything involving the data appears through the pair (\sum x_i, \sum x_i^2), so

$T = \left(\sum X_i, \; \sum X_i^2\right)$

is jointly sufficient — equivalently (\bar X, s^2).

Two parameters, two sufficient statistics. That's typical, and it's why \bar x and s^2 are the summary you always see reported for a Normal sample: together they are the entire information content of the data.

Worked example 4: Uniform

$f(x \mid \theta) = \frac{1}{\theta^n}\mathbb{1}\{x_{(n)} \le \theta\}\cdot\mathbb{1}\{x_{(1)} \ge 0\}$

Take g(T, \theta) = \theta^{-n}\mathbb{1}\{T \le \theta\} with T = x_{(n)}, and h(x) = \mathbb{1}\{x_{(1)} \ge 0\}.

The maximum is sufficient. This is the formal statement of what we kept noticing: X_{(n)} carries all the information about \theta, which is why it beat 2\bar X on MSE and why maximum likelihood found it.

The indicator function is doing the work here — that's typical for distributions whose support depends on the parameter.

Why sufficiency matters

Rao–Blackwell theorem. Given any unbiased estimator \hat\theta and a sufficient statistic T, define \tilde\theta = E[\hat\theta \mid T]. Then \tilde\theta is unbiased and

$\operatorname{Var}(\tilde\theta) \le \operatorname{Var}(\hat\theta)$

Conditioning on a sufficient statistic never makes an estimator worse. So good estimators should always be functions of sufficient statistics — anything else is discarding usable information, or keeping noise.

This is the theoretical justification for the whole approach: it tells you where to look for optimal estimators.

Minimal sufficiency. A sufficient statistic that is a function of every other sufficient statistic — the most compressed summary that loses nothing. The whole dataset is always trivially sufficient; minimal sufficiency is the useful end of the scale.

Exponential families. Distributions writable as

$f(x \mid \theta) = h(x)\exp\big(\eta(\theta)T(x) - A(\theta)\big)$

have T(x) sufficient by inspection. Normal, Bernoulli, Poisson, Gamma, Beta and more all belong, which is why they behave so well and why the same few summaries keep appearing.

Doing it in Python

Sufficiency is a claim about likelihoods, so check it on likelihoods: two datasets with the same sufficient statistic must give the same shape of likelihood.

import numpy as np

# Two very different Bernoulli sequences with the same number of successes
a = np.array([1, 1, 1, 0, 0, 0, 0, 0, 0, 0])       # all successes first
b = np.array([0, 1, 0, 1, 0, 0, 1, 0, 0, 0])       # scattered
print("same sum?", a.sum() == b.sum(), "  (T =", a.sum(), ")")

def loglik(x, p):
    return x.sum() * np.log(p) + (len(x) - x.sum()) * np.log(1 - p)

print(f"\n{'p':>6} {'loglik(a)':>12} {'loglik(b)':>12} {'difference':>12}")
for p in (0.1, 0.3, 0.5, 0.7):
    la, lb = loglik(a, p), loglik(b, p)
    print(f"{p:>6} {la:>12.5f} {lb:>12.5f} {la - lb:>12.2e}")

print("\nIdentical for every p -- the arrangement carries no information about p.")

A contrast: a statistic that is not sufficient loses information.

import numpy as np

rng = np.random.default_rng(0)
theta, n, trials = 10.0, 20, 100_000
x = rng.uniform(0, theta, size=(trials, n))

# Sufficient (the max) vs not sufficient (the mean)
from_max = (n + 1) / n * x.max(axis=1)
from_mean = 2 * x.mean(axis=1)

for name, est in [("from max (sufficient)", from_max), ("from mean (not)", from_mean)]:
    bias = est.mean() - theta
    print(f"{name:>24}  MSE {est.var() + bias**2:>8.4f}")
print("\nThe non-sufficient summary throws away information, and pays for it.")

Rao–Blackwell in action — conditioning a crude estimator on a sufficient statistic:

import numpy as np

rng = np.random.default_rng(1)
p, n, trials = 0.3, 20, 200_000
x = rng.random((trials, n)) < p

crude = x[:, 0].astype(float)        # "just use the first flip" -- unbiased, awful
T = x.sum(axis=1)
rao_black = T / n                    # E[X_1 | T] = T/n

for name, est in [("first flip", crude), ("E[first flip | T]", rao_black)]:
    print(f"{name:>18}  mean {est.mean():.4f}  variance {est.var():.6f}")
print(f"\nBoth unbiased for p = {p}; conditioning cut the variance by "
      f"{crude.var() / rao_black.var():.0f}x.")

Your turn

1. X_i \sim \text{Poisson}(\lambda). Find a sufficient statistic.

2. X_i \sim \text{Exponential}(\lambda). Find a sufficient statistic.

3. Is \bar X sufficient for \theta when $X_i \sim \text{Uniform}(0,\theta)$?

Solutions

$f(x \mid \lambda) = \prod\frac{\lambda^{x_i}e^{-\lambda}}{x_i!} = \underbrace{\lambda^{\sum x_i}e^{-n\lambda}}_{g(\sum x_i,\;\lambda)}\cdot\underbrace{\frac{1}{\prod x_i!}}_{h(x)}$

\lambda touches the data only through \sum x_i, and the factorial term is parameter-free. So T = \sum X_i is sufficient.

$f(x \mid \lambda) = \prod\lambda e^{-\lambda x_i} = \underbrace{\lambda^n e^{-\lambda\sum x_i}}_{g(\sum x_i,\;\lambda)}\cdot\underbrace{1}_{h(x)}$

Again T = \sum X_i is sufficient. (A pattern is emerging: for exponential families, the sufficient statistic is whatever the parameter multiplies in the exponent.)

3. No.

The likelihood is \theta^{-n}\mathbb{1}\{x_{(n)} \le \theta\}, and there's no way to write that with the data entering only through \bar x — the indicator depends on the maximum, which \bar x cannot recover.

The concrete failure: the samples \{1, 5\} and \{3, 3\} have the same mean of 3, but the first proves \theta \ge 5 while the second only proves \theta \ge 3. Same \bar x, genuinely different information about \theta — so \bar X cannot be sufficient.

This is the precise reason 2\bar X underperformed X_{(n)} back in lesson 1. Sufficiency explains why the MSE comparison came out as it did.

Check yourself in code

Demonstrate sufficiency for the Bernoulli: two different sequences with the same number of successes must produce identical log-likelihoods at every p.

Print exactly this:

same sum: True
max abs difference: 0.0
sufficient: True

Compare the two sequences at p \in \{0.1, 0.3, 0.5, 0.7, 0.9\} and report the largest absolute difference, rounded to 10 decimal places.

import numpy as np

a = np.array([1, 1, 1, 0, 0, 0, 0, 0, 0, 0])
b = np.array([0, 1, 0, 1, 0, 0, 1, 0, 0, 0])

print("same sum:", bool(a.sum() == b.sum()))

def loglik(x, p):
    return x.sum() * np.log(p) + (len(x) - x.sum()) * np.log(1 - p)

# Compare the log-likelihoods across p in {0.1, 0.3, 0.5, 0.7, 0.9},
# report the largest absolute difference, and whether it is zero.

import numpy as np

a = np.array([1, 1, 1, 0, 0, 0, 0, 0, 0, 0])
b = np.array([0, 1, 0, 1, 0, 0, 1, 0, 0, 0])

print("same sum:", bool(a.sum() == b.sum()))

def loglik(x, p):
    return x.sum() * np.log(p) + (len(x) - x.sum()) * np.log(1 - p)

diffs = [abs(loglik(a, p) - loglik(b, p)) for p in (0.1, 0.3, 0.5, 0.7, 0.9)]
worst = round(max(diffs), 10)
print("max abs difference:", worst)
print("sufficient:", worst == 0.0)

A sufficient statistic compresses the data without losing any information about the parameter. The factorization theorem turns checking that into a mechanical exercise: split the likelihood into a parameter part that sees the data only through T, and a data part free of the parameter. Rao–Blackwell then says good estimators should always be functions of sufficient statistics.

Next: how to quantify exactly how much information a sample carries — and the hard limit that places on any estimator.

29. Sufficiency and the factorization theorem