Stein’s Identity, Shrinkage, and Unbiased Risk Estimation
Estimating each coordinate of a Gaussian mean by the corresponding observation seems natural under squared error. In dimension at least three, however, shrinking the vector toward the origin can reduce total risk. Stein’s identity makes the comparison explicit: it replaces an unknown covariance term by a derivative. The same calculation leads to an unbiased estimate of risk and to the problem of choosing shrinkage weights.
The Gaussian sequence model observes \(Y_i=\theta_i+\sigma\varepsilon_i\), with independent standard normal errors. The loss is squared error. For the infinite-dimensional problem, summability and the parameter class matter; see the sequence-model setup.
We say an estimator \(\theta^*\) of a parameter \(\theta\) of a statistical model is inadmissible on \(\Theta \subset \mathbb R^d\) with respect to the squared risk if there exists another estimator \(\hat \theta\) such that
Furthermore, there exists some \(\theta_0 \in \Theta\) such that
Otherwise, we say the estimator \(\theta^*\) is admissible. Stein’s result shows that admissibility depends on the joint estimation problem: in the statistical model consisting of random vectors \(X_1, \dots, X_n\) where \(X_i = \theta + \eta_i\) for \(1 \leq i \leq n\), where \(\theta \in \mathbb R^d\) and \(\eta_i \overset{\rm i.i.d}{\sim} \mathcal{N}(0,I_d)\), then estimating \(\theta\) using the sample mean \(\overline X = \frac1n \sum_{1 \leq i \leq n} X_i\) when \(d \geq 3\) fails to be admissible: \(\overline X\) is inadmissible. Stein further proposed an estimator for \(\theta\) in this context, for which the risk is smaller than that of \(\overline X\) everywhere on \(\mathbb R^d\) if \(d \geq 3\). This estimator, as we will see, is based on a shrinkage of \(\overline X\) toward the origin with a shrinkage factor that depends on \(\|\overline X\|\). As is clear now, there is a potential for superefficiency phenomenon to occur in certain nonparametric and high-dimensional instances.
1. Stein’s Lemma
Let \(Y\sim N_d(\theta,\varepsilon^2I_d)\), with \(\varepsilon>0\). Suppose \(f\) is absolutely continuous on compact intervals in each coordinate for almost every fixed value of the others. A convenient sufficient integrability condition is \(\mathbb E_\theta f(Y)^2<\infty\) and \(\mathbb E_\theta|\partial_i f(Y)|<\infty\). Then
Here is a proof that justifies the boundary terms. First take \(\theta=0\), \(\varepsilon=1\), and let \(\chi_R\) be a smooth function of the \(i\)-th coordinate equal to one on \([-R,R]\), zero outside \([-2R,2R]\), bounded between zero and one, with \(|\chi_R'|\le C/R\). Integration by parts in that coordinate, followed by Fubini, gives
The Gaussian density has derivative \(-y_i\) times itself, and the cutoff removes boundary terms. The left side is dominated by \(|Y_i f(Y)|\), which is integrable by Cauchy–Schwarz. The first right-hand term converges by dominated convergence; the absolute value of the last is at most \((C/R)\mathbb E|f(Y)|\to0\). This proves the standard-normal identity.
For the general case put \(Z=(Y-\theta)/\varepsilon\) and \(\widetilde f(z)=f(\theta+\varepsilon z)\). The chain rule is \(\partial_i\widetilde f(z)=\varepsilon(\partial_i f)(\theta+\varepsilon z)\). Thus \(\mathbb E[(Y_i-\theta_i)f(Y)]=\varepsilon\mathbb E[Z_i\widetilde f(Z)]=\varepsilon^2\mathbb E[\partial_i f(Y)]\), as claimed. This sufficient version is all that the risk calculations below require. ◻
2. The Case of \(d \geq 3\)
For every \(\theta\in\mathbb R^d\) and \(d\ge3\),
The Gaussian density is bounded by \((2\pi\varepsilon^2)^{-d/2}\). On the unit ball the expectation is therefore bounded by a constant times \(\int_0^1r^{d-3}dr\), which is finite exactly when \(d>2\). Outside the ball, \(\|Y\|^{-2}\le1\). Positivity follows because \(Y\ne0\) almost surely. This separates the local singularity from the harmless tail; a bound with a positive quadratic exponential would instead diverge.
Rotational invariance also shows this expectation depends on \(\theta\) only through \(\|\theta\|\): for an orthogonal \(U\), \(U(Y-\theta)\) has the same centered Gaussian law and \(\|UY\|=\|Y\|\). At \(\theta=0\), radial integration gives \(\mathbb E\|Y\|^{-2}=1/[\varepsilon^2(d-2)]\), since integration by parts relates \(\int_0^\infty r^{d-1}e^{-r^2/2}dr\) to \((d-2)\int_0^\infty r^{d-3}e^{-r^2/2}dr\).
3. Stein’s Theorem
Write an estimator as \(\widehat\theta=Y+h(Y)\). Assume each \(h_i\) satisfies Stein's lemma and \(\mathbb E\|h(Y)\|^2<\infty\). Expanding the square and applying the identity coordinatewise gives
The term \(d\varepsilon^2\) is the risk of the unshrunk observation \(Y\). Consequently the expression inside the following display is Stein's unbiased risk estimate:
For \(\widehat\theta=g(Y)Y\), use \(h_i=(g-1)Y_i\) and \(\operatorname{div}h=d(g-1)+\sum_iY_i\partial_i g\). This yields
To prove the James–Stein improvement, take \(h_i(y)=-ay_i/\|y\|^2\), assigning any value at the null-probability point \(y=0\). For \(y\ne0\),
Each coordinate slice is smooth unless all other coordinates are zero, a Lebesgue-null set of slices. Moreover \(|\partial_i h_i|\le3|a|/\|y\|^2\) and \(\|h\|^2=a^2/\|y\|^2\); the previous section verifies the required integrability. Thus for \(d\ge3\),
Every \(0<a<2\varepsilon^2(d-2)\) strictly improves on \(Y\) at every \(\theta\). The coefficient is minimized at \(a=\varepsilon^2(d-2)\), giving
For \(Y=\overline X\) from \(n\) independent \(N_d(\theta,I_d)\) observations, \(\varepsilon^2=1/n\), so this becomes \((1-(d-2)/(n\|\overline X\|^2))\overline X\). These exact risk comparisons use the Gaussian location model and squared Euclidean loss. Medians and quantiles require their own distributional and risk arguments; normal approximations alone do not give the same finite-sample dominance theorem.
4. Example
The following archived baseball calculation illustrates shrinkage empirically. Its career batting averages are reference estimates, not known population parameters, and the sampled seasons also contribute to those averages. The observations have unequal binomial sampling variances, the common variance is estimated, and the shrinkage target is an estimated grand mean. Thus the printed error comparison is one descriptive calculation, not a verification of the exact Gaussian dominance theorem. Even in an equal-variance Gaussian model, shrinking the \(d-1\) dimensional residual about an estimated grand mean leads to a \(d-3\) coefficient; the archived code's \(N-2\) is a heuristic choice.
# real data ------------------
set.seed(422)
library(Lahman)
# 1. Select 30 players that don't have missing values and have a good amount of data (total hits > 500)
Batting <- as_tibble(Batting)
players <- Batting %>%
group_by(playerID) %>%
summarise(AB_total = sum(AB),
H_total = sum(H)) %>%
na.omit() %>%
filter(H_total>500) %>% # collect players with a lot of data
sample_n(size=30) %>% # randomly sample 30 players
pull(playerID) # `select` will produce a dataframe,
# `pull` gives a list which will be easier
# to work with when we `filter` by players
# 2. Get "truth"
TRUTH <- Batting %>%
filter(playerID %in% players) %>%
group_by(playerID) %>%
summarise(AB_total = sum(AB),
H_total = sum(H)) %>%
mutate(TRUTH = H_total/AB_total) %>%
dplyr::select(playerID,TRUTH)
# 3. Collect 5 observations per player
# 4.Predict truth with MLE and James-Stein Estimator
obs <- Batting %>%
filter(playerID %in% players) %>%
group_by(playerID) %>%
do(sample_n(., 5)) %>%
group_by(playerID) %>%
summarise(AB_total = sum(AB),
H_total = sum(H)) %>%
mutate(MLE = H_total/AB_total) %>%
dplyr::select(playerID,MLE,AB_total) %>%
inner_join(TRUTH,by="playerID")
# 5. Define variables for James-Stein estimator
p_=mean(obs$MLE)
N = length(obs$MLE)
truth_median = median(obs$AB_total)
df <- obs %>%
mutate(sigma2 = (p_*(1-p_))/truth_median,
JS=p_+(1-((N-2)*sigma2/(sum((MLE-p_)^2))))*(MLE-p_)
) %>%
dplyr::select(-AB_total,-sigma2)
head(df)
errors <- df %>%
mutate(mle_pred_error_i = (MLE-TRUTH)^2,
js_pred_error_i = (JS-TRUTH)^2) %>%
summarise(js_pred_error = sum(js_pred_error_i),
mle_pred_error = sum(mle_pred_error_i))
errors
# A tibble: 1 × 2
# js_pred_error mle_pred_error
# <dbl> <dbl>
#1 0.00200 0.002985. Unbiased Estimation of Risk
Return to \(y_j=\theta_j+\varepsilon\xi_j\), \(j\ge1\), with \(\theta\in\ell^2\). For a class \(\Lambda\subset\ell^2\), the oracle weights minimize \(R(\lambda,\theta)\) over \(\lambda\in\Lambda\), when a minimizer exists. They depend on unknown \(\theta\), so this is a benchmark rather than an estimator. Achieving oracle risk requires an additional theorem for the chosen class and selection rule.
For the first \(J\) coordinates and fixed deterministic weights, an unbiased risk estimate is
Indeed \(\mathbb E y_j^2=\theta_j^2+\varepsilon^2\), and \((1-\lambda_j)^2+2\lambda_j-1=\lambda_j^2\). This agrees with the SURE formula for a diagonal linear estimator.
One must not sum the unweighted noise terms over infinitely many coordinates. Instead subtract the common unknown quantity \(\|\theta\|^2\) from the target risk and estimate the differences:
For fixed \(\lambda\in\ell^2\), put \(c_j=\lambda_j^2-2\lambda_j\). This is bounded and square-summable. The mean series \(\sum_jc_j\theta_j^2\) converges absolutely. Independence and \(\operatorname{Var}(y_j^2)=2\varepsilon^4+4\varepsilon^2\theta_j^2\) imply
Thus the random series converges in mean square after centering, proving the stated expectation identity. A finite collection of candidate weights can be compared by minimizing \(\widehat U\). However, unbiasedness for each fixed candidate does not make the selected minimum an unbiased estimate of the selected estimator's risk. When weights depend smoothly on \(y\), the derivatives of those weights enter Stein's divergence; nonsmooth selection requires separate analysis. The oracle question therefore remains a selection-risk problem, rather than following from unbiasedness alone.