Chapter 18
25 min read
Section 119 of 175

Empirical Bayes

Bayesian Inference

Learning Objectives

By the end of this section, you will be able to:

📐 Core Mathematical Concepts

  • • Define Empirical Bayes and distinguish it from fully Bayesian methods
  • • Derive the marginal likelihood and explain hyperparameter estimation
  • • Explain the James-Stein paradox and its implications
  • • Calculate shrinkage estimators and their optimal shrinkage factors

🔧 Practical Skills

  • • Implement EB estimators for grouped data
  • • Apply shrinkage to improve predictions
  • • Use EB for multiple testing correction
  • • Choose when EB is appropriate vs. full Bayes

🧠 AI/ML Connections

  • • Hyperparameter Optimization: Evidence maximization (Type II ML) for learning regularization strength
  • • Transfer Learning: Using pretrained models as data-driven priors
  • • Mixed Effects Models: Random effect estimation in hierarchical models
  • • Sparse Coding: Automatic relevance determination for feature selection
Where You'll Apply This: Gene expression analysis, sports analytics (player rating), insurance premium estimation, meta-analysis combining multiple studies, false discovery rate control, recommendation systems with cold-start users, and any problem with multiple related estimation tasks.

The Big Picture: The Prior Problem

Full Bayesian inference requires specifying a prior distribution before seeing data. But where does the prior come from? In practice, this is often the hardest part of Bayesian analysis:

The Prior Problem
  • Where do hyperparameters come from?
  • How do we know the prior is "right"?
  • Results can be sensitive to prior choices
  • Critics accuse Bayesians of subjectivity
The EB Solution
  • Estimate hyperparameters from data
  • Prior is data-driven, not subjective
  • "Let the data speak" about the prior
  • Compromise between Bayesian and frequentist

Empirical Bayes (EB) offers an elegant solution: estimate the prior from the data itself. Instead of subjectively specifying hyperparameters, we use the marginal distribution of observations to infer what prior must have generated them.

The Empirical Bayes Philosophy

"When you have many similar problems, use the ensemble to learn the prior, then apply it to each individual problem."

Key Insight: EB requires multiple related estimation problems. If you only have one parameter to estimate, there's no way to learn the prior. But with many groups, schools, genes, players, or trials, the collective behavior reveals the underlying prior distribution.

Historical Development

📜
Herbert Robbins (1956)

Introduced the term "Empirical Bayes" and the foundational idea: use past data to construct a prior for future inference. His compound decision theory showed how to exploit the structure of multiple similar problems.

🎯
Charles Stein & Willard James (1961)

Discovered the James-Stein paradox: for estimating a multivariate normal mean in 3+ dimensions, shrinking all coordinates toward any fixed point dominates the MLE. This shocking result legitimized shrinkage estimation and connected to EB.

📊
Bradley Efron & Carl Morris (1973-1977)

Made EB practical with their famous baseball batting average analysis. Showed that EB estimates of early-season performance predict season-end averages better than raw averages. Developed parametric EB methods and connected them to shrinkage estimation.

🧬
Modern Era (1990s-Present)

Genomics revolutionized EB with massive multiple testing problems. Methods like limma, DESeq2, and Storey's q-value use EB for variance estimation and false discovery rate control. EB is now standard in high-dimensional statistics and machine learning.

Mathematical Framework

The Hierarchical Model

Empirical Bayes assumes a hierarchical structure where parameters themselves come from a distribution:

Level 2 (Prior):
θi∼iidG(η)\theta_i \overset{iid}{\sim} G(\eta)

Parameters are drawn from a prior GG with hyperparameters η\eta

Level 1 (Likelihood):
Xi∣θi∼F(θi)X_i | \theta_i \sim F(\theta_i)

Observations are generated conditionally on parameters

The key distinction from full Bayesian analysis:

ApproachHyperparametersPosterior
Full BayesSpecified subjectively OR given a hyperpriorP(theta|X, eta) with eta fixed or integrated
Empirical BayesEstimated from data: eta_hat = argmax p(X|eta)P(theta|X, eta_hat) with eta_hat plugged in

Marginal Likelihood

