Stein's Paradox & Empirical Bayes

1. The sample mean in high dimensions

Suppose you see $X \sim \mathcal{N}(\theta, 1)$ in one dimension and you want a guess for $\theta$. The natural guess is just $X$ itself and no other estimator beats it across the board.

In higher dimensions the picture shifts hard. Suppose you see $X \sim \mathcal{N}(\theta, I)$ in dimension $d = 100$ and you try the same guess $\hat{\theta} = X$. It has a built-in flaw because it overshoots $\|\theta\|$ every time.

By the law of large numbers, as $d \to \infty$:

Meanwhile, the squared norm of $X$ satisfies:

So $X$ sits on a shell of radius about $\sqrt{\|\theta\|^2 + d}$ and that is way farther from the origin than $\theta$ is. The noise piece $\|Z\|^2$ piles up around $d$ and the cross term $\theta \cdot Z$ is only $O(\sqrt{d})$ and the noise piece lands in directions orthogonal to $\theta$ because of how concentration of measure works in high dimensions.

The obvious fix is to pull $X$ back toward the origin and Stein proved in 1956 that this fix is not just nice to have but actually required. $d \ge 3 \implies \delta_0(X) = X \text{ is inadmissible.}$

2. The James-Stein estimator

Set up risk under squared error loss.

For $\delta_0(X) = X$ the risk is $d$ no matter what $\theta$ is.

The James-Stein estimator uses a shrinkage factor that depends on the data.

Theorem (Stein’s Paradox). For all $d \ge 3$ and all $\theta \in \mathbb{R}^d$,

The improvement is largest near the origin, and at $\theta = 0$ we have

In dimension 100 the MLE has risk 100 and James-Stein has risk 2 and that is a 98 percent cut.

To see why shrinkage is the natural move geometrically, think about the ideal or oracle shrinkage factor. Given $X = \theta + \epsilon$ the scalar $c$ that makes $\|\theta - cX\|^2$ smallest is

The oracle shrinkage factor matches the James-Stein form and $d$ gets swapped for $d - 2$ to handle the estimation error in the denominator.

3. Stein’s Lemma

To work out the risk of a nonlinear estimator $\delta(X) = X + g(X)$ without knowing $\theta$ we need the following identity.

Lemma (Stein, 1981). Let $X \sim \mathcal{N}_d(\theta, I)$ and let $g: \mathbb{R}^d \to \mathbb{R}^d$ be weakly differentiable with $\mathbb{E} | \nabla \cdot g(X) | < \infty$. Then

Proof. We go component by component. Because the coordinates are independent the joint density factors as $p(x) = \prod \phi(x_j - \theta_j)$ and each component can be handled on its own.

For the $i$-th component,

Pull out the integral over $x_i$,

The Gaussian kernel obeys $\phi'(z) = -z \phi(z)$ and so $(x_i - \theta_i)\phi(x_i - \theta_i) = -\frac{\partial}{\partial x_i} \phi(x_i - \theta_i)$. Plugging this in gives

Integrating by parts with $u = g_i(x)$ and $dv = \phi'(x_i - \theta_i) dx_i$,

The boundary term dies whenever $g_i$ grows slower than $e^{x^2/2}$ and that covers every polynomial-growth estimator and we are left with

Adding over $i = 1, \ldots, d$ gives

$\square$

4. Deriving the James-Stein estimator

Write $\delta(X) = X + g(X)$ and open up the risk,

The first term equals $d$ and by Stein’s Lemma the cross term equals $2 \mathbb{E}[\nabla \cdot g(X)]$ and so

The quantity $\|g\|^2 + 2 \nabla \cdot g$ is Stein’s Unbiased Risk Estimate and people call it SURE. To beat the MLE we need a nonzero $g$ that pushes this expression below zero.

Try $g(X) = - \frac{c}{\|X\|^2} X$. It blows up at $X = 0$ but that event has probability zero under a continuous distribution.

To work out the divergence note that since $g_i(X) = - c X_i (X_1^2 + \dots + X_d^2)^{-1}$ the partial derivative is

Adding over all the components gives

The squared norm of $g$ is

Plugging into SURE gives

The risk improvement is quadratic in $c$ and its minimum sits at $c^* = d-2$ and that gives

For $d \ge 3$ this is strictly negative and so

When $\theta = 0$ the distribution $\|X\|^2 \sim \chi^2_d$ gives $\mathbb{E}[1/\chi^2_d] = 1/(d-2)$ and so

The paradox is that a single estimator beats the MLE at the same time for every single value of the parameter.

5. The Bayesian derivation

SURE shows you the risk reduction is real and the Bayesian view explains why the estimator has this exact shape.

Put a Gaussian prior on $\theta$.

Level 1: $X | \theta \sim \mathcal{N}_d(\theta, I)$

Level 2: $\theta \sim \mathcal{N}_d(0, A \cdot I)$

By conjugacy the posterior is Gaussian and the precisions add,

The posterior mean is

Writing $B = \frac{1}{1+A}$ the Bayes estimate is $(1-B)X$ and that is linear shrinkage toward the origin.

