Maximum Likelihood Estimation

Learning Objectives

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

Where You'll Apply This: Every neural network training, logistic regression, Naive Bayes classifiers, language models (GPT, BERT), VAEs, normalizing flows, and essentially ALL probabilistic machine learning.

The Big Picture: A Historical Journey

It's 1912. You're Ronald Aylmer Fisher, a young mathematician at Cambridge University, frustrated with the ad-hoc methods statisticians use to estimate parameters. Different researchers use different methods, with no clear way to compare them. You ask: Is there a principled, unified approach?

Fisher's Revolutionary Insight

Consider a fundamental shift in perspective. Instead of asking "what is the probability of the data given parameters?" (the forward problem), Fisher asked:

This seemingly simple question revolutionized statistics and, a century later, forms the mathematical foundation for training every neural network.

The Core Intuition

Imagine you flip a coin 10 times and get 7 heads. What's the "best guess" for the coin's bias (probability of heads)?

1. If the coin is fair (p = 0.5), getting 7 heads has probability $\binom{10}{7}(0.5)^{10} \approx 0.117$
2. If p = 0.7, getting 7 heads has probability $\binom{10}{7}(0.7)^7(0.3)^3 \approx 0.267$
3. If p = 0.9, getting 7 heads has probability $\binom{10}{7}(0.9)^7(0.1)^3 \approx 0.057$

The value p = 0.7 makes our observed data (7 heads in 10 flips) most probable. This is the Maximum Likelihood Estimate!

Interactive: Likelihood Explorer

Explore how the likelihood function changes as you vary parameters. Drag the slider to find the MLE - the parameter value at the peak of the curve.

Mathematical Foundation

The Likelihood Function

Given i.i.d. observations $X_1, X_2, \ldots, X_n$ from a distribution with parameter $\theta$, the likelihood function is:

The Log-Likelihood Trick

In practice, we almost always work with the log-likelihood:

Why log-likelihood?

1. Numerical stability: Products of small probabilities cause underflow (10⁻¹⁰⁰). Sums of log-probabilities stay manageable (-230.5).
2. Same maximum: Since log is monotonically increasing, maximizing ℓ(θ) gives the same θ̂ as maximizing L(θ).
3. Easier calculus: Derivatives of sums are simpler than derivatives of products.
4. Information theory: Log-probabilities connect to entropy and information measures.

The MLE Recipe

Finding the MLE follows a systematic procedure:

1. Write the likelihood function:
   $L(\theta) = \prod_{i=1}^{n} f(x_i | \theta)$
2. Take the logarithm:
   $\ell(\theta) = \sum_{i=1}^{n} \log f(x_i | \theta)$
3. Differentiate (score function):
   $s(\theta) = \frac{d\ell}{d\theta} = \sum_{i=1}^{n} \frac{\partial}{\partial\theta} \log f(x_i | \theta)$
4. Set to zero and solve:
   $s(\hat{\theta}) = 0$
5. Verify maximum: Check that $\frac{d^2\ell}{d\theta^2}\bigg|_{\hat{\theta}} < 0$

Worked Examples

Numerical MLE: When Calculus Fails

Many distributions don't have closed-form MLEs. When we can't solve $\nabla \ell(\theta) = 0$analytically, we use numerical optimization - the same algorithms that train neural networks!

Interactive: Numerical Optimizer

Watch optimization algorithms climb the log-likelihood surface. The Gamma distribution (estimating α with known β) has no closed-form MLE - perfect for demonstrating numerical methods.

MLE in Deep Learning

Here's the profound connection that every AI/ML practitioner should understand: Training neural networks with standard loss functions IS maximum likelihood estimation!

Cross-Entropy = Classification MLE

When you train a classifier with cross-entropy loss, you are performing MLE for a categorical (multinomial) distribution over class labels.

# This code performs Maximum Likelihood Estimation!
import torch
import torch.nn as nn

model = nn.Linear(2, 1)  # θ = (w1, w2, bias)
sigmoid = nn.Sigmoid()
criterion = nn.BCELoss()  # Binary Cross-Entropy = -log P(y|x,θ)

# Training loop
for epoch in range(50):
    probs = sigmoid(model(X))     # P(y=1|x,θ)
    loss = criterion(probs, y)    # = -Σ log P(yᵢ|xᵢ,θ) / n

    # Gradient ascent on log-likelihood (descent on -log-likelihood)
    optimizer.zero_grad()
    loss.backward()               # ∇θ[-ℓ(θ)]
    optimizer.step()              # θ ← θ - α·∇[-ℓ(θ)] = θ + α·∇ℓ(θ)

MSE = Regression MLE

Similarly, MSE loss corresponds to MLE assuming Gaussian noise in your regression model. If $y = f(x; \theta) + \varepsilon$ where $\varepsilon \sim N(0, \sigma^2)$:

