Stein's Paradox & Empirical Bayes

1. The Singular Geometry of $d \ge 3$

The Motivation: In dimension $d=1$, observing $X \sim \mathcal{N}(\theta, 1)$, the intuitive best guess for $\theta$ is $X$. In dimension $d=100$, observing $X \sim \mathcal{N}(\theta, I)$, the guess $\hat{\theta} = X$ is catastrophically wrong. Specifically, for any fixed $\theta$, as $d \to \infty$:

However, the squared norm of the estimator itself behaves pathologically:

Taking expectations:

The estimator $X$ is systematically “too long”. It places the estimate on a shell of radius $\sqrt{\|\theta\|^2 + d}$, significantly further from the origin than the true parameter $\theta$. In high dimensions, volume concentrates on the surface of the hypersphere. The independent noise vectors $Z_i$ accumulate Euclidean length orthogonal to $\theta$. To reduce risk, we must shrink the vector $X$ back towards the origin to compensate for this phantom length.

The Object: We consider the problem of estimating a mean vector $\theta \in \mathbb{R}^d$ from a single observation $X \sim \mathcal{N}_d(\theta, I)$. Is the sample mean $\delta_0(X) = X$ admissible under squared error loss? Stein (1956) proved: $d \ge 3 \implies \delta_0 \text{ is Inadmissible.}$

2. Definitions and The Stein Effect

Let the risk function be the expected squared error:

For the standard estimator $\delta_0(X) = X$:

We define the James-Stein Estimator $\delta_{JS}$:

This estimator applies a non-linear scaling factor strictly less than 1 (for $d > 2$).

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

The improvement is most dramatic near the origin. If $\theta = 0$:

In Dimension $d=100$, the MLE has risk 100. James-Stein has risk 2. This is a 98% reduction in error “for free”.

The Geometric Intuition of Shrinkage: Consider the geometry. $X = \theta + \epsilon$. $\theta \cdot \epsilon$ is sum of zero-mean variables, so it is $O(\sqrt{d})$. $\|\epsilon\|^2$ is sum of $d$ squared Gaussians, so it is $\approx d$. Thus $X$ lies roughly on a sphere around $\theta$ of radius $\sqrt{d}$. However, from the perspective of the origin, $\|X\|^2 \approx \|\theta\|^2 + d$. The triangle inequality is strict in high dimensions because (random) high-dimensional vectors are nearly orthogonal to each other. We know strictly that $\|\theta\| < \|X\|$. It is logically incoherent to estimate $\theta$ with $X$ when we know, with probability approaching 1, that $\theta$ is strictly smaller in norm than $X$. The “optimal” shrinkage factor $c$ minimizes $\|\theta - cX\|^2$.

This heuristic aligns perfectly with the rigorous form derived below.

3. The Machinery: Stein’s Lemma

To prove dominance, we need a tool to compute the risk of a non-linear estimator $\delta(X) = X + g(X)$ without explicit integration against the unknown $\theta$. Charles Stein (1981) provided the Unbiased Risk Estimate (SURE).

Lemma (Stein’s Lemma): Let $X \sim \mathcal{N}_d(\theta, I)$. Let $g: \mathbb{R}^d \to \mathbb{R}^d$ be a weakly differentiable function such that $\mathbb{E} | \nabla \cdot g(X) | < \infty$. Then:

Proof (The Integration by Parts): We proceed component-wise. Consider the $i$-th component. Even in $d$-dimensions, the independence allows us to factor the density: $p(x) = \prod \phi(x_j - \theta_j)$.

Focus on the integration with respect to $x_i$:

Recall the property of the Gaussian kernel: $\phi'(z) = -z \phi(z)$. Thus $(x_i - \theta_i)\phi(x_i - \theta_i) = -\frac{\partial}{\partial x_i} \phi(x_i - \theta_i)$.