The marginal likelihood (or evidence) is obtained by integrating out the parameters:

p(X∣η)=∫p(X∣θ) p(θ∣η) dθp(X | \eta) = \int p(X | \theta) \, p(\theta | \eta) \, d\theta

The probability of data after integrating out the parameters, depending only on hyperparameters

Empirical Bayes estimates η\eta by maximizing this marginal likelihood:

η^=arg⁡max⁡η p(X∣η)=arg⁡max⁡η ∏i=1n∫p(Xi∣θi) p(θi∣η) dθi\hat{\eta} = \arg\max_\eta \, p(X | \eta) = \arg\max_\eta \, \prod_{i=1}^n \int p(X_i | \theta_i) \, p(\theta_i | \eta) \, d\theta_i
Also Called: This approach is known by many names: Type II Maximum Likelihood, Marginal Maximum Likelihood, Evidence Maximization, or Maximum Marginal Likelihood (MML). They all refer to the same idea: estimate hyperparameters by maximizing the marginal likelihood.

Interactive: Marginal Likelihood

See how the marginal likelihood varies with the prior variance hyperparameter. The optimal value balances model complexity with data fit.

Marginal Likelihood for Hyperparameter Estimation
Observations (5):[1.50, 2.30, 0.80, 3.10, 2.00]
0246810tau sq hat = 3.40Prior Variance (tau squared)Log Marginal Likelihood

Empirical Bayes Procedure: The marginal likelihood integrates out the parameter theta, leaving only the hyperparameter tau squared. We find the tau squared hat that maximizes this marginal likelihood, then use it as if it were the true prior variance. This "estimates the prior from the data" - the hallmark of Empirical Bayes.

Shrinkage Estimation

The most important consequence of Empirical Bayes is shrinkage: individual estimates are pulled toward a common value (usually the grand mean). This "borrowing of strength" across groups reduces overall estimation error.

Normal-Normal Shrinkage Formula

For Xi∣θi∼N(θi,σ2)X_i | \theta_i \sim N(\theta_i, \sigma^2) and θi∼N(μ,tau2)\theta_i \sim N(\mu, \\tau^2):

θ^iEB=B⏟shrinkage⋅μ^+(1−B)⋅Xi\hat{\theta}_i^{EB} = \underbrace{B}_{\text{shrinkage}} \cdot \hat{\mu} + (1 - B) \cdot X_i

where B=σ2σ2+tau2B = \frac{\sigma^2}{\sigma^2 + \\tau^2} is the shrinkage factor

Understanding the shrinkage factor B:

  • When tau2→0\\tau^2 \to 0 (no between-group variation): B → 1, complete shrinkage to mean
  • When σ2→0\sigma^2 \to 0 (no noise): B → 0, no shrinkage (MLE is optimal)
  • When tau2=σ2\\tau^2 = \sigma^2: B = 0.5, equal weight to prior and observation

The James-Stein Paradox

In 1961, Charles Stein proved one of statistics' most surprising results:

The James-Stein Theorem

For X∼N(θ,Ip)X \sim N(\theta, I_p) with p≥3p \geq 3, the MLE θ^=X\hat{\theta} = X is inadmissible.

θ^JS=(1−p−2∣∣X∣∣2)+X\hat{\theta}^{JS} = \left(1 - \frac{p-2}{||X||^2}\right)_+ X

has strictly lower risk than MLE for every possible θ\theta

Why is this shocking? The James-Stein estimator shrinks all coordinates toward zero. If the true mean is θ=(100,0,0)\theta = (100, 0, 0), shrinking the first coordinate toward zero seems absurd - yet the overall risk is still lower than MLE!

The Resolution: The paradox is resolved by understanding that we're evaluating total risk across all coordinates. The improvement in estimating near-zero coordinates can outweigh the harm done to coordinates far from zero. This only works in 3+ dimensions because we need enough "borrowing" to offset the losses.

Interactive: James-Stein Demo

Observe the James-Stein paradox in action. Run multiple trials to see how the shrinkage estimator consistently achieves lower total squared error than the MLE.

