Empirical Bayes

Take $N = 1000$ configurations whose true scores follow $g = 0.9\,\mathcal{N}(0, 0.3^2) + 0.1\,\mathcal{N}(2.5, 0.5^2)$ — 90% baseline configurations clustered near 0 and 10% genuinely better ones near 2.5 — with validation noise $\sigma = 1$. Rank by observed score and keep the top 20. On average, their observed scores overstate their true scores by +1.41.

Two corrections, previewed here and justified in §33.4–§33.5:

• The oracle Bayes rule $E[\theta \mid x]$, which knows $g$, overstates the top 20 by only +0.08.
• Tweedie’s formula with an estimated marginal — a Gaussian kernel density estimate (bandwidth $h = 0.35$) fit to the 1,000 observed scores, with $g$ unknown — overstates them by +0.21.

The leaderboard also hides an impostor: the configuration at rank 13 has true score 0.47. It is a baseline configuration that got lucky on the validation draw.

Topic 28 learned the two hyperparameters $(\mu, \tau^2)$ of a Normal prior $\theta_k \sim \mathcal{N}(\mu, \tau^2)$ by maximizing the marginal likelihood — Type-II maximum likelihood, §28.7. James–Stein (§28.5) was the striking special case: shrinking $d \geq 3$ Normal means toward a common point beats the MLE for every configuration of the means. Nothing here re-proves any of that. What changes is that $g$ is no longer forced to be Normal — it can be spiky, skewed, or multimodal, which is what real ensembles of effects look like. One notational change: Topic 28’s $K$ groups become this topic’s $N$ parallel problems, because $N$ will grow large and asymptotics in $N$ are the point.

Robbins (ROB1956) coined the phrase. A single problem’s Bayes rule needs a prior $g$, and with one observation there is no way to check it. But a large ensemble of parallel problems, each drawing its parameter from the same $g$, carries its own information about $g$: the histogram of the $x_i$ is a blurred picture of $g$. The prior stops being a statement of belief and becomes an estimable population distribution — the distribution of true effects across genes, of true quality across configurations, of true accident rates across drivers.

Fix real numbers $\theta = (\theta_1, \dots, \theta_N)$ — no prior, just $N$ unknown constants — and observe independent $X_i \sim p(\cdot \mid \theta_i)$, $i = 1, \dots, N$. A separable rule applies one function $t$ to every coordinate: $\hat\theta_i = t(X_i)$. Its compound risk is the average mean squared error

where $E_{\theta_i}$ averages over $X_i \sim p(\cdot \mid \theta_i)$ with $\theta_i$ held fixed.

For a prior $G$ on the parameter, the Bayes rule is the posterior mean $t_G(x) = E_G[\vartheta \mid X = x]$, computed with $\vartheta \sim G$ and $X \mid \vartheta \sim p(\cdot \mid \vartheta)$. With $G_N = \frac1N\sum_i \delta_{\theta_i}$ — the empirical distribution that puts mass $1/N$ on each true parameter — $t_{G_N}$ is the oracle Bayes rule (EFR2019’s “oracle Bayes”): it would be available only to someone who knew the histogram of the $\theta_i$. An empirical-Bayes rule $\hat t$ is any rule built from the data $X_1, \dots, X_N$ to estimate $t_{G_N}$, and its regret is $R_N(\hat t; \theta) - R_N(t_{G_N}; \theta)$.

Let $\theta_1, \dots, \theta_N$ be fixed real numbers and let $X_1, \dots, X_N$ be independent with $X_i \sim p(\cdot \mid \theta_i)$. For a separable rule $\hat\theta_i = t(X_i)$ with finite compound risk, and with $G_N = \frac1N\sum_{i=1}^N \delta_{\theta_i}$ the empirical distribution of the parameters:

(i) the compound risk equals the Bayes risk of $t$ under the prior $G_N$,

(ii) $R_N(t; \theta) \ge R_N(t_{G_N}; \theta)$ for every separable $t$, where

with equality if and only if $t = t_{G_N}$ almost everywhere under the marginal $m_{G_N}(x) = \frac1N\sum_i p(x \mid \theta_i)$.

f-modeling estimates the marginal density $f$ of the observations (or a functional of it, such as its log-derivative) and plugs the estimate into a formula for the Bayes rule $t_G$ that depends on $G$ only through $f$. g-modeling estimates the prior $G$ itself and computes the Bayes rule $t_{\hat G}$ of the estimate.