We integrate by parts with $u = g_i(x)$ and $dv = \phi'(x_i - \theta_i) dx_i$:

The Boundary Condition: For the boundary term to vanish, we require $g_i(x)$ to grow slower than $e^{x^2/2}$ at infinity. This creates a “Growth Condition” on the estimator. All polynomial estimators satisfy this. Assuming the boundary vanishes:

Summing over $i=1 \dots d$:

$\square$

4. Deriving the Dominator

We propose an estimator of the form $\delta(X) = X + g(X)$. The Risk is:

The first term is $d$ (MLE risk). The third term is handled by Stein’s Lemma: $2 \mathbb{E}[\nabla \cdot g(X)]$. Thus, we have an expression for Risk that does not depend on $\theta$:

This quantity $SURE(g) = \|g\|^2 + 2 \nabla \cdot g$ acts as an unbiased limit of the risk. To dominate the MLE, we need to find a non-zero $g$ such that the quantity in the bracket is strictly negative.

Candidate: The Inverse Shrinkage Let us test $g(X) = - \frac{c}{\|X\|^2} X$ for some constant $c$. Does this satisfy the regularity conditions? Yes, except at $X=0$, which is a set of measure zero for $d \ge 1$.

Step 4.1: Compute the Divergence We need $\nabla \cdot g(X) = \sum_i \partial_i g_i$.

Use the product rule:

Summing over $i$:

Step 4.2: Compute the Norm Squared

Step 4.3: Minimize the Risk Substitute back into the SURE formula:

This is a quadratic in $c$ of the form $Ac^2 - Bc$. The minimum occurs at $c^* = B / 2A$.

The value at the minimum is:

Step 4.4: The Conclusion For the improvement to be valid, we need $d-2 > 0 \implies d > 2$. Assuming $d \ge 3$, the expected risk change is strictly negative everywhere.

The term $\mathbb{E}[1/\|X\|^2]$ is large when $X$ is near 0. Thus, the improvement is massive when the true parameter $\theta$ is near 0. At $\theta=0$, $\|X\|^2 \sim \chi^2_d$. $\mathbb{E}[1/\chi^2_d] = 1/(d-2)$.

We have proven the paradox.

5. The Bayesian Perspective (Gaussian Prior)

The frequentist derivation looks like magic. Why does the estimator take the form $(1 - c/\|X\|^2)X$? It arises naturally from a hierarchical Bayesian model.

Level 1 (Data): $X | \theta \sim \mathcal{N}_d(\theta, I)$ Level 2 (Prior): $\theta \sim \mathcal{N}_d(0, A \cdot I)$

We seek the posterior mean $\mathbb{E}[\theta | X]$. Since the prior and likelihood are Gaussian, the posterior is Gaussian. The precision (inverse variance) adds:

The posterior mean is a precision-weighted average of the prior mean (0) and the data ($X$):

Let $B = \frac{1}{1+A}$. The optimal estimator is a linear shrinkage $(1-B)X$.

The Empirical Step: If we knew $A$ (the signal-to-noise ratio), we would use it. We don’t. However, consider the marginal distribution of $X$:

This implies that the squared norm $\|X\|^2$ scales with $1+A$:

We need to estimate the unknown shrinkage factor $B = \frac{1}{1+A}$. From the properties of the Chi-Square distribution:

Therefore:

The term $\frac{d-2}{\|X\|^2}$ is an unbiased estimator of the optimal frequentist shrinkage factor $B$. Plugging this estimate $\hat{B}$ into the Bayes rule gives precisely the James-Stein estimator.

Interpretation: James-Stein “learns” the prior variance $A$ from the data itself. If $\|X\|^2$ is large, it assumes $A$ is large (signal dominates noise) and shrinks little. If $\|X\|^2$ is small, it assumes $A$ is small (noise dominates signal) and shrinks aggressively.

6. The Geometry of Diffusions: Admissibility and Recurrence

