Probability and Statistics (6): Estimation — MLE, MAP, and the Bias-Variance Story

Everything we’ve built so far — distributions, expectations, limit theorems — assumed we knew the parameters. The Gaussian has mean $\mu$ and variance $\sigma^2$ . The Binomial has $n$ trials with success probability $p$ . But in practice, you don’t know $\mu$ or $p$ . You observe data and try to figure them out.

This is estimation theory: the bridge between probability (where parameters are given) and statistics (where parameters are inferred). It’s also where the foundations of machine learning live. Every time you train a model, you are estimating parameters from data. The quality of that estimation determines whether your model generalizes or overfits.

The Setup: Estimators vs Estimates#

Suppose we observe data $x_1, x_2, \ldots, x_n$ drawn i.i.d. from some distribution $p(x|\theta)$ , where $\theta$ is an unknown parameter (or vector of parameters).

An estimator $\hat{\theta}$ is a function of the data: $\hat{\theta} = g(X_1, \ldots, X_n)$ . It is a random variable (since the data are random). An estimate is the numerical value $\hat{\theta}(x_1, \ldots, x_n)$ obtained from a specific dataset.

Notation convention: $\hat{\theta}$ (hat) denotes an estimator/estimate. $\theta$ (no hat) denotes the true parameter.

Properties of Estimators#

Bias#

The bias of an estimator is

An estimator is unbiased if $\text{Bias}(\hat{\theta}) = 0$ , meaning $E[\hat{\theta}] = \theta$ for all values of $\theta$ .

So $E\left[\sum(X_i - \bar{X})^2\right] = n\sigma^2 + n\mu^2 - \sigma^2 - n\mu^2 = (n-1)\sigma^2$ . Dividing by $n$ gives $\frac{n-1}{n}\sigma^2$ . Dividing by $n-1$ gives $\sigma^2$ . $\blacksquare$

Consistency#

An estimator is consistent if $\hat{\theta}_n \xrightarrow{P} \theta$ as $n \to \infty$ .

By the WLLN, the sample mean $\bar{X}$ is consistent for $\mu$ . Both the biased ($1/n$ ) and unbiased ($1/(n-1)$ ) sample variances are consistent for $\sigma^2$ .

Efficiency#

Among all unbiased estimators of $\theta$ , the one with the smallest variance is the most efficient. The Cramer-Rao lower bound (derived below) tells us how small the variance can possibly be.

Sufficiency#

A statistic $T(X_1, \ldots, X_n)$ is sufficient for $\theta$ if the conditional distribution of the data given $T$ does not depend on $\theta$ . Informally, $T$ captures all the information in the data about $\theta$ — knowing $T$ , you can throw away the raw data without losing anything.

Example. For $X_i \sim \text{Bernoulli}(p)$ , the sum $T = \sum X_i$ is sufficient for $p$ . Knowing that you observed 7 successes out of 10 trials is all you need; the specific sequence (e.g., 1101110100 vs 1110101100) carries no additional information about $p$ .

Method of Moments#

The simplest estimation method: match sample moments to population moments.

Recipe:

1. Express the parameters in terms of population moments: $\theta = h(\mu_1, \mu_2, \ldots)$ where $\mu_k = E[X^k]$ .
2. Replace population moments with sample moments: $\hat{\mu}_k = \frac{1}{n}\sum X_i^k$ .
3. Plug in: $\hat{\theta}_{\text{MoM}} = h(\hat{\mu}_1, \hat{\mu}_2, \ldots)$ .

Example: Gamma distribution. $X \sim \text{Gamma}(\alpha, \beta)$ with $E[X] = \alpha/\beta$ and $E[X^2] = \alpha(\alpha+1)/\beta^2$ .

Method of moments estimators are easy to compute but generally not the most efficient. They serve as good starting points and as a sanity check.

Maximum Likelihood Estimation#

The Likelihood Function#

The log-likelihood is easier to work with (sums are simpler than products), and maximizing $\ell$ is equivalent to maximizing $L$ since $\ln$ is monotonically increasing.

The MLE#

The maximum likelihood estimator is the value of $\theta$ that maximizes the likelihood:

Example 1: Bernoulli#

The MLE is the sample proportion. Natural and intuitive.

Example 2: Gaussian (Both Parameters)#

Note: the MLE for $\sigma^2$ divides by $n$ , not $n-1$ . It is biased (by a factor of $(n-1)/n$ ) but consistent and efficient. For large $n$ , the difference is negligible.

Example 3: Poisson#

Again, the sample mean — which makes sense since $E[X] = \lambda$ for Poisson.

Example 4: Uniform#

where $x_{(n)} = \max(x_1, \ldots, x_n)$ . For $\theta \geq x_{(n)}$ , $L(\theta) = 1/\theta^n$ is decreasing in $\theta$ , so $L$ is maximized at $\hat{\theta}_{\text{MLE}} = x_{(n)} = \max_i x_i$ .

This is an interesting case where the MLE is biased: $E[X_{(n)}] = \frac{n}{n+1}\theta < \theta$ . The MLE consistently underestimates because the maximum of $n$ samples can never exceed $\theta$ but can be much less. The unbiased estimator is $\frac{n+1}{n} X_{(n)}$ .

This example also illustrates that MLE doesn’t always satisfy a smooth score equation — here the likelihood is discontinuous at $\theta = x_{(n)}$ .

Properties of the MLE#

Under regularity conditions (smooth, identifiable model):

1. Consistency: $\hat{\theta}_{\text{MLE}} \xrightarrow{P} \theta_0$ (the true parameter).
2. Asymptotic normality: $\sqrt{n}(\hat{\theta}_{\text{MLE}} - \theta_0) \xrightarrow{d} \mathcal{N}(0, I(\theta_0)^{-1})$ where $I(\theta_0)$ is the Fisher information.
3. Asymptotic efficiency: The MLE achieves the Cramer-Rao lower bound asymptotically.
4. Invariance: If $\hat{\theta}$ is the MLE of $\theta$ , then $g(\hat{\theta})$ is the MLE of $g(\theta)$ .

Proof of invariance. If $\hat{\theta} = \arg\max_\theta L(\theta)$ , and $\phi = g(\theta)$ is a one-to-one function, then $\hat{\phi} = g(\hat{\theta}) = \arg\max_\phi L(g^{-1}(\phi))$ . For non-injective $g$ , define $\hat{\phi} = g(\hat{\theta})$ directly; the profile likelihood is maximized at this value. $\blacksquare$

Example of invariance. If $\hat{\sigma}^2_{\text{MLE}}$ is the MLE of $\sigma^2$ , then $\hat{\sigma}_{\text{MLE}} = \sqrt{\hat{\sigma}^2_{\text{MLE}}}$ is the MLE of $\sigma$ . You don’t need to re-derive — just apply the transformation.

Warning: Invariance is a property of the MLE, not of unbiased estimators in general. If $\hat{\theta}$ is unbiased for $\theta$ , $g(\hat{\theta})$ is generally not unbiased for $g(\theta)$ (by Jensen’s inequality, unless $g$ is linear).

Fisher Information and the Cramer-Rao Bound#

Fisher Information#

The score function is $s(\theta) = \frac{\partial}{\partial\theta} \ln p(X|\theta)$ .

Under regularity conditions, $E[s(\theta)] = 0$ .

The second equality (derived via interchanging differentiation and expectation) gives a more convenient formula. For $n$ i.i.d. observations, $I_n(\theta) = n \cdot I_1(\theta)$ .

The Cramer-Rao Lower Bound#

Proof sketch. Apply the Cauchy-Schwarz inequality to $\text{Cov}(\hat{\theta}, s(\theta))$ , where $s$ is the total score. Since $E[s] = 0$ , the covariance equals $E[\hat{\theta} \cdot s]$ , which equals 1 (using the unbiasedness condition and differentiation under the integral). Then Cauchy-Schwarz gives $1 \leq \text{Var}(\hat{\theta}) \cdot \text{Var}(s) = \text{Var}(\hat{\theta}) \cdot I_n(\theta)$ . $\blacksquare$

Example. For Bernoulli, the CRLB is $\frac{p(1-p)}{n}$ . The MLE $\hat{p} = \bar{X}$ has $\text{Var}(\hat{p}) = \frac{p(1-p)}{n}$ , achieving the bound exactly. The MLE is efficient.

MAP Estimation#

From MLE to MAP#