Restrict attention to linear rules $t(x) = a + bx$ with Normal noise, $X_i \sim \mathcal{N}(\theta_i, \sigma^2)$. For each coordinate, $t(X_i) - \theta_i = b(X_i - \theta_i) + (a + (b - 1)\theta_i)$ is a mean-zero noise term plus a constant, so its mean square is $b^2\sigma^2 + (a + (b-1)\theta_i)^2$. Averaging over $i$,

For fixed $b$, the second term is the average squared distance between $a$ and the points $(1-b)\theta_i$, minimized at their mean: $a = (1-b)\bar\theta$, where $\bar\theta = \frac1N\sum_i\theta_i$. Substituting, with $s^2 = \frac1N\sum_i(\theta_i - \bar\theta)^2$ the spread of the true parameters,

Setting the derivative in $b$ to zero, $2b\sigma^2 - 2(1-b)s^2 = 0$, gives $b^* = s^2/(s^2 + \sigma^2)$. The best linear rule is therefore $\bar\theta + b^*(x - \bar\theta) = x - (1 - b^*)(x - \bar\theta)$, with shrinkage fraction $1 - b^* = \sigma^2/(s^2 + \sigma^2)$.

Now estimate the two unknowns. Replace $\bar\theta$ by $\bar X$, and estimate $1 - b^*$ by $(N-3)\sigma^2/\sum_i (X_i - \bar X)^2$. Under Normal-prior sampling, $\theta_i \sim \mathcal{N}(\mu, \tau^2)$, the observations are marginally $\mathcal{N}(\mu, \sigma^2 + \tau^2)$, so $\sum_i (X_i - \bar X)^2 \sim (\sigma^2 + \tau^2)\chi^2_{N-1}$; since $E[1/\chi^2_{N-1}] = 1/(N-3)$, the estimator has expectation exactly $\sigma^2/(\sigma^2 + \tau^2)$. The result is the Efron–Morris rule (EFR1973),

James–Stein is not a paradox: it is the empirical-Bayes estimate of the best linear separable rule. Nonparametric empirical Bayes drops the word “linear.”

No prior was assumed anywhere in Theorem 1 — the $\theta_i$ are fixed — and yet the optimal separable rule is a posterior mean under $G_N$. Robbins called procedures that approach this optimum without knowing $G_N$ “asymptotically subminimax” (ROB1951): for large $N$ they do nearly as well as if the histogram of the true parameters had been revealed. Separable rules are the natural class because the problems are exchangeable before we see the data — nothing distinguishes problem 7 from problem 700. A non-separable rule can beat $t_{G_N}$ only by exploiting the ordering of the indices, which carries no information when the labels are arbitrary.

Figure 2 shows why the two strategies are genuinely different. A prior $g$ can be spiky — here, three point masses — but the data never see it directly. They see the marginal $f = g \ast \varphi_\sigma$: the prior convolved with (smeared by) the Normal noise, which is always smooth. f-modeling works with that smooth object; g-modeling tries to undo the smearing. f-modeling is Robbins’ formula (§33.3) and Tweedie’s formula (§33.4–§33.5); g-modeling is the nonparametric MLE of the prior (§33.6); §33.7 compares them head to head.

Let $\theta \sim g$, any probability measure on $(0, \infty)$, and $X \mid \theta \sim \mathrm{Poisson}(\theta)$, with marginal

Then for every $x$ with $f(x) > 0$,

CASI Table 6.1 (EFR2016) records one year of claims for 9,461 policy holders of a European auto insurer:

These are identical to CASI’s published row. The first entry is the cold-start answer: customers with no claims this year should expect $1317/7840 = 0.168$ claims next year — not 0, which is what the MLE $\hat\theta = x$ says. Customers with one claim should expect 0.363, far less than 1: one claim is weak evidence of a high rate, because most of the book has low rates and some of them file a claim by chance.

The parametric alternative assumes a shape for $g$: $\theta \sim \mathrm{Gamma}(\nu, s)$ with shape $\nu$ and scale $s$. Then the posterior given $x$ is $\mathrm{Gamma}(\nu + x,\, s/(1+s))$, with mean

and the marginal of $X$ is negative binomial. Its maximum-likelihood fit preserves the sample mean, $\hat\nu\hat s = \bar x = 0.2144$, so the fit reduces to a one-dimensional profile in $\nu$, which gives $\hat\nu = 0.7015$ and $\hat s = 0.3056$. The fitted rule for $x = 0, \dots, 6$ is

