Statistical Decision Theory

Statistical decision theory provides a framework for choosing between actions (estimators, tests, decisions) by formalizing the consequences of each choice through loss functions and risk.

The Decision Problem

Loss Functions

Risk Function

Admissibility

Bayes Risk

Minimax Estimators

The James-Stein Estimator

Pareto Optimality

Python Implementation

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

np.random.seed(42)

# --- Risk calculation for normal means problem ---
def squared_error_risk(theta_hat, theta_true):
    return np.mean((theta_hat - theta_true) ** 2)

def james_stein_estimator(y, sigma2):
    J = len(y)
    shrinkage = max(0, 1 - (J - 2) * sigma2 / np.sum(y**2))
    return shrinkage * y

def simulate_risk(true_theta, sigma2=1.0, n_sims=5000):
    J = len(true_theta)
    y = np.random.randn(n_sims, J) * np.sqrt(sigma2) + true_theta
    
    risk_mle = np.zeros(n_sims)
    risk_js = np.zeros(n_sims)
    
    for i in range(n_sims):
        risk_mle[i] = squared_error_risk(y[i], true_theta)
        risk_js[i] = squared_error_risk(james_stein_estimator(y[i], sigma2), true_theta)
    
    return np.mean(risk_mle), np.mean(risk_js)

# --- Scenario 1: Small theta (near zero) ---
J = 10
theta_small = np.random.randn(J) * 0.3
risk_mle_small, risk_js_small = simulate_risk(theta_small)
print(f"Small θ: MLE risk = {risk_mle_small:.3f}, JS risk = {risk_js_small:.3f}")

# --- Scenario 2: Large theta (far from zero) ---
theta_large = np.random.randn(J) * 3.0
risk_mle_large, risk_js_large = simulate_risk(theta_large)
print(f"Large θ: MLE risk = {risk_mle_large:.3f}, JS risk = {risk_js_large:.3f}")

# --- Risk as function of ||θ||² ---
norms = np.linspace(0.1, 50, 100)
risks_mle = []
risks_js = []
for norm_sq in norms:
    theta = np.random.randn(J) * np.sqrt(norm_sq / J)
    rm, rj = simulate_risk(theta)
    risks_mle.append(rm)
    risks_js.append(rj)

fig, axes = plt.subplots(1, 3, figsize=(16, 5))

axes[0].plot(norms, risks_mle, 'b-', linewidth=2, label='MLE (y)')
axes[0].plot(norms, risks_js, 'r-', linewidth=2, label='James-Stein')
axes[0].axhline(J, color='blue', linestyle='--', alpha=0.5, label=f'J={J}')
axes[0].set_xlabel('||θ||²')
axes[0].set_ylabel('Risk (MSE)')
axes[0].set_title(f'Risk Comparison (J={J})')
axes[0].legend()

# --- Loss function comparison ---
theta_range = np.linspace(-3, 3, 200)
axes[1].plot(theta_range, theta_range**2, 'b-', linewidth=2, label='Squared error: $(θ-a)²$')
axes[1].plot(theta_range, np.abs(theta_range), 'r-', linewidth=2, label='Absolute error: $|θ-a|$')
axes[1].plot(theta_range, (theta_range != 0).astype(float), 'g-', linewidth=2, label='0-1 loss: $1(θ≠a)$')
axes[1].set_xlabel('θ - a (estimation error)')
axes[1].set_ylabel('Loss')
axes[1].set_title('Loss Functions')
axes[1].legend()

# --- Bias-variance tradeoff ---
lambdas = np.linspace(0, 5, 100)
bias_sq = (lambdas * 0.5)**2  # Squared bias increases with λ
variance = 1.0 / (1 + lambdas)  # Variance decreases with λ
total_risk = bias_sq + variance

axes[2].plot(lambdas, bias_sq, 'r--', linewidth=2, label='Bias²')
axes[2].plot(lambdas, variance, 'b--', linewidth=2, label='Variance')
axes[2].plot(lambdas, total_risk, 'k-', linewidth=2, label='Total risk')
min_idx = np.argmin(total_risk)
axes[2].axvline(lambdas[min_idx], color='green', linestyle=':', alpha=0.7, 
                label=f'Optimal λ={lambdas[min_idx]:.2f}')
axes[2].set_xlabel('Regularization parameter λ')
axes[2].set_ylabel('Risk')
axes[2].set_title('Bias-Variance Tradeoff')
axes[2].legend()

plt.tight_layout()
plt.savefig('decision_theory.png', dpi=150)
plt.show()

# --- Minimax vs Bayes comparison ---
print("\n=== Minimax vs Bayes ===")
theta_grid = np.linspace(-5, 5, 200)
risk_grid_mle = theta_grid**2  # Risk of MLE (constant = J)
# Bayes rule with N(0, τ²) prior
tau2 = 2.0
bayes_shrink = tau2 / (tau2 + 1.0)
risk_grid_bayes = bayes_shrink**2 * theta_grid**2 + (1 - bayes_shrink**2)

print(f"Minimax value (MLE risk): {J}")
print(f"Bayes risk (uniform prior): {np.mean(risk_grid_bayes):.3f}")

Related Topics

• See Bayesian Linear Regression for Bayes rules applied to regression
• See Ridge Regression for shrinkage in practice
• See Hypothesis Testing for decision-theoretic view of hypothesis tests

Key Takeaways