MLE treats $\theta$ as a fixed unknown. Maximum a posteriori (MAP) estimation treats $\theta$ as a random variable with a prior distribution $p(\theta)$ .

Connection to Regularization#

The MAP estimate is the MLE with a penalty term — the log-prior acts as a regularizer.

Example. Gaussian prior on the mean of a Gaussian with known variance $\sigma^2$ .

This is a weighted average of the sample mean $\bar{X}$ and the prior mean $\mu_0$ . When $n$ is large, the data dominates and MAP approaches MLE. When $n$ is small, the prior pulls the estimate toward $\mu_0$ .

The Bias-Variance Decomposition#

The Mean Squared Error#

Connection to Machine Learning#

• High bias = underfitting. The model is too simple to capture the true pattern.
• High variance = overfitting. The model fits noise in the training data.
• The tradeoff: Increasing model complexity reduces bias but increases variance. The optimal model balances both.

Regularization (= MAP with a prior) deliberately introduces bias to reduce variance, often lowering the total MSE.

Python: Bias-Variance Visualization#

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

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

np.random.seed(42)

def f_true(x):
    return np.sin(2 * np.pi * x)

n_train = 20
n_datasets = 200
x_test = np.linspace(0, 1, 100)

fig, axes = plt.subplots(2, 3, figsize=(16, 10))

degrees = [1, 3, 5, 10, 15, 19]