0.164, 0.398, 0.632, 0.866, 1.100, 1.334, 1.568

(CASI: .164 .398 .633 .87 1.10 1.34 1.57). Figure 3 puts the two side by side: they agree where the counts are large, and Robbins goes erratic where they are small.

At $x = 5$, Robbins’ estimate is $6 \cdot y_6/y_5 = 6 \cdot 4/4 = 6.00$: two tiny counts drive it. A single extra policy holder with six claims would change it to $6 \cdot 5/4 = 7.5$. The parametric fit borrows strength across every $x$; Robbins borrows only from the neighboring cell $x + 1$. The explorer’s perturbation buttons make the difference tactile — perturbing $y_5$ moves Robbins at $x = 4$ and $x = 5$ substantially and leaves the Gamma rule essentially unchanged.

For the Poisson family, $E[\theta \mid x]$ is nondecreasing in $x$ for every prior. The posterior at $x + 1$ is the posterior at $x$ tilted by a factor of $\theta$, so

because $E[\theta^2 \mid x] - (E[\theta \mid x])^2 = \mathrm{Var}(\theta \mid x) \ge 0$. This is the discrete twin of §33.4 Theorem 3(iii). Robbins’ drop from 6.00 at $x = 5$ to 1.75 at $x = 6$ is therefore an estimation artifact, not a feature of any prior. We’ll turn this observation into a diagnostic in §33.7.

Let $g$ be any probability measure on $\mathbb{R}$, $\theta \sim g$, and $X \mid \theta \sim \mathcal{N}(\theta, \sigma^2)$ with $\sigma^2 > 0$ known. Write $\varphi_\sigma(u) = (2\pi\sigma^2)^{-1/2}e^{-u^2/(2\sigma^2)}$ and

Then $f$ is strictly positive and twice differentiable, and for every $x$:

(i) $E[\theta \mid X = x] = x + \sigma^2\,\ell'(x)$;

(ii) $\mathrm{Var}(\theta \mid X = x) = \sigma^2\big(1 + \sigma^2\,\ell''(x)\big)$;

(iii) $\dfrac{d}{dx}E[\theta \mid X = x] = \mathrm{Var}(\theta \mid X = x)/\sigma^2 \ge 0$ — the Bayes rule is nondecreasing in $x$.

Take $g = \mathcal{N}(\mu, \tau^2)$. The marginal is $f = \mathcal{N}(\mu, \sigma^2 + \tau^2)$ — a Normal prior smeared by Normal noise is Normal with the variances added — so

Tweedie’s formula returns $x - \sigma^2(x - \mu)/(\sigma^2 + \tau^2) = (1 - B)x + B\mu$ with $B = \sigma^2/(\sigma^2 + \tau^2)$, and part (ii) gives the posterior variance $\sigma^2(1 - B) = (1 - B)\sigma^2$. This is exactly Topic 28’s partial-pooling formula (§28.6, Thm 4), now a one-line special case. The rule is linear because $\ell$ is quadratic; every non-Normal prior bends it.

Take $g = 0.9\,\delta_0 + 0.1\,\mathcal{N}(0, 3^2)$ — 90% of the effects are exactly zero, the rest are spread widely — with $\sigma = 1$ and $N = 2000$ draws. The oracle rule $E[\theta \mid x]$ at $x = 1, 2, 3, 5$ is 0.047, 0.316, 1.805, 4.498: near-total shrinkage near 0, where an observation is most likely a null plus noise, and near-identity in the tails, where it is most likely a real effect.

For the f-modeling estimate, fit a Gaussian KDE (Topic 30) with Silverman’s bandwidth $h = 0.32$: $\hat f_h(x) = \frac{1}{Nh}\sum_j \varphi((x - X_j)/h)$. Differentiating term by term and dividing by $\hat f_h$, the KDE score has the closed form

a softmax-weighted pull toward nearby data points. Per-coordinate mean squared error on the sample: the MLE $\hat\theta_i = x_i$ scores 1.033, Tweedie–KDE 0.244, the oracle 0.193. With $g$ unknown, Tweedie–KDE captures about 94% of the oracle’s improvement over the MLE.

Part (i) says the correction $\sigma^2\ell'(x)$ points toward higher marginal density. Observations in sparse regions get pulled toward where the data are dense — Figure 4’s arrows all point at the spike. The noise variance $\sigma^2$ sets the step size: noisier measurements are pulled harder, because they say less about their own $\theta$. That is the whole content of shrinkage, stated without a prior.