James-Stein Paradox Demonstration
NaNNaNNaNNaNNaNMLEJSTrial NumberSquared Error Loss
NaN
MLE Average Loss
Expected: 5.0
NaN
James-Stein Average Loss
NaN% improvement

The Paradox: For p ≥ 3 dimensions, the James-Stein estimator (which shrinks all coordinates toward zero) has strictly lower risk than the MLE for every possible true mean vector. This was shocking when discovered - how can shrinking toward an arbitrary point (zero) help estimate means that may be far from zero?

Interactive: EB Shrinkage Explorer

Visualize how Empirical Bayes shrinks individual group estimates toward the grand mean. Each group's MLE (red) is pulled toward the center, producing the EB estimate (blue). Compare both to the true values (green).

Empirical Bayes Shrinkage Explorer
-3-2-10123Grand Mean123456789101112131415Parameter ValueGroup
MLE (Observed)
EB Estimate
True Value
31%
Shrinkage Factor
1.546
MLE MSE
1.089
EB MSE
30% reduction

Observe: The EB estimates (blue squares) are shrunk toward the grand mean relative to the MLE (red circles). This shrinkage reduces overall MSE by borrowing strength across groups.

EB vs MLE vs Full Bayes

Understanding when to use Empirical Bayes requires comparing it to the alternatives:

MLE

  • Uses only individual-level data
  • No borrowing across groups
  • Optimal with infinite data
  • Can overfit with small samples

Use when: Lots of data per group, groups are truly unrelated

Empirical Bayes

  • Estimates prior from data
  • Borrows strength across groups
  • Data-driven, less subjective
  • Underestimates uncertainty

Use when: Many related groups, moderate data per group

Full Bayes

  • Specifies (or puts prior on) hyperparameters
  • Full uncertainty quantification
  • Requires prior specification
  • Computationally intensive (MCMC)

Use when: Uncertainty matters, have prior knowledge, need proper intervals

Interactive: Method Comparison

MLE vs Full Bayes vs Empirical Bayes
-4-2024MLEBayes/EBParameter Value (theta)
Prior
Likelihood
Posterior
3.00
MLE
Uses only data
2.40
Empirical Bayes
Data-driven prior
0.20
Shrinkage Factor
toward prior mean

Key insight: EB uses the same formula as Full Bayes but estimates hyperparameters from data instead of specifying them subjectively. The shrinkage factor determines how much weight is given to the prior vs. the data.

Real-World Applications

AI/ML Applications

Empirical Bayes ideas permeate modern machine learning, often under different names:

🎛️ Hyperparameter Learning

Gaussian Process regression learns kernel hyperparameters by maximizing marginal likelihood - exactly Type II ML. Bayesian optimization uses this to balance exploration and exploitation. The "evidence" in evidence maximization is the marginal likelihood.

🧠 Neural Network Regularization

Learning the regularization strength via validation is implicitly EB. More sophisticated: Automatic Relevance Determination (ARD) learns a separate regularization weight per feature, enabling automatic feature selection.

📚 Transfer Learning

Using a pretrained model as initialization is EB thinking: the pretrained weights serve as a "prior" learned from a large dataset. Fine-tuning with small learning rates implements shrinkage toward this prior.

🎰 Multi-Armed Bandits

Thompson Sampling with EB-estimated priors: instead of specifying a prior on arm rewards, learn it from observed rewards across all arms. This accelerates exploration-exploitation in recommendation systems.

🔢 Matrix Factorization

Collaborative filtering with EB: learn the prior on user/item factors from the data. This addresses the cold-start problem - new users are shrunk toward population averages until enough data accumulates.

📊 Variational Autoencoders

VAEs learn the prior on latent variables jointly with the encoder/decoder. The standard N(0,I) prior can be replaced with a learned prior - a form of empirical Bayes on the latent space distribution.

Python Implementation

Here's a complete implementation of Empirical Bayes methods. Click on highlighted lines to see detailed explanations of the key concepts.

Empirical Bayes Implementation
🐍empirical_bayes.py
10Hierarchical Model Structure

Empirical Bayes uses a two-level hierarchical model: (1) Parameters come from a prior distribution with unknown hyperparameters. (2) Data is generated from parameters. We estimate hyperparameters from the marginal distribution of data.

21Known Within-Group Variance

In many applications, the observation variance is known or well-estimated (e.g., from replicate measurements). This is crucial for separating within-group noise from between-group variation.

29Marginal Maximum Likelihood

After integrating out the parameters theta_i, observations follow X_i ~ N(mu, sigma^2 + tau^2). We estimate (mu, tau^2) by maximizing this marginal likelihood.

39Method of Moments for Variance

The sample variance S^2 estimates sigma^2 + tau^2. Subtracting the known sigma^2 gives the between-group variance tau^2. We take max(0, ...) because variance cannot be negative.

45The Shrinkage Formula

The EB posterior mean shrinks observations toward the grand mean. The shrinkage factor B = sigma^2/(sigma^2 + tau^2) controls the amount of shrinkage. More noise (larger sigma^2) or less between-group variation (smaller tau^2) means more shrinkage.

55Complete Shrinkage Case

When tau^2 = 0, all groups are estimated to have the same mean (no real between-group variation). All observations are shrunk completely to the grand mean.

64James-Stein Estimator

The JS estimator shrinks all coordinates of X toward zero. The shrinkage factor (1 - (p-2)*sigma^2/||X||^2) depends on the dimension p and the observed norm. It dominates MLE for p >= 3.

74Positive Part Shrinkage

We take max(0, shrinkage) to avoid shrinking past zero. When ||X||^2 is small, the naive formula could give negative shrinkage (expansion), which is not sensible.

84Baseball Batting Averages

The classic Efron-Morris (1975) example: early-season batting averages are noisy estimates of true ability. EB shrinks extreme averages toward the league mean, improving prediction of season-end performance.

97Binomial to Normal Approximation

For batting averages, the binomial variance is p(1-p)/n. With average around 0.25 and n=45 at-bats, this gives sigma^2 around 0.25/45 for the normal approximation.

121Multiple Testing and FDR

In genomics, we test thousands of hypotheses. EB estimates the proportion of true nulls (pi_0) from the distribution of p-values. This is used in Benjamini-Hochberg and Storey's q-value method for controlling false discovery rate.

130Estimating Pi_0

Under null hypothesis, p-values are uniform. So the proportion of p-values > lambda estimates the proportion of nulls times the probability a uniform exceeds lambda. Solving gives the EB estimate of pi_0.

166 lines without explanation
1import numpy as np
2from scipy import stats
3from scipy.optimize import minimize_scalar
4import matplotlib.pyplot as plt
5
6# ============================================
7# CORE: Empirical Bayes Framework
8# ============================================
9
10class EmpiricalBayes:
11    """
12    Empirical Bayes estimator for normal means.
13
14    Hierarchical model:
15      theta_i ~ N(mu, tau^2)        # Prior
16      X_i | theta_i ~ N(theta_i, sigma^2)  # Likelihood
17
18    We estimate (mu, tau^2) from the marginal distribution of X,
19    then compute posterior means as point estimates.
20    """
21
22    def __init__(self, sigma_sq=1.0):
23        """Initialize with known within-group variance."""
24        self.sigma_sq = sigma_sq
25        self.mu_hat = None
26        self.tau_sq_hat = None
27
28    def estimate_hyperparameters(self, observations):
29        """
30        Estimate prior hyperparameters using marginal MLE.
31
32        Marginal distribution: X_i ~ N(mu, sigma^2 + tau^2)
33        """
34        n = len(observations)
35        X_bar = np.mean(observations)
36        S_sq = np.var(observations, ddof=1)
37
38        # MLE for mu (grand mean)
39        self.mu_hat = X_bar
40
41        # Method of moments estimate for tau^2
42        # E[S^2] = tau^2 + sigma^2, so tau^2 = S^2 - sigma^2
43        self.tau_sq_hat = max(0, S_sq - self.sigma_sq)
44
45        return self.mu_hat, self.tau_sq_hat
46
47    def posterior_mean(self, x_i):
48        """
49        Compute EB posterior mean (shrinkage estimator).
50
51        E[theta_i | X_i] = B * mu + (1-B) * X_i
52        where B = sigma^2 / (sigma^2 + tau^2) is shrinkage factor
53        """
54        if self.tau_sq_hat == 0:
55            return self.mu_hat  # Complete shrinkage if no between-group variance
56
57        B = self.sigma_sq / (self.sigma_sq + self.tau_sq_hat)
58        return B * self.mu_hat + (1 - B) * x_i
59
60    def shrinkage_factor(self):
61        """Return the estimated shrinkage factor B."""
62        if self.tau_sq_hat == 0:
63            return 1.0
64        return self.sigma_sq / (self.sigma_sq + self.tau_sq_hat)
65
66
67def james_stein_estimator(X, sigma_sq=1.0):
68    """
69    James-Stein estimator for multivariate normal mean.
70
71    X ~ N(theta, sigma^2 * I_p)
72
73    JS estimator: theta_hat = (1 - (p-2)*sigma^2 / ||X||^2)_+ * X
74
75    Dominates MLE for p >= 3 dimensions!
76    """
77    p = len(X)
78    if p < 3:
79        return X  # JS not applicable for p < 3
80
81    X_norm_sq = np.sum(X**2)
82
83    # Shrinkage factor (positive part)
84    shrinkage = max(0, 1 - (p - 2) * sigma_sq / X_norm_sq)
85
86    return shrinkage * X
87
88
89# ============================================
90# Example 1: Baseball Batting Averages
91# ============================================
92
93# Simulated batting averages (similar to Efron-Morris 1975)
94np.random.seed(42)
95n_players = 18
96true_abilities = np.random.beta(80, 220, n_players)  # True abilities
97at_bats = 45  # Each player has 45 at-bats in early season
98
99# Observed batting averages (binomial -> normal approximation)
100hits = np.random.binomial(at_bats, true_abilities)
101observed_avg = hits / at_bats
102
103# Within-player variance (binomial variance / n)
104sigma_sq = 0.25 / at_bats  # Approx variance of batting average
105
106# Fit Empirical Bayes
107eb = EmpiricalBayes(sigma_sq=sigma_sq)
108mu_hat, tau_sq_hat = eb.estimate_hyperparameters(observed_avg)
109
110print("=== Baseball Batting Averages (Efron-Morris) ===")
111print(f"Grand mean (mu_hat): {mu_hat:.4f}")
112print(f"Between-player variance (tau^2_hat): {tau_sq_hat:.6f}")
113print(f"Shrinkage factor: {eb.shrinkage_factor():.3f}")
114
115# Compute EB estimates
116eb_estimates = np.array([eb.posterior_mean(x) for x in observed_avg])
117
118# Compare MSE
119mle_mse = np.mean((observed_avg - true_abilities)**2)
120eb_mse = np.mean((eb_estimates - true_abilities)**2)
121
122print(f"\nMLE MSE:  {mle_mse:.6f}")
123print(f"EB MSE:   {eb_mse:.6f}")
124print(f"MSE reduction: {100*(1 - eb_mse/mle_mse):.1f}%")
125
126# ============================================
127# Example 2: James-Stein for Normal Means
128# ============================================
129
130print("\n=== James-Stein Estimator ===")
131p = 10  # Dimension
132theta_true = np.array([1, 2, 0, -1, 3, 0.5, -2, 1.5, 0, 2.5])
133X = theta_true + np.random.normal(0, 1, p)
134
135mle = X
136js = james_stein_estimator(X, sigma_sq=1.0)
137
138mle_loss = np.sum((mle - theta_true)**2)
139js_loss = np.sum((js - theta_true)**2)
140
141print(f"MLE squared error loss: {mle_loss:.3f}")
142print(f"JS squared error loss:  {js_loss:.3f}")
143print(f"JS improvement: {100*(1 - js_loss/mle_loss):.1f}%")
144
145# ============================================
146# Example 3: Multiple Testing (FDR estimation)
147# ============================================
148
149def estimate_pi0_empirical_bayes(p_values, lambda_val=0.5):
150    """
151    Estimate proportion of true nulls (pi_0) using EB.
152
153    Under null: p ~ Uniform(0,1)
154    Under alternative: p tends to be small
155
156    pi_0 estimate = #{p > lambda} / (n * (1 - lambda))
157    """
158    n = len(p_values)
159    n_above_lambda = np.sum(p_values > lambda_val)
160    pi_0_hat = n_above_lambda / (n * (1 - lambda_val))
161    return min(1.0, pi_0_hat)  # Cap at 1
162
163# Simulate multiple testing scenario
164n_tests = 1000
165true_nulls = np.random.choice([True, False], n_tests, p=[0.9, 0.1])
166
167p_values = np.where(
168    true_nulls,
169    np.random.uniform(0, 1, n_tests),  # Null: uniform
170    np.random.beta(1, 10, n_tests)      # Alternative: skewed toward 0
171)
172
173pi_0_hat = estimate_pi0_empirical_bayes(p_values)
174pi_0_true = np.mean(true_nulls)
175
176print("\n=== Multiple Testing (Empirical Bayes FDR) ===")
177print(f"True null proportion: {pi_0_true:.3f}")
178print(f"EB estimated pi_0:    {pi_0_hat:.3f}")