Real-World Applications

MLE vs Method of Moments

Properties of MLE

Python Implementation

1import numpy as np
2from scipy.optimize import minimize
3from scipy import stats
4
5# ============================================
6# Closed-Form MLEs
7# ============================================
8
9def mle_normal(data):
10    """MLE for Normal distribution parameters.
11
12    Returns
13    -------
14    mu_hat : float - estimated mean
15    sigma_hat : float - estimated std (biased, uses n)
16    """
17    n = len(data)
18    mu_hat = np.mean(data)
19    # MLE uses n (biased), not n-1 (unbiased)
20    sigma2_hat = np.sum((data - mu_hat)**2) / n
21    return mu_hat, np.sqrt(sigma2_hat)
22
23
24def mle_bernoulli(data):
25    """MLE for Bernoulli probability p."""
26    return np.mean(data)  # Sample proportion
27
28
29def mle_exponential(data):
30    """MLE for Exponential rate parameter lambda."""
31    return 1.0 / np.mean(data)
32
33
34def mle_poisson(data):
35    """MLE for Poisson rate parameter lambda."""
36    return np.mean(data)
37
38
39# ============================================
40# Numerical MLE (when closed-form doesn't exist)
41# ============================================
42
43def mle_gamma_numerical(data):
44    """MLE for Gamma(alpha, beta) via numerical optimization."""
45    def neg_log_likelihood(params):
46        alpha, beta = params
47        if alpha <= 0 or beta <= 0:
48            return np.inf
49        # Sum of log PDFs
50        return -np.sum(stats.gamma.logpdf(data, a=alpha, scale=1/beta))
51
52    # Initialize with Method of Moments
53    mean_x = np.mean(data)
54    var_x = np.var(data)
55    alpha_init = mean_x**2 / var_x
56    beta_init = mean_x / var_x
57
58    result = minimize(
59        neg_log_likelihood,
60        x0=[alpha_init, beta_init],
61        method='L-BFGS-B',
62        bounds=[(1e-6, None), (1e-6, None)]
63    )
64    return result.x  # [alpha_hat, beta_hat]
65
66
67# ============================================
68# MLE = Deep Learning Training!
69# ============================================
70
71# PyTorch example (pseudocode)
72"""
73import torch
74import torch.nn as nn
75
76# Classification: Cross-Entropy = Negative Log-Likelihood
77model = nn.Sequential(nn.Linear(10, 32), nn.ReLU(), nn.Linear(32, 2))
78criterion = nn.CrossEntropyLoss()  # This IS -log P(y|x,theta)!
79
80for epoch in range(100):
81    logits = model(X)
82    loss = criterion(logits, y)  # Negative log-likelihood
83
84    optimizer.zero_grad()
85    loss.backward()  # Gradient of -log L
86    optimizer.step()  # Gradient ascent on log-likelihood
87
88
89# Regression: MSE = MLE with Gaussian noise assumption
90model = nn.Linear(10, 1)
91criterion = nn.MSELoss()  # Proportional to -log P(y|x,theta) for Gaussian
92
93for epoch in range(100):
94    predictions = model(X)
95    loss = criterion(predictions, y)
96
97    optimizer.zero_grad()
98    loss.backward()
99    optimizer.step()
100"""
101
102
103# ============================================
104# Example Usage
105# ============================================
106
107if __name__ == "__main__":
108    np.random.seed(42)
109
110    # Normal MLE
111    data_normal = np.random.normal(loc=5, scale=2, size=100)
112    mu_hat, sigma_hat = mle_normal(data_normal)
113    print(f"Normal MLE: μ̂={mu_hat:.3f} (true=5), σ̂={sigma_hat:.3f} (true=2)")
114
115    # Exponential MLE
116    data_exp = np.random.exponential(scale=2, size=100)  # scale = 1/lambda
117    lambda_hat = mle_exponential(data_exp)
118    print(f"Exponential MLE: λ̂={lambda_hat:.3f} (true=0.5)")
119
120    # Gamma MLE (numerical)
121    data_gamma = np.random.gamma(shape=3, scale=2, size=200)  # scale = 1/beta
122    alpha_hat, beta_hat = mle_gamma_numerical(data_gamma)
123    print(f"Gamma MLE: α̂={alpha_hat:.3f} (true=3), β̂={beta_hat:.3f} (true=0.5)")

Common Issues and Debugging

Knowledge Check

Test your understanding of Maximum Likelihood Estimation with this interactive quiz.

Summary

Looking Ahead: In the next section, we'll explore the deeper properties of MLE - including Fisher Information, the Cramér-Rao Lower Bound, and what makes MLE so special among all estimators.