The rule depends on $g$ only through $f$, which is what makes f-modeling possible: estimate $f$ from the $x_i$, differentiate its logarithm, done. The price is that we need a derivative of a density estimate, and derivative estimates are noisy where data are thin — Example 6’s KDE rule wiggles in the tails. Bandwidth controls the trade-off exactly as in Topic 30: see the MISE/AMISE analysis in §30.5 and data-driven bandwidth selection in §30.8.

Corrupt a clean sample $x_0$ with Gaussian noise, $x_t = x_0 + \sigma_t\varepsilon$ with $\varepsilon \sim \mathcal{N}(0, I)$. Theorem 3 in $\mathbb{R}^d$ — the proof runs coordinate by coordinate — gives

where $p_t$ is the density of the noisy samples. A denoiser trained by squared error, minimizing $E\|D(x_t) - x_0\|^2$, converges to the posterior mean $E[x_0 \mid x_t]$, so $(D(x_t) - x_t)/\sigma_t^2$ estimates the score $\nabla\log p_t$. That is denoising score matching (VIN2011), the training objective behind score-based diffusion models: every denoising network is an f-modeling estimate of Tweedie’s correction. formalML’s representation-learning topic develops the identity (§3.4 there, Thm 3.2): formalML: representation-learning .

Let $p(x \mid \eta) = e^{\eta x - \psi(\eta)}h(x)$ be a one-parameter exponential family of densities on an open interval $\mathcal{X}$, in Topic 7’s notation (natural parameter $\eta$, cumulant function $\psi$, carrier $h > 0$ twice differentiable on $\mathcal{X}$). Let $\eta \sim g$, a probability measure on the natural parameter space with $\int (1 + |\eta| + \eta^2)\, e^{\eta x - \psi(\eta)}\, dg(\eta) < \infty$ for every $x \in \mathcal{X}$. With $f(x) = \int p(x \mid \eta)\,dg(\eta)$ and $\ell = \log f$:

(i) $E[\eta \mid X = x] = \ell'(x) - (\log h)'(x)$; (ii) $\mathrm{Var}(\eta \mid X = x) = \ell''(x) - (\log h)''(x)$.

For an integer-valued family $p(x \mid \eta) = e^{\eta x - \psi(\eta)}h(x)$, derivatives in $x$ are replaced by shifts. Dividing the marginal by the carrier,

and the ratio of the two is a posterior expectation, because the posterior given $x$ is proportional to $e^{\eta x - \psi(\eta)}\,dg(\eta)$:

For the Poisson family, $p(x \mid \theta) = e^{-\theta}\theta^x/x! = \exp(x\log\theta - \theta)/x!$, so $\eta = \log\theta$, $e^\eta = \theta$, and $h(x) = 1/x!$, which makes $h(x)/h(x+1) = x + 1$. The display becomes $E[\theta \mid x] = (x+1)f(x+1)/f(x)$ — Theorem 2.

Bin the data into $K$ equal-width bins with centers $c_k$, width $\Delta$, and counts $y_k$. Fit the Poisson regression

where $u_k$ is the standardized bin center, and set $\hat f(x) = \exp\big(\sum_j \beta_j u(x)^j\big)/(N\Delta)$. The estimated score $\hat\ell'$ is the derivative of a polynomial — smooth everywhere. The fit is exactly the Poisson generalized linear model of Topic 22, with histogram counts as the response.

Let $(\theta_i, X_i)$, $i = 1, \dots, N$, be independent pairs with $\theta_i \sim g$ and $X_i \mid \theta_i \sim p(\cdot \mid \theta_i)$, and let $t_g(x) = E[\theta \mid X = x]$. Let $S \subseteq \{1, \dots, N\}$ be any selection rule that is a function of $\mathbf{X} = (X_1, \dots, X_N)$ — the indices of the $k$ largest $X_i$, or all $i$ with $X_i > c$. Then for every $i$ with $P(i \in S) > 0$,

Fit Lindsey’s method to the 6,033 prostate $z$-scores with 90 bins on $[\min z, \max z]$ and degree $J = 6$. The leading coefficient is negative, so the fitted density is integrable, and $1 + \hat\ell''(z) \ge 0.030$ across the data range, so the estimated rule $z + \hat\ell'(z)$ is monotone, as Theorem 3(iii) says a Bayes rule must be. For the ten largest $z$-scores: