Week 10: Bias, Variance, and Uncertainty

Bias and Variance#

Every estimator $\hat \theta$ is a random variable (it depends on the random sample). How do we judge whether an estimator is "good"? Two orthogonal qualities:

• Bias: $\text{Bias}(\hat \theta) = \mathbb{E}[\hat \theta] - \theta$. Does the estimator systematically over- or under-shoot the truth on average?
• Variance: $\mathrm{Var}(\hat \theta)$. How much does the estimate jump around from sample to sample?

Mean Squared Error (MSE) combines both:

Proof: Add and subtract $\mathbb{E}[\hat \theta]$:

The cross-term vanishes because $\mathbb{E}[\hat \theta - \mathbb{E}[\hat \theta]] = 0$.

This is the bias–variance decomposition: to reduce MSE, you can reduce bias or reduce variance. But estimators that reduce one often increase the other — the bias–variance tradeoff.

Estimator terminology:

• Unbiased: $\mathbb{E}[\hat \theta] = \theta$ (bias = 0).
• Consistent: $\hat \theta \xrightarrow{p} \theta$ as $n \to \infty$ (MSE → 0).
• Efficient: achieves the smallest possible variance among a class (e.g., among unbiased estimators).

Worked Estimators — Sample Variance#

The most famous bias–variance story. Given i.i.d. $X_i \sim \mathcal{N}(\mu, \sigma^2)$:

• MLE (denominator $n$): $\hat \sigma^2_{\text{MLE}} = \frac{1}{n}\sum (X_i - \bar X)^2$.
• Unbiased (denominator $n-1$): $S^2 = \frac{1}{n-1}\sum (X_i - \bar X)^2$.

$\mathbb{E}[\hat \sigma^2_{\text{MLE}}] = \frac{n-1}{n}\sigma^2 < \sigma^2$ — biased downward. $\mathbb{E}[S^2] = \sigma^2$ — unbiased.

Why $n-1$? The sample mean $\bar X$ is the value that minimises the sum of squared deviations. This "uses up" one degree of freedom, making the sum slightly too small on average. Dividing by $n-1$ corrects for this.

For large $n$, the difference is negligible. For small samples, the unbiased estimator has larger variance (its MSE can actually be higher than the biased MLE — the bias–variance tradeoff in action).

Confidence Interval (lite)#

A confidence interval gives a range of plausible values for $\theta$, calibrated by the sampling distribution of the estimator.

For the mean (CLT-based, large $n$): By the CLT, $\bar X_n \approx \mathcal{N}(\mu, \sigma^2 / n)$. Standardising:

Rearranging gives the $100(1-\alpha)\%$ confidence interval:

For a 95% CI, $z_{0.025} \approx 1.96$. In practice, $\sigma$ is unknown and replaced by the sample standard deviation $s$, which adds uncertainty; for large $n$ this is fine, for small $n$ use the $t$-distribution (not covered here).

Critical interpretation — what a CI is NOT:

A 95% CI does not mean "there is a 95% probability that $\mu$ lies in this interval." $\mu$ is a fixed constant (in the frequentist view). The correct interpretation: if we repeated the experiment many times and computed a 95% CI each time, approximately 95% of those intervals would contain the true $\mu$. The procedure works 95% of the time.

This is a subtle but important distinction from the Bayesian credible interval, where you can say "the probability that $\mu$ is in this interval is 95%," because $\mu$ is treated as a random variable with a posterior distribution.

Worked example. $n = 100$ observations, $\bar x = 52.3$, $s = 8.0$. A 95% CI for $\mu$:

Bootstrap#

The bootstrap is a computational method to estimate the sampling distribution of any statistic without assuming a parametric model. It replaces theory with computation.

Algorithm (nonparametric bootstrap):

1. From the original sample $x_1, \ldots, x_n$, draw a bootstrap sample $x_1^*, \ldots, x_n^*$ by sampling $n$ times with replacement from the original data.
2. Compute the statistic of interest $\hat \theta^* = g(x_1^*, \ldots, x_n^*)$ on the bootstrap sample.
3. Repeat steps 1–2 $B$ times (typically $B = 1000$ or more) to get $\hat \theta_1^*, \ldots, \hat \theta_B^*$.
4. The empirical distribution of the $\hat \theta_b^*$ approximates the sampling distribution of $\hat \theta$.

What the bootstrap gives you:

• Standard error: $\text{SE}_{\text{boot}} = \text{std}(\hat \theta_1^*, \ldots, \hat \theta_B^*)$.
• Confidence interval (percentile): the $(\alpha/2)$ and $(1-\alpha/2)$ quantiles of the bootstrap distribution.
• Bias estimate: $\mathbb{E}^*[\hat \theta^*] - \hat \theta$ (mean of bootstrap estimates minus original estimate).

The bootstrap works when: the empirical distribution of the sample is a good approximation of the true population distribution. It is especially useful when no closed-form formula for the standard error exists (e.g., median, correlation coefficient, or a complex model output).

When bootstrap fails: dependent data (need block bootstrap), heavy tails (need subsampling bootstrap), small $n$ (empirical distribution is too coarse).

What We Are Not Doing#

This course deliberately limits classical statistics to what supports ML practice:

• No Neyman–Pearson hypothesis testing (null vs alternative, rejection regions, p-values, Type I/II errors).
• No multiple testing corrections (Bonferroni, FDR).
• No ANOVA or experimental design.
• No power analysis.

These are important topics in classical applied statistics, but they are not prerequisites for understanding MLE, KL divergence, or probabilistic generative models — which is what this course targets.

Knowledge check#

Browser lab#

Bootstrap the sampling distribution of the sample mean and compare the bootstrap SE to the CLT formula.

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)

true_mu, true_sigma = 50.0, 10.0
n = 30
data = rng.normal(true_mu, true_sigma, size=n)

B = 2000
boot_means = np.empty(B)
for b in range(B):
    boot_sample = rng.choice(data, size=n, replace=True)
    boot_means[b] = boot_sample.mean()

se_boot = boot_means.std()
se_clt = data.std(ddof=1) / np.sqrt(n)

ci_lower = np.percentile(boot_means, 2.5)
ci_upper = np.percentile(boot_means, 97.5)

fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))

axes[0].hist(boot_means, bins=40, density=True, alpha=0.7, color="tab:blue")
axes[0].axvline(data.mean(), color="tab:red", linewidth=2, label=f"Sample mean = {data.mean():.2f}")
axes[0].axvline(ci_lower, color="black", linestyle="--", linewidth=1)
axes[0].axvline(ci_upper, color="black", linestyle="--", linewidth=1, label=f"95% CI: [{ci_lower:.2f}, {ci_upper:.2f}]")
axes[0].set_title("Bootstrap distribution of the mean")
axes[0].set_xlabel("Mean"); axes[0].legend()

x_grid = np.linspace(boot_means.min(), boot_means.max(), 200)
normal_approx = (1/(np.sqrt(2*np.pi)*se_clt)) * np.exp(-0.5*((x_grid - data.mean())/se_clt)**2)
axes[1].hist(boot_means, bins=40, density=True, alpha=0.5, label="Bootstrap")
axes[1].plot(x_grid, normal_approx, "tab:red", linewidth=2, label=f"CLT approx (SE={se_clt:.2f})")
axes[1].set_title("Bootstrap vs CLT normal approximation")
axes[1].set_xlabel("Mean"); axes[1].legend()

plt.tight_layout()
plt.show()

print(f"Sample: n={n}, mean={data.mean():.3f}, std={data.std(ddof=1):.3f}")
print(f"SE (CLT):  {se_clt:.4f}")
print(f"SE (boot): {se_boot:.4f}")
print(f"95% CI (boot): [{ci_lower:.3f}, {ci_upper:.3f}]")
print(f"95% CI (CLT):  [{data.mean()-1.96*se_clt:.3f}, {data.mean()+1.96*se_clt:.3f}]")