The snag is that $A$ is unknown but the marginal distribution $X \sim \mathcal{N}(0, (1+A)I)$ gives us enough to pin down $A$ from the data,

Since $\mathbb{E}[1/\chi^2_d] = 1/(d-2)$ the quantity $\frac{d-2}{\|X\|^2}$ is an unbiased estimator of $B = \frac{1}{1+A}$ and

Plugging this empirical estimate of the shrinkage factor into the Bayes rule gives you the James-Stein estimator exactly.

The empirical Bayes story makes what the estimator is doing clear because it learns the signal-to-noise ratio straight from the data. When $\|X\|^2$ is big compared to $d$ the data says there is strong signal and the estimator barely shrinks and when $\|X\|^2$ is small the data is mostly noise and the estimator shrinks hard.

6. Admissibility and Brownian motion

There is a deep link between whether the MLE is admissible and how Brownian motion recurs and Brown worked this out in 1971.

Think about estimators of the form $\delta(X) = X + \nabla \phi(X)$. Through Tweedie’s formula any such estimator can be read as a formal Bayes estimator against $\pi(\theta) = e^{\phi(\theta)}$.

The risk improvement of $\delta$ over $X$ is controlled by the Laplacian of $\sqrt{m}$ where $m$ is the marginal density and

If $\sqrt{m}$ is superharmonic (meaning $\Delta \sqrt{m} \le 0$) then the risk improves. The flat prior $\pi = 1$ gives $m = 1$ and $\Delta 1 = 0$ and that is the boundary case that matches the MLE.

The harmonic prior $\pi(\theta) \propto \|\theta\|^{-(d-2)}$ gives a superharmonic marginal for $d \ge 3$.

The algebra reflects something topological. In $d = 1$ and $d = 2$ Brownian motion is recurrent and keeps coming back to every neighborhood over and over and in $d \ge 3$ Brownian motion is transient and runs off to infinity.

The inadmissibility of the MLE in $d \ge 3$ lines up with probability mass leaking off to infinity under the flat prior and the James-Stein estimator makes up for that leak. Shrinkage is needed exactly because the surrounding space is high-dimensional enough for random walks to escape.

7. The positive-part correction

The standard James-Stein estimator has a failure mode. When $\|X\|^2 < d - 2$ the shrinkage factor $1 - \frac{d-2}{\|X\|^2}$ goes negative and the estimate flips direction and that is obviously bad because the estimator shoves $\hat\theta$ over to the opposite side of the origin from $X$.

Baranchik’s positive-part tweak clips the shrinkage factor,

This estimator dominates $\delta_{JS}$ but it is not itself admissible because it has non-analytic kinks at $\|X\|^2 = d - 2$ and Brown showed that admissible estimators for exponential families have to be analytic.

Nobody knows a simple analytic estimator that beats JS+ and so the positive-part estimator is what people actually use in practice.

8. Exact risk of the positive-part estimator

Working out the exact risk of $\delta_{JS+}$ takes some care because the estimator acts differently in two regions.

The difference between $\delta_{JS}$ and $\delta_{JS+}$ is only nonzero when the shrinkage overshoots and that happens when $\|X\|^2 < d-2$. In that region $\delta_{JS}(X) = (1 - \frac{d-2}{\|X\|^2})X$ and $\delta_{JS+}(X) = 0$.

Let $A = \{ x : \|x\|^2 < d-2 \}$. The risk gap between the two estimators is

Opening up the risk gap on $A$ using Stein’s lemma on the full space and handling the truncation gives

The cross term $-2(1 - \frac{d-2}{\|X\|^2})(X \cdot \theta)$ does not die on $A$ because the truncation breaks the symmetry that Stein’s lemma leans on over the full space. Even so you can check that $\Delta R_+ < 0$ for every $\theta$ by writing

because on $A$ the factor $(1 - \frac{d-2}{\|x\|^2})$ is negative and $|(1 - \frac{d-2}{\|x\|^2})| \cdot \|x\| < \|x\|$ and so the JS estimate shoots past the origin and makes a bigger error than just guessing $\theta$ is zero.

The exact risk of $\delta_{JS+}$ needs you to integrate against the non-central chi-square distribution $Y = \|X\|^2 \sim \chi^2_d(\|\theta\|^2)$ and that has a Poisson mixture form

where $\lambda = \|\theta\|^2$.

The risk of JS is

The extra improvement from JS+ means you integrate the full risk gap including the $\theta$-dependent cross term over the region $\|X\|^2 < d-2$ and you have to do this numerically for general $\theta$. The qualitative takeaway that $R_{JS+} < R_{JS}$ for every $\theta$ falls out of the geometric argument above.

9. Simulation

The simulation below compares the risk of the MLE and James-Stein and the positive-part estimators as a function of $\|\theta\|$. The MLE risk is a flat line at $d$ and the shrinkage estimators get their biggest win near the origin.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import noncentral_chisquare
 
def simulate_risk(d=10, n_trials=5000, theta_norms=None):
    if theta_norms is None:
        theta_norms = np.linspace(0, 10, 20)
        
    risk_mle = []
    risk_js = []
    risk_plus = []
    
    for r in theta_norms:
        # Create theta vector with norm r
        theta = np.zeros(d)
        theta[0] = r
        
        # Generate Data
        # X ~ N(theta, I). Shape (n_trials, d)
        X = np.random.randn(n_trials, d) + theta
        
        # Norms squared
        X_norm_sq = np.sum(X**2, axis=1)
        
        # MLE Error
        # ||X - theta||^2
        loss_mle = np.sum( (X - theta)**2, axis=1 )
        risk_mle.append(np.mean(loss_mle))
        
        # JS Estimator
        # (1 - (d-2)/||X||^2) * X
        shrinkage = 1 - (d-2)/X_norm_sq
        Theta_JS = X * shrinkage[:, np.newaxis]
        loss_js = np.sum( (Theta_JS - theta)**2, axis=1 )
        risk_js.append(np.mean(loss_js))
        
        # JS+ Estimator
        shrinkage_plus = np.maximum(0, shrinkage)
        Theta_Plus = X * shrinkage_plus[:, np.newaxis]
        loss_plus = np.sum( (Theta_Plus - theta)**2, axis=1 )
        risk_plus.append(np.mean(loss_plus))
        
    # Plot
    plt.figure(figsize=(10, 6))
    plt.plot(theta_norms, risk_mle, 'k--', label='MLE Risk (d)')
    plt.plot(theta_norms, risk_js, 'b-o', label='James-Stein Risk')
    plt.plot(theta_norms, risk_plus, 'r-', label='Positive Part JS')
    
    plt.axhline(d, color='gray', linestyle=':')
    plt.axhline(2, color='green', linestyle=':', label='Min Risk (at 0)')
    
    plt.xlabel('||Theta||')
    plt.ylabel('Risk (MSE)')
    plt.title(f'Stein\'s Paradox in {d} Dimensions')
    plt.legend()
    # plt.show()
    return theta_norms, risk_js
 
# Theoretical Minimum Risk:
# At theta=0, Risk = d - (d-2)^2 E[1/ChiSq_d].
# E[1/ChiSq_d] = 1/(d-2).
# Risk(0) = d - (d-2) = 2.
# Massive dependency! From d to 2.

10. The Efron-Morris baseball example

Efron and Morris (1975) ran shrinkage estimation on batting averages. They took 18 players’ first 45 at-bats and used them to predict how those players would do over the rest of the season.

The MLE uses each player’s own 45-at-bat average and the James-Stein estimator pulls every player toward the grand mean.

James-Stein cut the total squared prediction error by a factor of 3. Players who started unusually hot at batting .450 got pulled down toward the league average and players who started cold got pulled up and in both cases the shrunk estimates landed closer to the true end-of-season averages.

Shrinkage toward a common value is the same trick that underlies ridge regression which is Bayesian estimation with a Gaussian prior centered at zero. James-Stein can be read as ridge regression where the regularization parameter tunes itself from the data as $\lambda = \frac{d-2}{\|y\|^2}$.

11. Implications

Stein’s paradox shows something basic about estimation in high dimensions. Treating parameters one at a time is worse than it should be even when they really are unrelated and borrowing strength across the ensemble by shrinking toward a common value cuts total risk because it trades a small amount of bias for a big cut in variance.

In modern machine learning this idea shows up everywhere. Regularization and dropout and weight decay are all forms of shrinkage that bias models toward simpler solutions because high-dimensional parameter spaces blow up estimation error.

Background

Terminology

An estimator is admissible if no other estimator gets lower risk $R(\theta, \delta)$ for every $\theta$ at the same time. When estimator A gets $R(\theta, A) \le R(\theta, B)$ for every $\theta$ and beats it strictly somewhere then A dominates B.

The James-Stein estimator is a shrinkage estimator that beats the MLE for the multivariate normal mean in dimensions $d \ge 3$. It works by pulling estimates toward a central value and trading bias for lower variance. The positive-part tweak clips the shrinkage factor so it stays non-negative and that stops overshrinkage.

Stein’s Lemma is the integration-by-parts identity $\mathbb{E}[(X - \theta)^T g(X)] = \mathbb{E}[\nabla \cdot g(X)]$ for Gaussian $X$ and it lets us work out risk without knowing $\theta$. The empirical Bayes framework estimates prior hyperparameters straight from data and the James-Stein estimator drops out as an empirical Bayes procedure. A minimax estimator minimizes worst-case risk across all $\theta$.

References

1. Stein, C. (1956). “Inadmissibility of the usual estimator for the mean of a multivariate normal distribution”. The original paper establishing the paradox.

2. James, W., & Stein, C. (1961). “Estimation with quadratic loss”. Constructs the explicit estimator and computes its risk.

3. Efron, B., & Morris, C. (1973). “Stein’s estimation rule and its competitors, an empirical Bayes approach”. Develops the empirical Bayes interpretation and demonstrates practical applications.

4. Lehmann, E. L., & Casella, G. (2006). “Theory of Point Estimation”. Standard graduate reference for decision theory and point estimation.