There is a profound connection between Admissibility and the heat equation. Lawrence Brown (1971) proved that the admissibility of an estimator is equivalent to the recurrence of a related diffusion process.

Consider the generalized estimator $\delta(X) = X + \nabla \phi(X)$. Using Tweedie’s Formula, we can view any such estimator as a formal Bayes estimator against a prior measure $\pi(\theta) = e^{\phi(\theta)}$. The marginal distribution becomes $m(X)$. Ideally, $m(X)$ should be a superharmonic function.

The Laplacian Connection: The risk improvement of $\delta$ over $X$ is determined by the Laplacian of the square root of the marginal $g = \sqrt{m}$:

If $\sqrt{m}$ is superharmonic ($\Delta \sqrt{m} \le 0$), the risk is improved. Consider the prior $\pi(\theta) = 1$ (Flat prior). This yields $m(X) = 1$. $\Delta 1 = 0$. This is the boundary case (MLE).

Consider the Harmonic Prior $\pi(\theta) \propto \|\theta\|^{-(d-2)}$. The marginal is roughly $m(X) \propto \|X\|^{-(d-2)}$. Let $f(r) = r^{-(d-2)/2}$. We compute the Laplacian in spherical coordinates:

For $d \ge 3$, this function is strictly superharmonic.

Recurrence vs Transience: This algebraic fact reflects a topological one.

• In $d=1, 2$, Brownian motion is Recurrent. It visits every neighborhood infinitely often. The heat does not escape.
• In $d \ge 3$, Brownian motion is Transient. It wanders off to infinity.

The inadmissibility of the MLE in $d \ge 3$ corresponds to the “leakage” of probability mass to infinity. The James-Stein estimator essentially plugs this leak at infinity by enforcing a boundary condition that reflects the transient nature of the space. We shrink because the space is too vast to let the estimator wander freely.

7. The Positive Part Estimator

Standard JS can over-shrink. If $\|X\|^2$ is very small, $1 - \frac{d-2}{\|X\|^2}$ becomes negative. The vector flips sign. This is clearly wrong. If we observe $X$, the estimate should be in the same orthant. Baranchik’s Positive Part Estimator:

Theorem: $\delta_{JS+}$ dominates $\delta_{JS}$. Is $\delta_{JS+}$ admissible? No. It is non-analytic. (Kinks at the boundary). Admissible estimators for exponential families must be analytic (Brown, 1971). However, no simple analytic estimator is known to dominate JS+. The search for the Minimax Admissible estimator stays open practically.

8. The Exact Risk of the Positive Part Estimator

While the James-Stein estimator dominates the MLE, the Positive Part estimator $\delta_{JS+}$ dominates James-Stein.

Simulating the risk is easy, but can we derive the exact analytic expression for $R(\theta, \delta_{JS+})$? Yes, but it requires integrating over the non-central chi-square distribution.

Step 8.1: Decomposition The risk function can be written as:

The difference is non-zero only when shrinkage overshoots ($1 - \frac{d-2}{\|X\|^2} < 0$), i.e., $\|X\|^2 < d-2$. Inside this region:

• $\delta_{JS}(X) = (1 - \frac{d-2}{\|X\|^2}) X$
• $\delta_{JS+}(X) = 0$

Step 8.2: The Difference Term Let $A = \{ x : \|x\|^2 < d-2 \}$. The reduction in risk is:

Expand the first term:

Subtracting $\|\theta\|^2$:

This algebra is getting messy. Let’s use Stein’s Lemma in reverse. There is a cleaner identity (Kariya et al.):

The standard result states that the improvement is exactly:

Simpler path: $R_{JS+} = \mathbb{E} [ \| \hat{\theta}_{JS+} - \theta \|^2 ]$. Using the identity $X = \theta + Z$:

This expectation is over a Non-Central Chi-Square variable $Y = \|X\|^2 \sim \chi^2_d(\|\theta\|^2)$.