for idx, degree in enumerate(degrees):
    ax = axes[idx // 3, idx % 3]
    predictions = np.zeros((n_datasets, len(x_test)))

    for d in range(n_datasets):
        x_train = np.random.uniform(0, 1, n_train)
        y_train = f_true(x_train) + np.random.normal(0, 0.3, n_train)
        coeffs = np.polyfit(x_train, y_train, degree)
        predictions[d] = np.polyval(coeffs, x_test)

    # Plot individual fits (first 20)
    for d in range(min(20, n_datasets)):
        ax.plot(x_test, predictions[d], 'b-', alpha=0.1, linewidth=0.5)

    # Plot mean prediction
    mean_pred = predictions.mean(axis=0)
    ax.plot(x_test, mean_pred, 'r-', linewidth=2, label='E[$\\hat{f}$]')

    # Plot true function
    ax.plot(x_test, f_true(x_test), 'k--', linewidth=2, label='f(x)')

    # Compute bias^2 and variance
    bias_sq = np.mean((mean_pred - f_true(x_test))**2)
    variance = np.mean(predictions.var(axis=0))

    ax.set_title(f'Degree {degree}\nBias$^2$={bias_sq:.3f}, Var={variance:.3f}',
                 fontsize=11)
    ax.set_ylim(-2, 2)
    ax.legend(fontsize=9)

plt.suptitle('Bias-Variance Tradeoff: Polynomial Regression', fontsize=14, y=1.02)
plt.tight_layout()
plt.savefig('bias_variance.png', dpi=150, bbox_inches='tight')
plt.show()

The visualization shows the tradeoff directly. At degree 1 (underfitting), the mean prediction (red) is far from the true function (black dashed) — high bias — but the individual fits (blue) cluster tightly — low variance. At degree 19 (overfitting), the mean prediction matches the true function well — low bias — but individual fits are wildly different — high variance. The sweet spot is somewhere in between.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

degrees = range(1, 20)
biases = []
variances = []

for degree in degrees:
    predictions = np.zeros((n_datasets, len(x_test)))
    for d in range(n_datasets):
        x_train = np.random.uniform(0, 1, n_train)
        y_train = f_true(x_train) + np.random.normal(0, 0.3, n_train)
        coeffs = np.polyfit(x_train, y_train, degree)
        predictions[d] = np.polyval(coeffs, x_test)
    mean_pred = predictions.mean(axis=0)
    biases.append(np.mean((mean_pred - f_true(x_test))**2))
    variances.append(np.mean(predictions.var(axis=0)))

mse = np.array(biases) + np.array(variances) + 0.3**2  # + noise variance

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(list(degrees), biases, 'b-o', label='Bias$^2$', markersize=5)
ax.plot(list(degrees), variances, 'r-o', label='Variance', markersize=5)
ax.plot(list(degrees), mse, 'k-o', label='Total MSE', markersize=5)
ax.axhline(y=0.3**2, color='gray', linestyle=':', label='Noise ($\\sigma^2$)')
ax.set_xlabel('Polynomial Degree', fontsize=13)
ax.set_ylabel('Error', fontsize=13)
ax.set_title('Bias-Variance Tradeoff', fontsize=14)
ax.legend(fontsize=12)
ax.set_yscale('log')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('bias_variance_curve.png', dpi=150)
plt.show()

Sufficient Statistics and Data Reduction#

The Factorization Theorem#

where $g$ depends on the data only through $T$ , and $h$ does not depend on $\theta$ .

So $T = \sum X_i$ is sufficient for $\lambda$ . Once you know the total count, the individual observations add no further information about $\lambda$ .

The Rao-Blackwell Theorem#

This theorem says: always condition on a sufficient statistic. It can only help, never hurt.

Cross-Validation: The Practical Bias-Variance Tool#

The bias-variance decomposition is a theoretical framework. In practice, you can’t compute bias or variance directly (they depend on the unknown true function). Cross-validation is the practical alternative.

k-fold cross-validation:

1. Split data into $k$ equal folds.
2. For each fold $i$ : train on all other folds, evaluate on fold $i$ .
3. Average the $k$ test errors.

This estimates the test error (which includes both bias and variance effects) without needing a separate test set. Common choices: $k = 5$ or $k = 10$ .

Leave-one-out cross-validation (LOOCV): $k = n$ . Each observation is held out once. Low bias (training set is almost the full dataset) but high variance (the $n$ training sets are very similar, so the $n$ error estimates are correlated).

The bias-variance tradeoff applies to cross-validation itself:

• Small $k$ : high bias (training sets are smaller), low variance
• Large $k$ : low bias, high variance (correlated estimates)
• $k = 5$ or $10$ is a good compromise in most settings

Exponential Families: A Unifying Framework#

where $T(x)$ is the sufficient statistic, $\eta(\theta)$ is the natural parameter, and $A(\theta)$ is the log-partition function.

For exponential families, the MLE is uniquely determined by the sufficient statistic, and it equals the method of moments estimator based on $T(X)$ . Conjugate priors also have a natural form, and the MLE always exists and is unique (under mild conditions).

Numerical MLE: When Closed Forms Don’t Exist#

Many models don’t have closed-form MLE solutions. In such cases, we optimize the log-likelihood numerically.

Newton-Raphson Method#

Fisher Scoring#

Fisher scoring is more stable than Newton-Raphson because the Fisher information matrix is guaranteed to be positive semi-definite (ensuring we move in an ascent direction). For exponential families, Newton-Raphson and Fisher scoring coincide.

The EM Algorithm#

For models with latent (unobserved) variables, direct MLE is often intractable. The Expectation-Maximization (EM) algorithm alternates between:

1. E-step: Compute $Q(\theta | \theta^{(t)}) = E_{\mathbf{Z}|\mathbf{X}, \theta^{(t)}}[\ell(\theta; \mathbf{X}, \mathbf{Z})]$ — the expected complete-data log-likelihood.
2. M-step: $\theta^{(t+1)} = \arg\max_\theta Q(\theta | \theta^{(t)})$ — maximize the expected log-likelihood.

Each iteration is guaranteed to increase (or maintain) the log-likelihood: $\ell(\theta^{(t+1)}) \geq \ell(\theta^{(t)})$ . The algorithm converges to a local maximum.

EM is used for Gaussian mixture models, hidden Markov models, factor analysis, and any model where marginalizing over latent variables makes the likelihood intractable but conditioning on them makes it tractable.

Gradient Descent for MLE#

where $\eta$ is the learning rate. This is gradient ascent on the log-likelihood — identical to what deep learning calls “training.” The MLE framework gives the theoretical justification: we are maximizing the likelihood of the observed data under the model.

In stochastic settings (mini-batches), we use stochastic gradient ascent, and the CLT from Article 5 guarantees that the noise in the gradient estimate is approximately Gaussian with variance $O(1/B)$ .

Summary#

What’s Next#

Estimation gives us point estimates — single best guesses for parameters. But how confident should we be? The next article tackles hypothesis testing and confidence intervals: the framework for quantifying uncertainty, controlling error rates, and making principled decisions from data. We’ll also see why p-values are so often misunderstood and how to avoid the most common statistical pitfalls.