Common Pitfalls

❌

Underestimating Uncertainty

EB treats estimated hyperparameters as known truth, ignoring uncertainty in their estimation. This makes credible intervals too narrow. For proper uncertainty, use full Bayes with hyperpriors or bootstrap.

Fix: For critical uncertainty quantification, use hierarchical Bayes with MCMC to integrate over hyperparameter uncertainty.

❌

Using EB with Too Few Groups

With only 3-5 groups, the hyperparameter estimates are unstable. EB needs many groups to reliably estimate the prior. With few groups, results can be worse than MLE.

Rule of thumb: EB works best with 10+ groups. With fewer, consider using fixed informative priors based on domain knowledge.

❌

Ignoring Model Misspecification

Parametric EB assumes a specific prior family (usually normal). If the true prior is multimodal or heavy-tailed, shrinkage toward a single point can hurt.

Fix: Use nonparametric EB methods that estimate the prior distribution flexibly, or check residuals for non-normality.

❌

Circular Use of Data

Using the same data to estimate hyperparameters and then to compute individual posteriors can lead to over-shrinkage and overfitting to noise in extreme observations.

Mitigation: Leave-one-out cross-validation or regularized estimates of hyperparameters can help. The effect is usually minor with many groups.

💡

When EB Excels

EB is ideal when: (1) You have many related estimation problems, (2) Individual sample sizes are small to moderate, (3) You care about prediction accuracy more than uncertainty quantification, (4) You have no strong prior knowledge to specify hyperparameters.

Knowledge Check

Test your understanding of Empirical Bayes with this interactive quiz.

Knowledge Check1 of 5

What distinguishes Empirical Bayes from fully Bayesian inference?

Summary

Key Takeaways

  1. EB estimates the prior from data: Instead of subjectively specifying hyperparameters, EB uses the marginal likelihood to learn them. This addresses the "prior problem" in Bayesian inference.
  2. Shrinkage improves overall estimation: EB shrinks individual estimates toward a common value, borrowing strength across groups. This reduces total MSE even when some individual estimates get worse.
  3. The James-Stein paradox: For p ≥ 3 dimensions, shrinking toward any point dominates MLE. This counterintuitive result motivated the development of shrinkage estimation.
  4. EB requires multiple related problems: You need many groups/units to estimate the prior. With too few groups, hyperparameter estimates are unstable.
  5. EB underestimates uncertainty: By treating hyperparameter estimates as fixed, EB produces intervals that are too narrow. For proper uncertainty, use full hierarchical Bayes.
  6. Wide applications in ML: Hyperparameter learning, transfer learning, automatic relevance determination, and variance estimation in genomics all use EB ideas.
Chapter Complete: You've now covered the full spectrum of Bayesian Inference - from point estimation (MAP) to interval estimation (credible intervals), model comparison (Bayes factors), and now the data-driven approach of Empirical Bayes. In the next chapter, we'll explore Computational Bayesian Methods - MCMC, Metropolis-Hastings, and variational inference - the algorithms that make Bayesian inference practical for complex models.
Loading comments...
PreviousTable of ContentsNext