Step 8.3: Series Expansion The PDF of a non-central chi-square $f_Y(y; d, \lambda)$ is a Poisson mixture of central chi-squares:

where $\lambda = \|\theta\|^2$. We can compute the risk by term-wise integration:

This series converges rapidly. For the Positive Part correction, we need to integrate the chopped quadratic over the region $[0, d-2]$.

Since $2(d-2) - y - (d-2)^2/y = - \frac{1}{y} (y - (d-2))^2$, which is strictly negative. Our convention is reversed. The dominance implies $R_{JS+} < R_{JS}$. Calculating this integral analytically involves the Cumulative Distribution Function of the Central Chi-Square, $F_{\chi^2}$.

The integral of $(1 - \frac{d-2}{y})^2 y$ against the pdf $f_Y(y)$ over $[0, d-2]$ is strictly positive. This provides the rigorous justification for the red line in the simulation plot being strictly below the blue line.

9. Numerical Simulation: Risk vs Norm

Let’s visualize the risk function. We plot $R(\theta)$ as a function of $\|\theta\|$. For MLE, it’s a flat line at $d$. For JS, it starts low at $\|\theta\|=0$ and approaches $d$ asymptotically.

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

In 1975, Efron and Morris applied JS to the batting averages of 18 players. First 45 at-bats ($X$). Goal: Predict remainder of season ($p$). MLE prediction: Just use the first 45 average. JS prediction: Shrink towards the Grand Mean of all players. Result: The James-Stein estimator reduced the total squared prediction error by a factor of 3. The “Grand Mean” acts as the learned prior. Players with extreme lucky starts (e.g., hitting 450) were pulled down. Players with bad starts were pulled up. This is the essence of Regularization in Machine Learning (Ridge Regression). Ridge Regression is exactly Bayesian estimation with a Gaussian Prior centered at 0. James-Stein is Ridge Regression where the $\lambda$ is automatically tuned: $\lambda = \frac{d-2}{\|y\|^2}$.

11. Conclusion: The Price of Admissibility

Stein’s Paradox teaches us a counter-intuitive lesson about high-dimensional space: Isolation is expensive. Treating parameters independently forces us to ignore the shared structure that inevitably reduces variance. By borrowing strength from the ensemble—even when the ensemble members are seemingly unrelated—we protect ourselves from the extreme deviations that plague high-dimensional vectors. In modern Machine Learning, we don’t just estimate means; we train millions of parameters. Regularization, Dropout, and Weight Decay are all spiritual successors to Stein’s insight. We happily bias our models to buy back the stability that dimensions steal from us.

Historical Timeline

Appendix A: Glossary of Terms

• Admissible: An estimator that cannot be strictly dominated by another.
• Dominance: Estimator A dominates B if its Risk is lower for all $\theta$.
• Empirical Bayes: Using the data to estimate the prior hyperparameters.
• James-Stein Estimator: An explicit shrinkage estimator that dominates MLE in $d \ge 3$.
• Minimax: Minimizing the maximum possible risk.
• Positive Part: Truncating the shrinkage factor to be non-negative.
• Shrinkage: Pushing estimates towards a central value to reduce variance.
• Stein’s Lemma: An integration-by-parts identity for Gaussian expectations.

References

1. Stein, C. (1956). “Inadmissibility of the usual estimator for the mean of a multivariate normal distribution”. The bomb. The original paper that shocked the statistical world.

2. James, W., & Stein, C. (1961). “Estimation with quadratic loss”. Refined the estimator to the explicit form used today and computed the risk improvement.

3. Efron, B., & Morris, C. (1973). “Stein’s estimation rule and its competitors—an empirical Bayes approach”. Demystified the paradox by grounding it in Empirical Bayes. Made it palatable to practitioners.

4. Lehmann, E. L., & Casella, G. (2006). “Theory of Point Estimation”. The standard graduate text on statistical decision theory.