8. Statistical Estimation

Estimation is the step from data to a statement about the unknown numbers that describe a population: a mean, a proportion, a rate, a regression slope. This chapter covers how to judge an estimator (bias, variance, consistency, efficiency), the general methods that produce estimators (method of moments, maximum likelihood, least squares), confidence intervals, which turn a single number into a range with a known error rate, and resampling methods that estimate uncertainty when no formula is available. It builds directly on the sampling distributions of Chapter 7, and it supplies the machinery that hypothesis testing (Chapter 9) and regression (Chapter 10) rely on.

<a id="81-point-estimation"></a>

8.1 Point Estimation

A point estimator turns a sample into a single best guess for a parameter; this subchapter sets out the criteria (bias, variance, consistency, mean squared error and sufficiency) by which one estimator is judged better than another.

Estimators and Estimates

Why it matters. "The average order value is 47.30" is an answer, but it was produced by a rule applied to one particular sample. A different rule (the median, a trimmed mean) applied to the same data gives a different answer, and the same rule applied to a different sample gives a different answer too. To decide which rule to trust, we have to evaluate the rule rather than the single number it happened to produce. That separation, between the procedure and its output, is the starting point of estimation theory.

Vocabulary.

Term	Meaning	Example
Parameter (estimand) $\theta$	the unknown, fixed population quantity we want	the mean order value $\mu$ of all customers
Estimator $\hat\theta = T(X_1, \dots, X_n)$	a rule: a function of the random sample, hence a random variable	$\bar X = \frac{1}{n}\sum X_i$
Estimate $\hat\theta = T(x_1, \dots, x_n)$	the number the rule produces on the observed data	$\bar x = 47.30$
Sampling distribution	the distribution of the estimator over repeated samples (Section 7.3)	$\bar X \approx N(\mu, \sigma^2/n)$

The hat notation is used for both the estimator and the estimate; the context says which. Capital letters ($\bar X$) emphasize the random variable, lower case ($\bar x$) the observed value.

Analogy. An estimator is like an archer's technique, and an estimate is one arrow. A single arrow that lands near the bullseye might be luck; a coach judges the technique by where many arrows land: how far their center is from the bullseye (bias) and how widely they scatter (variance). The rest of this subchapter makes those two ideas precise.

  estimator = a rule                       sampling distribution of the rule
                                            (many samples, many estimates)
  sample 1 --> T(sample 1) = 49.4                    .
  sample 2 --> T(sample 2) = 51.1                   ...
  sample 3 --> T(sample 3) = 50.3                  .....      <- judged by its center
     ...                                          .......        and its spread
                                           -------------|-------
                                                    true theta

Any statistic can be an estimator. Nothing in the definition requires the estimator to be sensible: $T = 42$ (ignore the data) is a legitimate estimator of $\mu$, just a terrible one in most situations. Estimation theory does not tell you which rules are allowed; it tells you how to compare them.

The code applies four estimators of the center of a distribution to one sample, then to 20,000 samples from three populations with the same center.

import numpy as np

rng = np.random.default_rng(801)

# Four estimators of the center of a distribution: each is a RULE (a function of the sample)
estimators = {
    "mean":         np.mean,
    "median":       np.median,
    "midrange":     lambda x: (x.min() + x.max()) / 2,
    "10% trimmed":  lambda x: np.mean(np.sort(x)[len(x) // 10 : -(len(x) // 10)]),
}

# ONE sample gives one ESTIMATE per estimator (true center = 50)
sample = rng.normal(50, 10, size=30)
print("estimates from one sample: " + ", ".join(f"{k} = {f(sample):.2f}" for k, f in estimators.items()))

# MANY samples reveal each estimator's sampling distribution, which is what we judge them by
for name, draw in [("normal", lambda: rng.normal(50, 10, 30)),
                   ("uniform", lambda: rng.uniform(50 - 17.32, 50 + 17.32, 30)),   # same mean and SD
                   ("t, 2 df", lambda: 50 + 10 * rng.standard_t(2, 30))]:
    samples = [draw() for _ in range(20_000)]
    row = "  ".join(f"{k} {np.std([f(s) for s in samples]):5.2f}" for k, f in estimators.items())
    print(f"{name:8s} spread (SD) of the estimates:  {row}")

estimates from one sample: mean = 49.42, median = 48.51, midrange = 47.66, 10% trimmed = 49.02
normal   spread (SD) of the estimates:  mean  1.82  median  2.22  midrange  3.55  10% trimmed  1.87
uniform  spread (SD) of the estimates:  mean  1.82  median  3.00  midrange  0.78  10% trimmed  2.12
t, 2 df  spread (SD) of the estimates:  mean  7.18  median  2.60  midrange 87.70  10% trimmed  2.73

On one sample the four estimates differ, and nothing in that single sample says which is closest to the truth. Over many samples the picture is clear, and it depends on the population: for normal data the mean has the smallest spread, for uniform data the midrange wins by a wide margin (the extremes pin down the edges of a bounded distribution), and for heavy-tailed $t_2$ data the mean and midrange are badly unstable while the median and trimmed mean stay reliable. There is no estimator that is best everywhere; choosing one means making assumptions about the population.

Key takeaways

• A parameter is fixed and unknown; an estimator is a random variable (a rule applied to random data); an estimate is the number the rule produces on your data.
• Estimators are judged by their sampling distribution, never by a single estimate.
• The best estimator depends on the population model: the mean is ideal for normal data, poor for heavy tails.

🧪 Practice

1. Classify each as a parameter, an estimator or an estimate: (a) the proportion of all voters who support a law; (b) $\hat p = X/n$; (c) "41% of our 1,200 respondents support the law".
2. In the code, add the estimator $T = 50$ (a constant). What is its spread? Why is it still a bad estimator?
3. Explain why the midrange is so precise for uniform data and so poor for $t_2$ data.
4. Interview: Your team reports the mean session length as a KPI, but a few sessions last for days because of bots and open browser tabs. Which estimator would you propose instead, and how would you convince the team it is better? Hint: compare the estimators' sampling distributions by bootstrapping (Section 8.4) the real data; think about what parameter each one estimates.

Unbiasedness

Why it matters. An estimator that is systematically too high or too low misleads in the same direction every time, and collecting more samples of the same size does not average the error away. Bias is the first property we check: a scale that always reads 200 g too heavy is useless no matter how consistent it is. This topic also explains, as promised in Section 2.3, why the sample variance divides by $n - 1$.

Definition. The bias of an estimator $\hat\theta$ of $\theta$ is

$$\operatorname{Bias}(\hat\theta) = E[\hat\theta] - \theta,$$

and $\hat\theta$ is unbiased if its bias is zero for every possible value of $\theta$. The expectation is over the sampling distribution: "on average over all possible samples, the estimator hits the target". An estimator is asymptotically unbiased if its bias tends to 0 as $n \to \infty$.

Classic results (i.i.d. sample with mean $\mu$, variance $\sigma^2$):

Estimator	Target	Unbiased?
$\bar X$	$\mu$	yes: $E[\bar X] = \frac{1}{n}\sum E[X_i] = \mu$
$\hat p = X/n$	$p$	yes (a sample mean of 0/1 values)
$S^2 = \frac{1}{n-1}\sum (X_i - \bar X)^2$	$\sigma^2$	yes
$\hat\sigma^2 = \frac{1}{n}\sum (X_i - \bar X)^2$	$\sigma^2$	no: $E = \frac{n-1}{n}\sigma^2$
$S = \sqrt{S^2}$	$\sigma$	no: biased low
$X_{(n)} = \max X_i$, Uniform$(0, \theta)$	$\theta$	no: $E = \frac{n}{n+1}\theta$

Proof that $S^2$ is unbiased. Expand the sum of squares around $\mu$:

$$\sum_{i=1}^n (X_i - \bar X)^2 = \sum_{i=1}^n (X_i - \mu)^2 - n(\bar X - \mu)^2.$$

Take expectations. Each $E[(X_i - \mu)^2] = \sigma^2$, and $E[(\bar X - \mu)^2] = \operatorname{Var}(\bar X) = \sigma^2/n$. So

$$E\left[\sum (X_i - \bar X)^2\right] = n\sigma^2 - n \cdot \frac{\sigma^2}{n} = (n - 1)\sigma^2,$$

and dividing by $n - 1$ gives exactly $\sigma^2$. The lost $\sigma^2$ is the variance of $\bar X$ itself: measuring deviations from the estimated center rather than the true one makes them systematically smaller. The proof uses only independence and a finite variance, not normality.

Unbiasedness is not preserved by nonlinear transformations. If $E[\hat\theta] = \theta$, then $E[g(\hat\theta)] \ne g(\theta)$ in general. Since $\sqrt{\cdot}$ is concave, Jensen's inequality (Section 6.5) gives $E[S] \le \sqrt{E[S^2]} = \sigma$. The same happens with $1/\bar X$ as an estimator of a rate, or $\log \hat p$. This is one reason unbiasedness, while useful, is not the only criterion.

Removing bias. When the bias has a known form, rescale: $\frac{n+1}{n}X_{(n)}$ is unbiased for the uniform $\theta$, and $\frac{n}{n-1}\hat\sigma^2 = S^2$. When it does not, the jackknife (Section 8.4) and the bootstrap can estimate and subtract it.

import numpy as np

rng = np.random.default_rng(802)
sigma2, n, reps = 4.0, 5, 200_000
x = rng.normal(10, np.sqrt(sigma2), size=(reps, n))       # 200,000 samples of size 5

# Average each estimator over all samples: an unbiased estimator averages to the true value
s2_n1 = x.var(axis=1, ddof=1)        # divide by n - 1
s2_n = x.var(axis=1, ddof=0)         # divide by n
print(f"E[S^2] (n - 1) ~ {s2_n1.mean():.3f}   E[divide by n] ~ {s2_n.mean():.3f}   "
      f"theory: {sigma2}, {sigma2 * (n - 1) / n}")

# S^2 is unbiased for sigma^2, but S = sqrt(S^2) is biased low for sigma (Jensen's inequality)
s = np.sqrt(s2_n1)
print(f"E[S] ~ {s.mean():.3f}  vs sigma = {np.sqrt(sigma2):.3f}  (bias {s.mean() - np.sqrt(sigma2):+.3f})")

# Uniform(0, theta): the sample maximum is biased low; scaling by (n + 1)/n removes the bias
theta = 100
u = rng.uniform(0, theta, size=(reps, n))
mx = u.max(axis=1)
print(f"E[max] ~ {mx.mean():.2f}   E[(n+1)/n * max] ~ {((n + 1) / n * mx).mean():.2f}   "
      f"E[2 * mean] ~ {(2 * u.mean(axis=1)).mean():.2f}   (theta = {theta})")

E[S^2] (n - 1) ~ 3.997   E[divide by n] ~ 3.197   theory: 4.0, 3.2
E[S] ~ 1.879  vs sigma = 2.000  (bias -0.121)
E[max] ~ 83.36   E[(n+1)/n * max] ~ 100.03   E[2 * mean] ~ 100.02   (theta = 100)

With $n = 5$, dividing by $n$ underestimates $\sigma^2 = 4$ by 20% (the factor $(n-1)/n = 0.8$), while $S^2$ averages to the truth. Taking the square root reintroduces a bias of about 6% for $\sigma$. For the uniform distribution, the maximum is always below $\theta$, so it is biased low; both the rescaled maximum and $2\bar X$ are unbiased, and the next topic shows which of the two is better.

Unbiased is not the same as good. An unbiased estimator can have a huge variance (using only $X_1$ to estimate $\mu$ is unbiased), and for some problems no unbiased estimator exists at all, or the only unbiased one is absurd. Bias must be weighed together with variance, which leads to the mean squared error later in this subchapter.

Key takeaways

• Bias is $E[\hat\theta] - \theta$; unbiased means correct on average over repeated samples, for every value of $\theta$.
• $\bar X$, $\hat p$ and $S^2$ (divisor $n - 1$) are unbiased; the divisor-$n$ variance and $S$ are biased low.
• The $n - 1$ corrects for measuring deviations from $\bar X$ instead of $\mu$; the proof needs no normality.
• Nonlinear functions of unbiased estimators are generally biased (Jensen); unbiasedness alone does not make an estimator good.

🧪 Practice

1. Show that $\frac{1}{2}(X_1 + X_n)$ is an unbiased estimator of $\mu$ for any $n \ge 2$. Why would you still prefer $\bar X$?
2. For Uniform$(0, \theta)$, show that $E[X_{(n)}] = \frac{n}{n+1}\theta$ using the CDF $P(X_{(n)} \le x) = (x/\theta)^n$.
3. If $\hat p$ is unbiased for $p$, is $\hat p(1 - \hat p)$ unbiased for $p(1 - p)$? Compute $E[\hat p(1 - \hat p)]$ and find a correction factor.
4. Interview: Someone claims that the average of daily conversion rates over 30 days is an unbiased estimate of the monthly conversion rate. When is that true and when is it not? Hint: the monthly rate is total conversions over total visits; compare a simple average of ratios with a ratio of totals when daily traffic varies (Section 2.2, weighted means).

Efficiency and Minimum Variance

Why it matters. Many estimators are unbiased: for the uniform distribution both $2\bar X$ and $\frac{n+1}{n}X_{(n)}$ are. If both are right on average, the one whose estimates cluster more tightly around the truth is better, because any single estimate is then likely to be closer. Efficiency compares estimators by their variance, and a more efficient estimator is equivalent to free extra data: it reaches the same precision with fewer observations.

Relative efficiency. For two unbiased estimators $\hat\theta_1$ and $\hat\theta_2$ of $\theta$,

$$\operatorname{eff}(\hat\theta_1, \hat\theta_2) = \frac{\operatorname{Var}(\hat\theta_2)}{\operatorname{Var}(\hat\theta_1)}.$$

A value below 1 means $\hat\theta_1$ is less efficient. Since variances of typical estimators shrink like $1/n$, an efficiency of $0.64$ means $\hat\theta_1$ needs about $1/0.64 \approx 1.57$ times as many observations to match $\hat\theta_2$.

Minimum variance unbiased estimators. An unbiased estimator whose variance is the smallest among all unbiased estimators, for every value of $\theta$, is called the uniformly minimum variance unbiased estimator (UMVUE). Two tools find it:

• the Cramér-Rao lower bound (Section 8.2) gives a floor that no unbiased estimator can beat; an estimator that reaches it is efficient and therefore UMVUE;
• the Rao-Blackwell and Lehmann-Scheffé theorems (last topic of this subchapter) show that the UMVUE is a function of a sufficient statistic, and construct it.

Known comparisons.

Problem	Estimators	Efficiency of the first
Center of a normal	median vs. mean	$2/\pi \approx 0.64$ (large $n$)
Center of a Laplace (double exponential)	mean vs. median	$0.5$: the median is better
Uniform$(0, \theta)$	$2\bar X$ vs. $\frac{n+1}{n}X_{(n)}$	$\frac{3}{n+2}$: falls to 0 as $n$ grows
Normal mean, $n$ vs. $m$ observations	$\bar X_n$ vs. $\bar X_m$	$n/m$

The uniform case is striking. $\operatorname{Var}(2\bar X) = \theta^2/(3n)$ shrinks like $1/n$, but $\operatorname{Var}\left(\frac{n+1}{n}X_{(n)}\right) = \frac{\theta^2}{n(n+2)}$ shrinks like $1/n^2$. The maximum is a sufficient statistic for $\theta$, and the sample mean throws away most of the information.

import numpy as np

rng = np.random.default_rng(803)
theta, reps = 100, 50_000

# Two unbiased estimators of theta for Uniform(0, theta)
print(f"{'n':>4s} {'Var(2*mean)':>12s} {'Var(max adj.)':>14s} {'efficiency of 2*mean':>21s}")
for n in [5, 20, 100]:
    u = rng.uniform(0, theta, size=(reps, n))
    v_mom = (2 * u.mean(axis=1)).var()                   # theory theta^2 / (3n)
    v_max = ((n + 1) / n * u.max(axis=1)).var()          # theory theta^2 / (n (n + 2))
    print(f"{n:4d} {v_mom:12.2f} {v_max:14.2f} {v_max / v_mom:21.3f}")

# Mean vs median for the center of a normal distribution: efficiency -> 2/pi
n = 101
x = rng.normal(0, 1, size=(reps, n))
eff = x.mean(axis=1).var() / np.median(x, axis=1).var()
print(f"normal data, n = {n}: Var(mean) / Var(median) = {eff:.3f}   (2/pi = {2 / np.pi:.3f})")

   n  Var(2*mean)  Var(max adj.)  efficiency of 2*mean
   5       672.90         288.86                 0.429
  20       166.45          23.33                 0.140
 100        33.03           0.99                 0.030
normal data, n = 101: Var(mean) / Var(median) = 0.639   (2/pi = 0.637)

The simulated efficiencies match $3/(n+2)$: $3/7 = 0.43$, $3/22 = 0.14$ and $3/102 = 0.03$. With 100 observations, the estimator based on the mean is about 33 times noisier. For normal data, the median's efficiency relative to the mean is $2/\pi$: using the median wastes about 36% of the data. That is the price of its robustness, worth paying when outliers are likely (Section 2.2).

Efficiency is model-dependent. The mean is the most efficient estimator of the center for normal data, but the comparison flips for heavy-tailed data (the $t_2$ row of the first topic). A robust estimator such as the trimmed mean or the Huber estimator gives up a little efficiency at the normal model in exchange for protection elsewhere.

Key takeaways

• Among unbiased estimators, prefer the one with smaller variance; relative efficiency is the variance ratio.
• An efficiency of $e$ means the less efficient estimator needs about $1/e$ times the sample size to match.
• The UMVUE has the smallest variance among all unbiased estimators; Cramér-Rao and Rao-Blackwell are the tools to find it.
• Efficiency depends on the model: the mean beats the median for normal data and loses for heavy tails.

🧪 Practice

1. Using $\operatorname{Var}(X_{(n)}) = \frac{n\theta^2}{(n+1)^2(n+2)}$, derive $\operatorname{Var}\left(\frac{n+1}{n}X_{(n)}\right)$ and the efficiency $3/(n+2)$.
2. Survey A uses the median of 400 normal responses, survey B the mean. How many responses would B need to match A's precision?
3. Estimators $\hat\theta_1$ and $\hat\theta_2$ are unbiased and independent with variances $\sigma_1^2$ and $\sigma_2^2$. Find the weight $w$ that minimizes the variance of $w\hat\theta_1 + (1 - w)\hat\theta_2$ (inverse-variance weighting).
4. Interview: Two teams measure the same latency metric; one reports the mean of 1,000 requests, the other the median. Their results differ. How do you decide which estimate to rely on? Hint: first ask whether they estimate the same parameter (mean and median differ for skewed data); then compare their standard errors, for example by bootstrapping.

Consistency

Why it matters. Unbiasedness and variance describe an estimator at a fixed sample size. The most basic requirement, though, concerns what happens as data accumulate: with enough data, the estimator should home in on the truth. An estimator that stays wrong no matter how many observations you collect is broken in a way that no amount of data can fix. Consistency is that minimal guarantee.

Definition. A sequence of estimators $\hat\theta_n$ is consistent for $\theta$ if it converges in probability to $\theta$ (Section 6.5):

$$\lim_{n \to \infty} P\big(|\hat\theta_n - \theta| > \varepsilon\big) = 0 \quad \text{for every } \varepsilon > 0.$$

In words: for any tolerance, the probability that the estimate misses by more than that tolerance goes to zero as $n$ grows.

A practical sufficient condition. If both the bias and the variance tend to zero, the estimator is consistent. This follows from Markov's inequality applied to the squared error:

$$P\big(|\hat\theta_n - \theta| > \varepsilon\big) \le \frac{E\big[(\hat\theta_n - \theta)^2\big]}{\varepsilon^2} = \frac{\operatorname{Var}(\hat\theta_n) + \operatorname{Bias}(\hat\theta_n)^2}{\varepsilon^2} \to 0.$$

Tools for proving consistency.

• Law of large numbers: sample means of i.i.d. variables with a finite mean are consistent for the population mean.
• Continuous mapping: if $\hat\theta_n \to \theta$ in probability and $g$ is continuous at $\theta$, then $g(\hat\theta_n) \to g(\theta)$. So $S = \sqrt{S^2}$ is consistent for $\sigma$ even though it is biased.
• Combining: sums, products and ratios of consistent estimators are consistent for the corresponding combinations (Slutsky's theorem).

Consistency vs. unbiasedness. The two properties are independent of each other:

	Unbiased	Biased
Consistent	$\bar X$ for $\mu$	$\hat\sigma^2$ (divisor $n$); $S$ for $\sigma$
Inconsistent	$X_1$ alone for $\mu$	the constant $42$

import numpy as np

rng = np.random.default_rng(804)
eps, reps = 0.1, 4_000

# P(|estimate - truth| > eps) for growing n: a consistent estimator drives this to 0
estimators = {
    "mean (normal)":        (lambda n: rng.normal(0, 1, (reps, n)).mean(axis=1), 0.0),
    "var ddof=0 (normal)":  (lambda n: rng.normal(0, 1, (reps, n)).var(axis=1), 1.0),
    "first obs. only":      (lambda n: rng.normal(0, 1, (reps, n))[:, 0], 0.0),
    "mean (Cauchy)":        (lambda n: rng.standard_cauchy((reps, n)).mean(axis=1), 0.0),
    "median (Cauchy)":      (lambda n: np.median(rng.standard_cauchy((reps, n)), axis=1), 0.0),
}
ns = [10, 100, 1_000, 5_000]
print(f"{'estimator':22s}" + "".join(f"{'n=' + str(n):>9s}" for n in ns))
for name, (est, truth) in estimators.items():
    probs = [np.mean(np.abs(est(n) - truth) > eps) for n in ns]
    print(f"{name:22s}" + "".join(f"{p:9.3f}" for p in probs))

estimator                  n=10    n=100   n=1000   n=5000
mean (normal)             0.773    0.320    0.002    0.000
var ddof=0 (normal)       0.847    0.486    0.024    0.000
first obs. only           0.913    0.927    0.926    0.917
mean (Cauchy)             0.933    0.940    0.942    0.933
median (Cauchy)           0.829    0.526    0.047    0.000

For the consistent estimators, the probability of missing by more than 0.1 falls toward 0. The divisor-$n$ variance is biased, but its bias $-\sigma^2/n$ vanishes, so it is consistent too. Using only the first observation never improves, because it ignores the extra data. The Cauchy distribution is the instructive failure: it has no mean, the law of large numbers does not apply, and the sample mean of $n$ Cauchy variables has the same Cauchy distribution as a single observation. The median of the same data is consistent for the center.

Key takeaways

• Consistent means $\hat\theta_n \to \theta$ in probability: more data eventually gives arbitrarily accurate estimates.
• Bias and variance both tending to 0 is enough (via the MSE and Markov's inequality).
• Consistency and unbiasedness are different properties; a biased estimator can be consistent and an unbiased one inconsistent.
• Continuous functions of consistent estimators are consistent, which makes plug-in estimators work.

🧪 Practice

1. Show that $\hat p = X/n$ is consistent for $p$ using Chebyshev's inequality, and find an $n$ that guarantees $P(|\hat p - p| > 0.05) \le 0.05$ for every $p$.
2. Is $\frac{n+1}{n}X_{(n)}$ consistent for the Uniform$(0, \theta)$ endpoint? Is $X_{(n)}$ itself consistent?
3. An estimator equals $\bar X$ with probability $1 - 1/n$ and $n^2$ with probability $1/n$. Show it is consistent for $\mu$ although its bias does not go to 0.
4. Interview: A recommender's click-through rate is estimated from logs that only include users who saw at least one recommendation. More data keeps arriving. Will the estimate converge to the CTR of all users? Hint: consistency is about converging to the right target; a selection bias (Section 7.4) makes the estimator consistent for the wrong parameter.

Mean Squared Error and Bias-Variance Decomposition

Why it matters. Insisting on unbiasedness can be costly. An estimator with a small bias and much less variance often lands closer to the truth than an unbiased but noisy one. To compare estimators that differ in both respects, we need one number that combines them. The mean squared error is that number, and its decomposition into bias and variance is the same trade-off that drives regularization and model selection in machine learning (Sections 10.4 and 12.5).

Definition.

$$\operatorname{MSE}(\hat\theta) = E\big[(\hat\theta - \theta)^2\big] = \operatorname{Var}(\hat\theta) + \operatorname{Bias}(\hat\theta)^2.$$

Proof. Add and subtract $E[\hat\theta]$:

$$E\big[(\hat\theta - \theta)^2\big] = E\Big[\big(\hat\theta - E[\hat\theta]\big)^2\Big] + 2,\big(E[\hat\theta] - \theta\big),\underbrace{E\big[\hat\theta - E[\hat\theta]\big]}_{=,0} + \big(E[\hat\theta] - \theta\big)^2.$$

The cross term vanishes, leaving variance plus squared bias. For an unbiased estimator, MSE equals variance. The root mean squared error $\sqrt{\operatorname{MSE}}$ is in the units of $\theta$ and is the typical size of the estimation error.

   low bias, high variance      high bias, low variance      low bias, low variance

          .      .
       .     .      .                  ...
         .  (+)  .              (+)   .....                       .:.
       .   .   .     .                 ...                       .(+).
          .     .                                                 ':'

   right on average, but       tightly grouped, but off     the goal: every estimate
   estimates scatter widely    center in one direction      lands close to (+)

Example: estimating a normal variance. For normal data, among estimators of the form $\sum (X_i - \bar X)^2 / d$, which divisor $d$ minimizes the MSE? Using $\operatorname{Var}(S^2) = 2\sigma^4/(n-1)$ (Section 7.3),

Divisor $d$	Bias	Variance	MSE
$n - 1$	$0$	$\frac{2\sigma^4}{n-1}$	$\frac{2\sigma^4}{n-1}$
$n$	$-\frac{\sigma^2}{n}$	$\frac{2(n-1)\sigma^4}{n^2}$	$\frac{(2n-1)\sigma^4}{n^2}$
$n + 1$	$-\frac{2\sigma^2}{n+1}$	$\frac{2(n-1)\sigma^4}{(n+1)^2}$	$\frac{2\sigma^4}{n+1}$ (smallest)

The unbiased choice is not the one with the smallest MSE: dividing by $n + 1$ shrinks the estimate toward zero, and the reduction in variance more than pays for the bias.

Shrinkage. The same idea in its simplest form: estimate $\mu$ by $c\bar X$ with $0 < c \le 1$. Its MSE is $c^2\sigma^2/n + (1 - c)^2\mu^2$, which is minimized at

$$c^* = \frac{\mu^2}{\mu^2 + \sigma^2/n}.$$

The optimal amount of shrinkage depends on the unknown $\mu$. When the signal is small relative to the noise, shrinking helps a lot; when it is large, shrinking hurts. Ridge regression, the James-Stein estimator and Bayesian priors (Chapter 11) are practical ways to choose the amount of shrinkage from the data.

import numpy as np

rng = np.random.default_rng(805)
sigma2, n, reps = 4.0, 10, 400_000
x = rng.normal(0, np.sqrt(sigma2), size=(reps, n))
ss = ((x - x.mean(axis=1, keepdims=True)) ** 2).sum(axis=1)     # sum of squared deviations

# Estimate sigma^2 by SS / d for different divisors d: MSE = Var + Bias^2
print(f"{'divisor':>8s} {'bias':>8s} {'variance':>9s} {'bias^2 + var':>13s} {'MSE (direct)':>13s}")
for d in [n - 1, n, n + 1]:
    est = ss / d
    bias, var = est.mean() - sigma2, est.var()
    print(f"{d:8d} {bias:8.3f} {var:9.3f} {bias**2 + var:13.3f} {np.mean((est - sigma2) ** 2):13.3f}")

# Shrinking the sample mean toward 0: c * xbar trades bias for variance (data SD 2, n = 10)
print("MSE of c * xbar:")
for mu in [0.5, 3.0]:
    xbar = rng.normal(mu, 2, size=(reps, n)).mean(axis=1)
    row = "  ".join(f"c={c:.1f}: {np.mean((c * xbar - mu) ** 2):.3f}" for c in [1.0, 0.9, 0.8, 0.5])
    print(f"  true mu = {mu}:  {row}")

 divisor     bias  variance  bias^2 + var  MSE (direct)
       9    0.002     3.562         3.562         3.562
      10   -0.398     2.885         3.044         3.044
      11   -0.726     2.384         2.911         2.911
MSE of c * xbar:
  true mu = 0.5:  c=1.0: 0.402  c=0.9: 0.328  c=0.8: 0.267  c=0.5: 0.163
  true mu = 3.0:  c=1.0: 0.400  c=0.9: 0.413  c=0.8: 0.615  c=0.5: 2.349

The simulated MSEs match the table: with $n = 10$ and $\sigma^2 = 4$, theory gives $32/9 = 3.56$, $76/25 = 3.04$ and $32/11 = 2.91$. For the shrinkage estimator, $c = 0.5$ cuts the MSE by 60% when $\mu = 0.5$, which is small compared to the standard error $2/\sqrt{10} \approx 0.63$. When $\mu = 3$, the same shrinkage multiplies the MSE by almost 6. Shrinkage is a bet on where the parameter lies.

Key takeaways

• $\operatorname{MSE} = \operatorname{Var} + \operatorname{Bias}^2$: one number that rewards both accuracy and precision.
• A slightly biased estimator can have a much smaller MSE than an unbiased one (for normal variance, divisor $n + 1$ beats $n - 1$).
• Shrinking toward a fixed value reduces variance at the cost of bias; whether it pays off depends on the unknown parameter.
• The same decomposition explains the bias-variance trade-off in prediction and the purpose of regularization.

🧪 Practice

1. Verify the MSE formula for divisor $n$ in the table, starting from $\operatorname{Bias} = -\sigma^2/n$ and $\operatorname{Var}(\hat\sigma^2) = \left(\frac{n-1}{n}\right)^2 \frac{2\sigma^4}{n-1}$.
2. Derive $c^*$ for the shrinkage estimator by differentiating its MSE with respect to $c$.
3. For Uniform$(0, \theta)$, find the constant $c$ that minimizes the MSE of $cX_{(n)}$. Is the resulting estimator unbiased?
4. Interview: A product manager says, "We should only use unbiased estimators; biased ones are wrong by definition." How would you respond? Hint: give the variance divisor or ridge regression as an example where accepting a known small bias reduces the typical error; also mention that consistency matters more for large samples.

Sufficiency and the Factorization Theorem

Why it matters. A dataset of a million coin flips takes a lot of storage, but for estimating the probability of heads you only need one number: how many heads there were. The order of the flips adds nothing. A sufficient statistic is a summary that keeps all the information the data contain about a parameter, so nothing is lost by discarding the rest. Sufficiency explains why good estimators look the way they do (they depend on the data only through sufficient statistics), it lets you improve any estimator (Rao-Blackwell), and it underlies data compression for streaming and distributed computation.

Definition. A statistic $T(X_1, \dots, X_n)$ is sufficient for $\theta$ if the conditional distribution of the sample given $T$ does not depend on $\theta$:

$$P(X_1 = x_1, \dots, X_n = x_n \mid T = t) \text{ is the same for every } \theta.$$

Interpretation: once you know $T$, the remaining detail of the data is pure noise as far as $\theta$ is concerned. You could regenerate a dataset with the same distribution from $T$ alone, without knowing $\theta$.

Example. For $n$ Bernoulli$(p)$ trials with $T = \sum X_i = t$, each specific sequence has probability $p^t(1 - p)^{n - t}$, and there are $\binom{n}{t}$ such sequences, so

$$P(\text{sequence} \mid T = t) = \frac{p^t (1 - p)^{n - t}}{\binom{n}{t} p^t (1 - p)^{n - t}} = \frac{1}{\binom{n}{t}},$$

which does not involve $p$. The number of successes is sufficient.

Factorization theorem (Fisher-Neyman). Checking the definition directly is tedious. The practical criterion: $T$ is sufficient for $\theta$ if and only if the joint PMF or PDF factors as

$$f(x_1, \dots, x_n; \theta) = g\big(T(x); \theta\big) \cdot h(x),$$

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

Model	Joint density (up to factors without $\theta$)	Sufficient statistic
Bernoulli$(p)$	$p^{\sum x_i}(1 - p)^{n - \sum x_i}$	$\sum X_i$
Poisson$(\lambda)$	$e^{-n\lambda}\lambda^{\sum x_i} \cdot \frac{1}{\prod x_i!}$	$\sum X_i$
Exponential$(\lambda)$	$\lambda^n e^{-\lambda \sum x_i}$	$\sum X_i$
Normal$(\mu, \sigma^2)$	$(2\pi\sigma^2)^{-n/2}\exp\left(-\frac{\sum x_i^2 - 2\mu\sum x_i + n\mu^2}{2\sigma^2}\right)$	$\left(\sum X_i, \sum X_i^2\right)$
Uniform$(0, \theta)$	$\theta^{-n},\mathbf{1}{\max x_i \le \theta}$	$X_{(n)} = \max X_i$

The normal row explains why a mean and a variance (or equivalently a running sum and sum of squares) summarize normal data completely. The uniform row explains the efficiency result of the previous topics: the maximum contains all the information about $\theta$, the mean does not.

Minimal sufficiency. The full dataset is always sufficient, trivially. The useful summary is the minimal sufficient statistic: the coarsest one, a function of every other sufficient statistic. The statistics in the table are minimal. Models in the exponential family (Section 10.5) always have a fixed-dimension sufficient statistic, however large $n$ is; most other models, such as the Cauchy, do not.

Rao-Blackwell theorem. If $\tilde\theta$ is any unbiased estimator and $T$ is sufficient, then

$$\hat\theta = E[\tilde\theta \mid T]$$

is also unbiased (law of iterated expectation) and $\operatorname{Var}(\hat\theta) \le \operatorname{Var}(\tilde\theta)$ (law of total variance, Section 6.2). Sufficiency guarantees that $\hat\theta$ does not depend on the unknown $\theta$, so it is a valid estimator. Conditioning on a sufficient statistic averages out noise that carries no information. If $T$ is also complete (a technical condition that holds for the exponential family), the Lehmann-Scheffé theorem states that the result is the unique UMVUE.

import numpy as np
from collections import Counter

rng = np.random.default_rng(806)

# Sufficiency: given the number of successes T = 3 in n = 5 Bernoulli trials,
# every arrangement is equally likely, whatever p is. The pattern carries no extra information about p.
for p in [0.2, 0.7]:
    x = rng.random((400_000, 5)) < p
    patterns = Counter("".join("1" if v else "0" for v in row) for row in x[x.sum(axis=1) == 3])
    total = sum(patterns.values())
    freqs = sorted(round(c / total, 3) for c in patterns.values())
    print(f"p = {p}: {len(patterns)} patterns with T = 3, conditional frequencies {freqs[0]} to {freqs[-1]} (1/10 each)")

# Rao-Blackwell: estimate P(X = 0) = exp(-lam) for Poisson data
lam, n, reps = 2.0, 20, 200_000
x = rng.poisson(lam, size=(reps, n))
naive = (x[:, 0] == 0).astype(float)            # uses only the first observation: unbiased but crude
t = x.sum(axis=1)                               # sufficient statistic
rb = ((n - 1) / n) ** t                         # E[naive | T] = ((n - 1)/n)^T: still unbiased
plug_in = np.exp(-x.mean(axis=1))               # MLE plug-in, for comparison: slightly biased
print(f"true exp(-lam) = {np.exp(-lam):.4f}")
for name, est in [("naive 1{X1 = 0}", naive), ("Rao-Blackwellized", rb), ("plug-in exp(-xbar)", plug_in)]:
    print(f"{name:20s} mean {est.mean():.4f}   SD {est.std():.4f}")

p = 0.2: 10 patterns with T = 3, conditional frequencies 0.094 to 0.104 (1/10 each)
p = 0.7: 10 patterns with T = 3, conditional frequencies 0.099 to 0.101 (1/10 each)
true exp(-lam) = 0.1353
naive 1{X1 = 0}      mean 0.1344   SD 0.3411
Rao-Blackwellized    mean 0.1354   SD 0.0438
plug-in exp(-xbar)   mean 0.1422   SD 0.0449

Whether $p = 0.2$ or $p = 0.7$, the 10 possible arrangements of 3 successes in 5 trials appear equally often: given $T$, the arrangement carries no information about $p$. For the Poisson problem, the naive estimator only records whether the first observation was zero; it is unbiased but extremely noisy (SD 0.34 for a target of 0.135). Conditioning it on the sufficient statistic $\sum X_i$ gives $\left(\frac{n-1}{n}\right)^{\sum X_i}$: still unbiased, with a standard deviation almost eight times smaller. The plug-in estimator $e^{-\bar X}$ has a similar spread but a small bias.

Key takeaways

• A sufficient statistic captures all the information in the sample about $\theta$; given it, the data's remaining detail is independent of $\theta$.
• Use the factorization $f(x; \theta) = g(T(x); \theta),h(x)$ to find sufficient statistics.
• Good estimators depend on the data only through a (minimal) sufficient statistic.
• Rao-Blackwell: conditioning an unbiased estimator on a sufficient statistic never increases the variance; with completeness it gives the UMVUE.

🧪 Practice

1. Use the factorization theorem to find a sufficient statistic for $\lambda$ in an i.i.d. Exponential$(\lambda)$ sample, and for $\theta$ in a Uniform$(\theta, \theta + 1)$ sample (two-dimensional).
2. Show that, given $T = \sum X_i = t$ for an i.i.d. Poisson sample, $(X_1, \dots, X_n)$ is Multinomial$(t; 1/n, \dots, 1/n)$. Use it to derive $E[\mathbf{1}{X_1 = 0} \mid T = t] = \left(\frac{n-1}{n}\right)^t$.
3. For a normal sample, explain why storing $n$, $\sum x_i$ and $\sum x_i^2$ is enough to compute the MLE of $(\mu, \sigma^2)$ for any amount of data. What numerical problem can the formula $\sum x_i^2 - (\sum x_i)^2/n$ cause (Section 2.3)?
4. Interview: You need the mean and variance of a metric across 10,000 servers without moving the raw data. What does each server send, and how are the results combined? Hint: each server sends a sufficient summary $(n_k, \text{mean}_k, M2_k)$; combine pairs with the parallel form of Welford's algorithm.

<a id="82-estimation-methods"></a>

8.2 Estimation Methods

The criteria of Section 8.1 say how to judge an estimator but not how to find one; this subchapter covers the general recipes that produce estimators for almost any model (the method of moments, maximum likelihood and least squares) and the Fisher information, which sets a limit on how precise any estimator can be.

Method of Moments

Why it matters. Suppose you model repair times with a gamma distribution and need its two parameters. There is no obvious "sample version" of a shape parameter. But the model says what the mean and variance should be as functions of the parameters, and the data tell you what the mean and variance actually are. Setting the two equal and solving gives estimates. The method of moments, introduced by Karl Pearson in the 1890s, is the oldest general estimation recipe. It is quick, often has a closed form, and provides good starting values for more refined methods.

The recipe. For a model with $k$ unknown parameters $\theta = (\theta_1, \dots, \theta_k)$:

1. Write the first $k$ population moments as functions of the parameters: $\mu_j(\theta) = E[X^j]$ (or central moments for $j \ge 2$).
2. Compute the corresponding sample moments: $m_j = \frac{1}{n}\sum_{i=1}^n X_i^j$.
3. Solve the $k$ equations $\mu_j(\hat\theta) = m_j$ for $\hat\theta$.

The justification is the law of large numbers: sample moments converge to population moments, so if the map from parameters to moments is continuous and invertible, the solution converges to the true parameters. Method-of-moments estimators are therefore consistent under mild conditions.

Examples.

Model	Moment equations	Estimators
Poisson$(\lambda)$	$\lambda = \bar X$	$\hat\lambda = \bar X$
Exponential$(\lambda)$	$1/\lambda = \bar X$	$\hat\lambda = 1/\bar X$
Normal$(\mu, \sigma^2)$	$\mu = \bar X$, $\sigma^2 = \hat\sigma^2$	$\bar X$ and the divisor-$n$ variance
Gamma(shape $k$, scale $s$)	$ks = \bar X$, $ks^2 = \hat\sigma^2$	$\hat k = \bar X^2/\hat\sigma^2$, $\hat s = \hat\sigma^2/\bar X$
Uniform$(0, \theta)$	$\theta/2 = \bar X$	$\hat\theta = 2\bar X$
Binomial$(N, p)$, both unknown	$Np = \bar X$, $Np(1-p) = \hat\sigma^2$	$\hat p = 1 - \hat\sigma^2/\bar X$, $\hat N = \bar X/\hat p$

Here $\hat\sigma^2 = \frac{1}{n}\sum (X_i - \bar X)^2$ is the second sample central moment.

import numpy as np

rng = np.random.default_rng(807)

# Gamma(shape k, scale s): mean = k s, variance = k s^2. Match them to the sample moments.
k_true, s_true = 2.5, 4.0
x = rng.gamma(k_true, s_true, size=200)              # e.g. repair times in hours
m, v = x.mean(), x.var()                             # first moment and second central moment (divisor n)
k_hat, s_hat = m**2 / v, v / m                       # solve  k s = m,  k s^2 = v
print(f"sample mean {m:.2f}, variance {v:.2f}  ->  shape {k_hat:.2f} (true {k_true}), scale {s_hat:.2f} (true {s_true})")

# Sampling distribution of the method-of-moments shape estimator
xs = rng.gamma(k_true, s_true, size=(10_000, 200))
k_hats = xs.mean(axis=1) ** 2 / xs.var(axis=1)
print(f"over 10,000 samples: mean of shape estimates {k_hats.mean():.3f}, SD {k_hats.std():.3f}")

# Uniform(0, theta): MoM gives 2 * mean, which can be IMPOSSIBLE given the data
u = rng.uniform(0, 10, size=(10_000, 10))
mom = 2 * u.mean(axis=1)
print(f"uniform, n = 10: 2*mean < sample max in {np.mean(mom < u.max(axis=1)):.1%} of samples")

sample mean 9.90, variance 32.46  ->  shape 3.02 (true 2.5), scale 3.28 (true 4.0)
over 10,000 samples: mean of shape estimates 2.552, SD 0.298
uniform, n = 10: 2*mean < sample max in 27.8% of samples

The gamma estimates are in the right region, and their average over 10,000 samples is close to the true shape 2.5. The uniform example shows the method's main weakness: in 28% of samples of size 10, the estimate $2\bar x$ is smaller than the largest observation, which is logically impossible because every observation is at most $\theta$. The method matches averages and ignores everything else the model says about the data.

Strengths and weaknesses.

Strengths	Weaknesses
simple, often closed-form	usually less efficient than maximum likelihood
consistent under mild conditions	can give impossible values (outside the parameter space or contradicting the data)
needs only moments, not the full density	depends on which moments you choose to match
good starting values for numerical MLE	higher sample moments are very sensitive to outliers

The generalized method of moments (GMM), widely used in econometrics, extends the idea to more moment conditions than parameters, combining them with optimal weights.

Key takeaways

• Equate the first $k$ population moments to the sample moments and solve for the $k$ parameters.
• The estimators are consistent by the law of large numbers, and often available in closed form.
• They are generally less efficient than maximum likelihood and can produce values that contradict the data.
• They make good starting values for numerical optimization of the likelihood.

🧪 Practice

1. Find the method-of-moments estimator for $p$ in a Geometric$(p)$ model (number of trials until the first success, mean $1/p$).
2. For Beta$(a, b)$ with mean $\mu = \frac{a}{a+b}$ and variance $\frac{\mu(1-\mu)}{a+b+1}$, derive the method-of-moments estimators of $a$ and $b$.
3. The counts $3, 5, 2, 6, 4, 4, 7, 3$ come from a Binomial$(N, p)$ with both unknown. Compute $\hat N$ and $\hat p$. What happens if the sample variance exceeds the sample mean?
4. Interview: You have click counts per user and suspect overdispersion relative to a Poisson model. How could you fit a negative binomial model quickly without an optimizer? Hint: match the mean $\mu$ and the variance $\mu + \mu^2/r$; solve for $r$, and check whether the estimate makes sense (it is negative if the variance is below the mean).

Maximum Likelihood Estimation

Why it matters. Maximum likelihood is the default estimation method of modern statistics: logistic regression, generalized linear models, most of machine learning's probabilistic models (where minimizing cross-entropy is maximizing a likelihood) and the fitting functions of every statistics package use it. The idea is intuitive, it works for almost any model you can write a density for, and in large samples it produces estimators that are as precise as possible (next topics).

The idea. Given observed data, consider each candidate value of the parameter and ask: if this were the true value, how probable would the data I actually observed be? Pick the value under which the observed data are most probable.

Analogy. A detective with several suspects asks, for each one, "if this person did it, how likely is the evidence we found?" and focuses on the suspect who best explains the evidence. The detective is not computing the probability that each suspect is guilty (that would need prior information, the Bayesian approach of Chapter 11), only how well each explains the facts.

Definition. For i.i.d. data $x_1, \dots, x_n$ from a density or PMF $f(x; \theta)$, the likelihood function is

$$L(\theta) = \prod_{i=1}^n f(x_i; \theta),$$

and the maximum likelihood estimate (MLE) is

$$\hat\theta_{\text{MLE}} = \arg\max_\theta L(\theta) = \arg\max_\theta \ell(\theta), \qquad \ell(\theta) = \log L(\theta) = \sum_{i=1}^n \log f(x_i; \theta).$$

The log does not change the maximizer (it is increasing) and turns the product into a sum, which is easier to differentiate and numerically stable (next topic).

The standard procedure.

1. Write the log-likelihood $\ell(\theta)$.
2. Differentiate and set the score to zero: $\ell'(\theta) = 0$ (the likelihood equation).
3. Solve for $\theta$; if there is no closed form, solve numerically (Newton-Raphson, gradient methods).
4. Check that it is a maximum ($\ell''(\hat\theta) < 0$) and look at the boundaries of the parameter space.

Worked examples.

Bernoulli. With $t = \sum x_i$ successes: $\ell(p) = t\log p + (n - t)\log(1 - p)$. Then $\ell'(p) = \frac{t}{p} - \frac{n - t}{1 - p} = 0$ gives $\hat p = t/n = \bar x$.

Exponential. $\ell(\lambda) = n\log\lambda - \lambda\sum x_i$, so $\ell'(\lambda) = \frac{n}{\lambda} - \sum x_i = 0$ gives $\hat\lambda = 1/\bar x$.

Normal. $\ell(\mu, \sigma^2) = -\frac{n}{2}\log(2\pi\sigma^2) - \frac{1}{2\sigma^2}\sum (x_i - \mu)^2$. The partial derivatives give $\hat\mu = \bar x$ and $\hat\sigma^2 = \frac{1}{n}\sum (x_i - \bar x)^2$, the divisor-$n$ variance, which is biased (Section 8.1).

Uniform$(0, \theta)$. $L(\theta) = \theta^{-n}$ if $\theta \ge \max x_i$ and $0$ otherwise, because a value of $\theta$ below the largest observation makes that observation impossible. $\theta^{-n}$ decreases in $\theta$, so the maximum is at the boundary: $\hat\theta = \max x_i$. Calculus would not find it; the derivative is never zero.

import numpy as np

rng = np.random.default_rng(808)

# Bernoulli: 7 conversions in 20 visits. Evaluate the likelihood on a grid of p values.
x = np.array([1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0])
p_grid = np.linspace(0.001, 0.999, 999)
lik = p_grid ** x.sum() * (1 - p_grid) ** (len(x) - x.sum())
print(f"grid maximizer p = {p_grid[lik.argmax()]:.3f}   closed form xbar = {x.mean():.3f}")

# Exponential: lambda_hat = 1 / xbar; check numerically on the log-likelihood
t = rng.exponential(scale=1 / 0.4, size=50)              # waiting times, true rate 0.4 per minute
lam_grid = np.linspace(0.05, 1.5, 2_901)
loglik = len(t) * np.log(lam_grid) - lam_grid * t.sum()   # n log(lambda) - lambda * sum(t)
print(f"grid maximizer lambda = {lam_grid[loglik.argmax()]:.3f}   closed form 1/tbar = {1 / t.mean():.3f}")

# Normal: both parameters, closed form
y = rng.normal(170, 8, size=30)
mu_hat, sigma2_hat = y.mean(), np.mean((y - y.mean()) ** 2)    # MLE of the variance divides by n
print(f"normal MLE: mu = {mu_hat:.2f}, sigma^2 = {sigma2_hat:.2f} (S^2 with n - 1 would give {y.var(ddof=1):.2f})")

# Uniform(0, theta): L(theta) = theta^-n for theta >= max(x), else 0 -> maximized at the boundary
u = rng.uniform(0, 10, size=8)
for theta in [u.max() - 0.01, u.max(), u.max() + 0.5, 12]:
    L = theta ** -len(u) if theta >= u.max() else 0.0
    print(f"  theta = {theta:6.3f}: L = {L:.3e}")

grid maximizer p = 0.350   closed form xbar = 0.350
grid maximizer lambda = 0.415   closed form 1/tbar = 0.416
normal MLE: mu = 168.68, sigma^2 = 49.70 (S^2 with n - 1 would give 51.42)
  theta =  9.224: L = 0.000e+00
  theta =  9.234: L = 1.892e-08
  theta =  9.734: L = 1.241e-08
  theta = 12.000: L = 2.326e-09

The grid searches agree with the closed forms (up to the grid spacing). The uniform case shows the shape of that likelihood: zero just below the sample maximum, largest exactly at it, and decreasing afterwards. Unlike the method-of-moments estimate $2\bar x$, the MLE can never contradict the data.

MLE vs. method of moments.

	Method of moments	Maximum likelihood
Uses	a few moments	the full model density
Computation	usually closed form	often numerical optimization
Large-sample efficiency	often below 1	optimal under regularity conditions
Respects the support	not always	yes
Sensitivity to model misspecification	lower (only moments must be right)	higher (the whole density is assumed)

Key takeaways

• The MLE is the parameter value under which the observed data are most probable.
• Maximize the log-likelihood: differentiate, set the score to zero, solve (numerically if needed), and check boundaries.
• Common results: $\hat p = \bar x$, $\hat\lambda = 1/\bar x$, $\hat\mu = \bar x$, $\hat\sigma^2$ with divisor $n$, uniform endpoint $= \max x_i$.
• When the support depends on the parameter, the maximum may be at a boundary where the derivative is not zero.

🧪 Practice

1. Derive the MLE of $\lambda$ for an i.i.d. Poisson sample. Does it coincide with the method-of-moments estimator?
2. Derive the MLE of $p$ for a Geometric$(p)$ sample (number of trials until the first success).
3. Observations $x_1, \dots, x_n$ come from the density $f(x; \theta) = \theta x^{\theta - 1}$ on $(0, 1)$. Find the MLE of $\theta$.
4. Interview: Explain why training a classifier with cross-entropy loss is maximum likelihood estimation. Hint: write the likelihood of the observed labels under a Bernoulli (or categorical) model whose probability is the model's output; take minus the log and average over examples.

Likelihood and Log-Likelihood Functions

Why it matters. The likelihood is more than a device for finding the maximum. The whole function describes what the data say about the parameter: its peak is the best estimate, its width tells you how precisely the data pin the parameter down, and comparing its height at two values measures how strongly the data favor one over the other. Reading likelihood curves correctly, and computing them without numerical disasters, is a basic skill for everything from confidence intervals to likelihood ratio tests (Section 9.6).

The likelihood is not a probability distribution of $\theta$. $L(\theta)$ is the probability (or density) of the fixed observed data evaluated as $\theta$ varies. As a function of $\theta$ it need not integrate to 1, and $L(0.3) = 2L(0.5)$ does not mean "$\theta = 0.3$ is twice as probable". It means "the observed data are twice as probable if $\theta = 0.3$ than if $\theta = 0.5$". Turning a likelihood into a probability statement about $\theta$ requires a prior (Bayes' theorem, Chapter 11).

Same formula $f(x; \theta)$ read as	Varies	Fixed	Sums/integrates to 1?
Probability model	$x$	$\theta$	yes, over $x$
Likelihood function	$\theta$	$x$ (the observed data)	not in general

Only ratios matter. Multiplying $L(\theta)$ by any constant that does not depend on $\theta$ changes nothing: the same maximizer and the same ratios. So the binomial coefficient can be dropped from the binomial likelihood, and 7 successes in 20 trials give the same likelihood whether you record the full sequence or only the count. This is why it is common to work with the relative likelihood

$$R(\theta) = \frac{L(\theta)}{L(\hat\theta)} \in [0, 1].$$

Why work on the log scale.

• Products become sums. $\ell(\theta) = \sum \log f(x_i; \theta)$: easier to differentiate, and each observation contributes an additive term.
• Numerical stability. A product of thousands of densities underflows to exactly 0 in floating point (the code shows this), while the sum of logs is an ordinary number. Every fitting routine works with log-likelihoods.
• Shape. For moderate to large $n$, the log-likelihood near its peak is close to a parabola, $\ell(\theta) \approx \ell(\hat\theta) - \frac{1}{2}J(\hat\theta)(\theta - \hat\theta)^2$. The curvature $J(\hat\theta) = -\ell''(\hat\theta)$, the observed information, measures how sharp the peak is (Fisher information, two topics ahead).

Likelihood intervals. The set of $\theta$ values with $R(\theta) \ge c$ is a likelihood interval: the parameter values that explain the data almost as well as the best one. With $c = e^{-1.92} \approx 0.147$, it is approximately a 95% confidence interval in large samples (because $2[\ell(\hat\theta) - \ell(\theta)]$ is approximately $\chi^2_1$, and $\chi^2_{1, 0.95} = 3.84$; Section 9.6).

import numpy as np

rng = np.random.default_rng(809)

# Products of many densities underflow; sums of logs do not
x = rng.normal(0, 1, size=2_000)
dens = np.exp(-x**2 / 2) / np.sqrt(2 * np.pi)
print(f"product of 2,000 densities: {np.prod(dens)}    sum of log densities: {np.log(dens).sum():.1f}")

# Relative likelihood L(p) / L(p_hat) for a proportion: same p_hat = 0.3, different n
def rel_lik(p, k, n):
    p_hat = k / n
    return np.exp(k * np.log(p / p_hat) + (n - k) * np.log((1 - p) / (1 - p_hat)))

p = np.linspace(0.01, 0.99, 99)
for k, n in [(3, 10), (30, 100)]:
    r = rel_lik(p, k, n)
    inside = p[r >= 0.147]                            # 0.147 = exp(-1.92): about a 95% likelihood interval
    print(f"\nk = {k}, n = {n}: likelihood interval p in [{inside.min():.2f}, {inside.max():.2f}]")
    for level in [0.9, 0.5, 0.147]:                   # crude ASCII plot: where is the curve above each level?
        print(f"{level:5.3f} |" + "".join("#" if v >= level else " " for v in r[::2]).rstrip())
    print("      +" + "-" * 50 + "\n      0" + " " * 48 + "1  p")

product of 2,000 densities: 0.0    sum of log densities: -2868.5

k = 3, n = 10: likelihood interval p in [0.09, 0.60]
0.900 |            ######
0.500 |        ################
0.147 |    ##########################
      +--------------------------------------------------
      0                                                1  p

k = 30, n = 100: likelihood interval p in [0.22, 0.39]
0.900 |              ##
0.500 |            ######
0.147 |           #########
      +--------------------------------------------------
      0                                                1  p

The product of 2,000 normal densities underflows to 0.0, while the log-likelihood is $-2868.5$. Both likelihoods peak at $\hat p = 0.3$, but with 100 observations the curve is much narrower: the 95% likelihood interval shrinks from $[0.09, 0.60]$ to $[0.22, 0.39]$, about $\sqrt{10}$ times narrower for ten times the data. The small-sample curve is also visibly asymmetric, which a symmetric interval $\hat p \pm 1.96,\text{SE}$ cannot capture (Section 8.3).

The likelihood principle. If two experiments produce proportional likelihood functions, they carry the same evidence about $\theta$. Example: "7 successes in 20 fixed trials" (binomial) and "trials until the 7th success took 20 trials" (negative binomial) give likelihoods proportional to $p^7(1-p)^{13}$. Bayesian and pure likelihood inference obey the principle; some frequentist procedures, such as p-values that depend on the stopping rule, do not. It is a point of genuine disagreement between schools of inference (Section 11.1).

Key takeaways

• $L(\theta)$ is the probability of the observed data viewed as a function of $\theta$; it is not a distribution over $\theta$.
• Only ratios of likelihoods carry meaning; constants not involving $\theta$ can be dropped.
• Always compute on the log scale: sums instead of products, no underflow, and a near-parabolic shape.
• The curvature at the peak measures precision; ${\theta : R(\theta) \ge 0.147}$ is an approximate 95% interval.

🧪 Practice

1. For 3 successes in 10 Bernoulli trials, compute $R(0.5)$. Is $p = 0.5$ inside the 0.147 likelihood interval?
2. Show that for an exponential sample the log-likelihood is $n\log\lambda - \lambda\sum x_i$, and compute the observed information $J(\hat\lambda)$. How does it depend on $n$ and $\hat\lambda$?
3. In the code, change the 0.147 cut-off to $e^{-0.5,z^2}$ with $z = 1.645$. What confidence level does the new interval approximate?
4. Interview: A logistic regression fit in your code returns a log-likelihood of $-\infty$ or NaN. What are the likely causes and fixes? Hint: predicted probabilities of exactly 0 or 1 for an observed outcome ($\log 0$), often from perfect separation or unscaled features; compute $\log$ probabilities with stable functions such as log-sum-exp or log-sigmoid, and consider regularization.

Properties of MLEs

Why it matters. Maximum likelihood dominates statistical practice because of a set of guarantees that hold under general conditions: it converges to the truth, its sampling distribution is approximately normal with a variance you can compute from the likelihood itself, no other well-behaved estimator does better in large samples, and it transforms consistently. Knowing these properties, and their limits, tells you when to trust the standard errors your software prints.

The main properties. Under regularity conditions (the model is identifiable, the support does not depend on $\theta$, the log-likelihood is smooth, and the true $\theta$ is in the interior of the parameter space):

Property	Statement
Consistency	$\hat\theta_n \to \theta$ in probability
Asymptotic normality	$\sqrt n(\hat\theta_n - \theta) \to N\big(0, 1/I(\theta)\big)$, i.e. $\hat\theta_n \approx N\big(\theta, \frac{1}{nI(\theta)}\big)$
Asymptotic efficiency	the asymptotic variance $\frac{1}{nI(\theta)}$ equals the Cramér-Rao lower bound
Invariance	the MLE of $g(\theta)$ is $g(\hat\theta)$, for any function $g$

Here $I(\theta)$ is the Fisher information per observation (next topic). In practice, the standard error of an MLE is computed as

$$\widehat{\operatorname{SE}}(\hat\theta) = \frac{1}{\sqrt{n I(\hat\theta)}} \quad \text{or} \quad \frac{1}{\sqrt{J(\hat\theta)}},$$

using the expected or observed information. This is what statistical software reports next to each coefficient, and it gives the Wald interval $\hat\theta \pm 1.96,\widehat{\operatorname{SE}}$.

Why it works (sketch). The score $\ell'(\theta)$ is a sum of $n$ i.i.d. terms with mean 0 and variance $I(\theta)$, so by the CLT it is approximately $N(0, nI(\theta))$. A first-order Taylor expansion of $\ell'(\hat\theta) = 0$ around the true $\theta$ gives $\hat\theta - \theta \approx \ell'(\theta)/(-\ell''(\theta))$, and $-\ell''(\theta)/n \to I(\theta)$ by the law of large numbers. Dividing the approximately normal score by $nI(\theta)$ gives the result.

What the properties do not promise.

• No small-sample guarantees. MLEs can be biased: $\hat\sigma^2$ has divisor $n$, and $\hat\lambda = 1/\bar X$ for the exponential has $E[\hat\lambda] = \frac{n}{n-1}\lambda$.
• Invariance does not preserve unbiasedness, but it does preserve maximum likelihood: the MLE of the exponential median $\ln 2/\lambda$ is $\ln 2/\hat\lambda$.
• Regularity failures. For Uniform$(0, \theta)$ the support depends on $\theta$: the MLE converges at rate $1/n$ (faster than usual) and its limit is not normal. When the number of parameters grows with $n$ (Neyman-Scott problem: one mean per pair of observations), the MLE of the variance stays inconsistent.
• Misspecification. If the model is wrong, the MLE converges to the parameter of the model closest to the truth (in Kullback-Leibler divergence), and the usual standard errors can be wrong; sandwich (robust) standard errors fix the latter.

import math
import numpy as np

rng = np.random.default_rng(810)

# 1. Small-sample bias of the exponential rate MLE 1/xbar: E = n/(n - 1) * lambda
lam = 2.0
for n in [5, 20, 200]:
    est = 1 / rng.exponential(1 / lam, size=(100_000, n)).mean(axis=1)
    print(f"n = {n:3d}: mean of lambda_hat = {est.mean():.3f} (theory {n / (n - 1) * lam:.3f})", end="  ")
    z = (est - lam) / (lam / np.sqrt(n))                 # standardize with the asymptotic SD
    print(f"P(|z| < 1.96) = {np.mean(np.abs(z) < 1.96):.3f}")

# 2. Invariance: the MLE of the median ln(2)/lambda is ln(2)/lambda_hat
t = rng.exponential(1 / lam, size=500)
print(f"MLE of median = ln2 / lambda_hat = {math.log(2) * t.mean():.4f}   (true {math.log(2) / lam:.4f})")

# 3. Efficiency: Gamma shape, method of moments vs maximum likelihood
def digamma(a, h=1e-5):                                  # numerical derivative of log-gamma
    return (math.lgamma(a + h) - math.lgamma(a - h)) / (2 * h)

def gamma_shape_mle(x):
    c = math.log(x.mean()) - np.log(x).mean()            # the MLE solves log(k) - digamma(k) = c
    lo, hi = 1e-3, 1e3
    for _ in range(60):                                  # bisection: log(k) - digamma(k) decreases in k
        mid = math.sqrt(lo * hi)
        lo, hi = (mid, hi) if math.log(mid) - digamma(mid) > c else (lo, mid)
    return math.sqrt(lo * hi)

k_true, n = 2.5, 50
samples = rng.gamma(k_true, 4.0, size=(3_000, n))
mom = samples.mean(axis=1) ** 2 / samples.var(axis=1)
mle = np.array([gamma_shape_mle(s) for s in samples])
print(f"Gamma shape, n = {n}: MoM mean {mom.mean():.3f} SD {mom.std():.3f}   "
      f"MLE mean {mle.mean():.3f} SD {mle.std():.3f}   efficiency of MoM {mle.var() / mom.var():.2f}")

n =   5: mean of lambda_hat = 2.501 (theory 2.500)  P(|z| < 1.96) = 0.868
n =  20: mean of lambda_hat = 2.105 (theory 2.105)  P(|z| < 1.96) = 0.926
n = 200: mean of lambda_hat = 2.010 (theory 2.010)  P(|z| < 1.96) = 0.947
MLE of median = ln2 / lambda_hat = 0.3478   (true 0.3466)
Gamma shape, n = 50: MoM mean 2.697 SD 0.628   MLE mean 2.639 SD 0.524   efficiency of MoM 0.69

With $n = 5$, the exponential rate MLE overestimates $\lambda = 2$ by 25%, exactly the factor $n/(n-1)$, and the normal approximation covers only 87% instead of 95%. By $n = 200$ the bias is half a percent and the normal approximation is accurate. The invariance property gives the MLE of the median without any new derivation. For the gamma shape, the MLE has to be computed numerically (bisection on the likelihood equation $\log k - \psi(k) = \log \bar x - \overline{\log x}$, where $\psi$ is the digamma function), but it is noticeably more precise: the method-of-moments estimator has an efficiency of about 0.69, so it wastes roughly 30% of the data.

Key takeaways

• Under regularity conditions, MLEs are consistent, asymptotically normal and asymptotically efficient.
• The standard error comes from the information: $\operatorname{SE} \approx 1/\sqrt{nI(\hat\theta)}$ or $1/\sqrt{J(\hat\theta)}$.
• Invariance: the MLE of $g(\theta)$ is $g(\hat\theta)$.
• In small samples MLEs can be biased and their normal approximation poor; the guarantees also fail when the support depends on $\theta$ or the number of parameters grows with $n$.

🧪 Practice

1. Using invariance, find the MLE of $P(X > 1)$ for an exponential sample with $\bar x = 0.8$.
2. Show that $\frac{n-1}{n}\hat\lambda$ is unbiased for the exponential rate. Is it still the MLE?
3. For Uniform$(0, \theta)$, simulate $n(\theta - \hat\theta)/\theta$ for large $n$ and compare its distribution with Exponential(1). Why is the scaling $n$ rather than $\sqrt n$?
4. Interview: Your logistic regression reports a coefficient with a huge standard error (e.g. 25 with SE 4,000). What might be going on? Hint: complete or quasi-complete separation sends the MLE to infinity, the log-likelihood becomes flat and the information near zero; check the data and consider penalized (Firth or ridge) estimation.

Fisher Information and the Cramer-Rao Bound

Why it matters. Before collecting data, it is useful to know the best precision any method could possibly achieve. If your estimator already reaches that limit, there is nothing left to improve; if it is far from it, a better method exists. The Fisher information quantifies how much a single observation tells you about a parameter, and the Cramér-Rao bound converts it into a lower limit on the variance of unbiased estimators. It also supplies the standard errors of MLEs and guides experimental design: choose the design that maximizes information.

The score. For one observation, the score is the derivative of the log-density with respect to the parameter:

$$U(\theta) = \frac{\partial}{\partial\theta}\log f(X; \theta).$$

It measures how sensitive the log-probability of the observation is to the parameter. Under regularity conditions, its mean is zero at the true parameter: $E[U(\theta)] = \int \frac{\partial f/\partial\theta}{f} f,dx = \frac{\partial}{\partial\theta}\int f,dx = \frac{\partial}{\partial\theta}1 = 0$.

Fisher information. The variance of the score:

$$I(\theta) = \operatorname{Var}\big(U(\theta)\big) = E\left[\left(\frac{\partial \log f(X; \theta)}{\partial\theta}\right)^2\right] = -E\left[\frac{\partial^2 \log f(X; \theta)}{\partial\theta^2}\right].$$

The second form (equal under regularity conditions) shows the geometric meaning: the expected curvature of the log-likelihood. A sharply peaked log-likelihood (high information) pins the parameter down well. For $n$ i.i.d. observations, information adds up: $I_n(\theta) = nI(\theta)$.

Model (one observation)	$I(\theta)$	Cramér-Rao bound for $n$ obs.	Attained by
Bernoulli$(p)$	$\frac{1}{p(1-p)}$	$\frac{p(1-p)}{n}$	$\hat p = \bar X$
Poisson$(\lambda)$	$\frac{1}{\lambda}$	$\frac{\lambda}{n}$	$\bar X$
Normal$(\mu, \sigma^2)$, $\mu$ unknown	$\frac{1}{\sigma^2}$	$\frac{\sigma^2}{n}$	$\bar X$
Normal, $\sigma^2$ unknown ($\mu$ known)	$\frac{1}{2\sigma^4}$	$\frac{2\sigma^4}{n}$	$\frac{1}{n}\sum (X_i - \mu)^2$
Exponential$(\lambda)$	$\frac{1}{\lambda^2}$	$\frac{\lambda^2}{n}$	not attained by an unbiased estimator for finite $n$

Cramér-Rao lower bound. If $T$ is an unbiased estimator of $\theta$ and regularity conditions hold, then

$$\operatorname{Var}(T) \ge \frac{1}{nI(\theta)}.$$

More generally, for an unbiased estimator of $g(\theta)$, $\operatorname{Var}(T) \ge \frac{g'(\theta)^2}{nI(\theta)}$, and for a biased estimator with bias $b(\theta)$, $\operatorname{Var}(T) \ge \frac{(1 + b'(\theta))^2}{nI(\theta)}$.

Proof (via Cauchy-Schwarz, Section 6.5). Let $U_n = \ell'(\theta)$ be the score of the whole sample, with mean 0 and variance $nI(\theta)$. Differentiating $E[T] = \theta$ with respect to $\theta$ under the integral sign gives $E[T,U_n] = 1$; since $E[U_n] = 0$, this is $\operatorname{Cov}(T, U_n) = 1$. Cauchy-Schwarz says $\operatorname{Cov}(T, U_n)^2 \le \operatorname{Var}(T)\operatorname{Var}(U_n)$, so $1 \le \operatorname{Var}(T) \cdot nI(\theta)$.

An unbiased estimator that attains the bound is called efficient; it is automatically the UMVUE. Equality in Cauchy-Schwarz requires $T$ to be a linear function of the score, which happens essentially only in exponential families. When no estimator attains the bound, the UMVUE has a strictly larger variance, and the bound is still a benchmark: MLEs reach it asymptotically.

import numpy as np

rng = np.random.default_rng(811)
lam, n, reps = 3.0, 25, 200_000

# Score of one Poisson observation: d/d lambda [x log(lambda) - lambda] = x / lambda - 1
x1 = rng.poisson(lam, size=reps)
score = x1 / lam - 1
print(f"E[score] ~ {score.mean():+.4f}   Var(score) ~ {score.var():.4f}   I(lambda) = 1/lambda = {1 / lam:.4f}")

# Two unbiased estimators of lambda: the sample mean and the sample variance
x = rng.poisson(lam, size=(reps, n))
crlb = lam / n                                            # 1 / (n I(lambda))
print(f"Cramer-Rao bound:  {crlb:.4f}")
print(f"Var(sample mean):  {x.mean(axis=1).var():.4f}   (mean {x.mean(axis=1).mean():.3f})")
print(f"Var(sample var.):  {x.var(axis=1, ddof=1).var():.4f}   (mean {x.var(axis=1, ddof=1).mean():.3f})")

# Observed information: curvature of the log-likelihood at the MLE gives the standard error
data = x[0]
loglik = lambda l: data.sum() * np.log(l) - n * l
l_hat, h = data.mean(), 1e-4
J = -(loglik(l_hat + h) - 2 * loglik(l_hat) + loglik(l_hat - h)) / h**2    # minus the second derivative
print(f"one sample: lambda_hat = {l_hat:.3f}, SE from curvature 1/sqrt(J) = {1 / np.sqrt(J):.4f}, "
      f"formula sqrt(lambda_hat / n) = {np.sqrt(l_hat / n):.4f}")

E[score] ~ -0.0015   Var(score) ~ 0.3339   I(lambda) = 1/lambda = 0.3333
Cramer-Rao bound:  0.1200
Var(sample mean):  0.1199   (mean 2.999)
Var(sample var.):  0.8671   (mean 2.999)
one sample: lambda_hat = 3.040, SE from curvature 1/sqrt(J) = 0.3487, formula sqrt(lambda_hat / n) = 0.3487

The simulated score has mean 0 and variance $1/\lambda$, as the theory says. Both the sample mean and the sample variance are unbiased for $\lambda$ (a Poisson's mean and variance are equal), but only the sample mean reaches the bound $\lambda/n = 0.12$; the sample variance has about seven times as much variance. The observed information, computed as the numerical curvature of the log-likelihood at its peak, gives the standard error without any formula, and here it matches $\sqrt{\hat\lambda/n}$ exactly, which is how software computes standard errors for models with no closed form.

Multiparameter version. With a parameter vector $\theta$, the Fisher information matrix has entries $I_{jk} = -E\left[\frac{\partial^2 \log f}{\partial\theta_j,\partial\theta_k}\right]$, and the covariance matrix of an unbiased estimator satisfies $\operatorname{Cov}(T) \succeq (nI)^{-1}$ (the difference is positive semidefinite). Standard errors of the MLE are the square roots of the diagonal of the inverse information matrix, so estimating extra nuisance parameters can only increase them.

Key takeaways

• The score has mean 0; its variance is the Fisher information $I(\theta)$, the expected curvature of the log-likelihood.
• Information adds over independent observations: $I_n = nI$.
• Cramér-Rao: no unbiased estimator has variance below $1/(nI(\theta))$; one that reaches it is efficient.
• The MLE attains the bound asymptotically, and $1/\sqrt{J(\hat\theta)}$ is its standard error.

🧪 Practice

1. Derive $I(p) = \frac{1}{p(1-p)}$ for a Bernoulli observation, using both the variance of the score and the expected second derivative.
2. For which $p$ is a Bernoulli observation most informative? Explain why that makes sense, and what it implies for the sample size needed to estimate a rare event rate with a given relative precision.
3. Compute the Cramér-Rao bound for unbiased estimators of $e^{-\lambda}$ from a Poisson sample, and compare it with the variance of the Rao-Blackwellized estimator of Section 8.1 (by simulation).
4. Interview: You can choose the input values $x$ at which to run an experiment for fitting $y = \beta_0 + \beta_1 x + \varepsilon$ within a range $[a, b]$. Where would you put them to estimate the slope most precisely, and what is the risk? Hint: $\operatorname{Var}(\hat\beta_1) = \sigma^2/\sum (x_i - \bar x)^2$, so information is maximized at the ends of the range; but you can no longer detect curvature in between.

Least Squares Estimation

Why it matters. Least squares is the oldest and most widely used method for fitting a model to data: Legendre and Gauss used it around 1805 to determine the orbits of comets and asteroids from noisy observations. It needs no full probability model, only a predicted value for each observation, and it has closed-form solutions for linear models. It is the engine of linear regression (Chapter 10), and it coincides with maximum likelihood when errors are normal.

The principle. For a model that predicts the value of observation $i$ as $g_i(\theta)$, choose $\theta$ to minimize the sum of squared residuals:

$$\hat\theta_{\text{LS}} = \arg\min_\theta \sum_{i=1}^n \big(y_i - g_i(\theta)\big)^2.$$

The simplest case is a single constant, $g_i(\mu) = \mu$: the minimizer of $\sum (y_i - \mu)^2$ is $\bar y$ (Section 2.2). Least squares generalizes the mean.

Simple linear regression. For $y_i = \beta_0 + \beta_1 x_i + \varepsilon_i$, setting the partial derivatives of $\sum (y_i - \beta_0 - \beta_1 x_i)^2$ to zero gives

$$\hat\beta_1 = \frac{\sum (x_i - \bar x)(y_i - \bar y)}{\sum (x_i - \bar x)^2}, \qquad \hat\beta_0 = \bar y - \hat\beta_1 \bar x.$$

In matrix form, with design matrix $X$ (a column of ones and a column of $x$ values), the normal equations are $X^\top X\hat\beta = X^\top y$, with solution $\hat\beta = (X^\top X)^{-1}X^\top y$ when $X^\top X$ is invertible. In code, never form the inverse explicitly; use a QR-based solver such as np.linalg.lstsq.

Geometric view. The fitted values $X\hat\beta$ are the orthogonal projection of $y$ onto the space spanned by the columns of $X$. Orthogonality means the residuals are uncorrelated with every predictor and (with an intercept) sum to zero.

              y
              *
              |\
     residual | \
     (minimal)|  \
              |   \
   -----------+----*--------------  plane spanned by the columns of X
         fitted value X beta_hat

Statistical justification.

Assumptions on the errors $\varepsilon_i$	What least squares gives
none (just a fitting criterion)	the best fit in squared error; a description of the data
mean 0	unbiased coefficients
mean 0, uncorrelated, equal variance	Gauss-Markov theorem: the best linear unbiased estimator (BLUE), with the smallest variance among unbiased estimators that are linear in $y$
i.i.d. $N(0, \sigma^2)$	also the MLE: the normal log-likelihood is $-\frac{1}{2\sigma^2}\sum (y_i - g_i(\theta))^2$ plus constants, so maximizing it is minimizing the sum of squares; exact $t$ and $F$ inference (Chapter 10)

The connection with maximum likelihood also tells you what to do when the assumptions change: with unequal error variances, weighted least squares (weights $1/\sigma_i^2$) is the MLE; with Laplace errors, the MLE minimizes the sum of absolute residuals (least absolute deviations, a robust alternative related to the median).

import numpy as np

rng = np.random.default_rng(812)

# Advertising spend (thousand EUR) vs weekly sales (thousand units): y = 3 + 0.8 x + noise
n = 40
x = rng.uniform(0, 20, n)
y = 3 + 0.8 * x + rng.normal(0, 2, n)

# Closed-form least squares
b1 = np.sum((x - x.mean()) * (y - y.mean())) / np.sum((x - x.mean()) ** 2)
b0 = y.mean() - b1 * x.mean()
print(f"closed form:   intercept {b0:.3f}, slope {b1:.3f}")

# Matrix form: solve the normal equations X^T X beta = X^T y (lstsq does it stably)
X = np.column_stack([np.ones(n), x])
beta, *_ = np.linalg.lstsq(X, y, rcond=None)
print(f"lstsq:         intercept {beta[0]:.3f}, slope {beta[1]:.3f}")
resid = y - X @ beta
print(f"residuals: sum {resid.sum():.1e}, correlation with x {np.corrcoef(resid, x)[0, 1]:.1e}, "
      f"sigma^2 estimate SSE/(n-2) = {resid @ resid / (n - 2):.2f}")

# Unbiasedness of the slope over many samples (same x values, new noise)
slopes = [np.polyfit(x, 3 + 0.8 * x + rng.normal(0, 2, n), 1)[0] for _ in range(20_000)]
print(f"slope over 20,000 samples: mean {np.mean(slopes):.4f}, SD {np.std(slopes):.4f} "
      f"(theory {2 / np.sqrt(np.sum((x - x.mean()) ** 2)):.4f})")

# Sensitivity to one gross error
y_bad = y.copy()
y_bad[np.argmax(x)] -= 40                                 # a data-entry error at the largest x
print(f"with one corrupted point: slope {np.polyfit(x, y_bad, 1)[0]:.3f}")

closed form:   intercept 3.145, slope 0.787
lstsq:         intercept 3.145, slope 0.787
residuals: sum -9.2e-14, correlation with x -7.7e-16, sigma^2 estimate SSE/(n-2) = 5.10
slope over 20,000 samples: mean 0.8004, SD 0.0567 (theory 0.0564)
with one corrupted point: slope 0.489

The closed form and the matrix solver agree. The residuals sum to zero and are uncorrelated with $x$, as the projection geometry requires, and $\text{SSE}/(n - 2)$ estimates the error variance $\sigma^2 = 4$ (two degrees of freedom are used by the two coefficients). Over 20,000 simulated datasets the slope averages to the true 0.8, with a spread that matches $\sigma/\sqrt{\sum (x_i - \bar x)^2}$. One corrupted point at the edge of the data cuts the slope from 0.79 to 0.49: squaring the residuals gives large errors a large influence, the same weakness the mean has (Section 10.3 covers diagnosing such points).

Key takeaways

• Least squares chooses parameters that minimize the sum of squared residuals; for a constant it gives the mean.
• Linear models have a closed-form solution via the normal equations; solve them with QR (lstsq), not an explicit inverse.
• Gauss-Markov: with uncorrelated, equal-variance, mean-zero errors, OLS is the best linear unbiased estimator.
• With normal errors, least squares is maximum likelihood; other error models lead to weighted or absolute-deviation fits.
• Like the mean, least squares is sensitive to outliers and high-leverage points.

🧪 Practice

1. Derive $\hat\beta_0$ and $\hat\beta_1$ by setting the two partial derivatives of the sum of squares to zero.
2. Fit the model $y = \beta x$ (no intercept) by least squares. Show that $\hat\beta = \sum x_i y_i / \sum x_i^2$, and explain why the residuals no longer need to sum to zero.
3. Show that if errors are i.i.d. Laplace, the MLE of a constant location parameter is the sample median.
4. Interview: Why might you prefer np.linalg.lstsq or a QR decomposition over computing $(X^\top X)^{-1}X^\top y$ directly, and when does it matter most? Hint: forming $X^\top X$ squares the condition number; with nearly collinear predictors (Section 10.3) the explicit inverse loses many digits of precision.

<a id="83-confidence-intervals"></a>

8.3 Confidence Intervals

A point estimate on its own hides how uncertain it is; a confidence interval reports a range of plausible parameter values, built by a procedure that captures the true value at a known long-run rate.

Interpretation of Confidence Level

Why it matters. "Average delivery time: 3.2 days, 95% CI [2.9, 3.5]" is one of the most common sentences in applied statistics, and one of the most commonly misread. Many readers take it as "there is a 95% chance the true mean is between 2.9 and 3.5". In the frequentist framework that produced the interval, that statement is not what the 95% means. Getting the interpretation right is essential before learning any of the formulas, because every interval in this subchapter carries the same meaning.

Definition. A $100(1 - \alpha)%$ confidence interval for $\theta$ is a pair of statistics $L(X)$ and $U(X)$, computed from the sample, such that

$$P_\theta\big(L(X) \le \theta \le U(X)\big) = 1 - \alpha \quad \text{for every value of } \theta.$$

The probability statement is about the random endpoints, before the data are seen. The parameter $\theta$ is fixed; the interval is what varies from sample to sample. $1 - \alpha$ is the confidence level (or nominal coverage), and $\alpha$ the error rate.

What 95% means. If you repeated the study many times, each time computing an interval by the same procedure, about 95% of those intervals would contain the true value. Once you have computed one specific interval, such as [2.9, 3.5], it either contains the true mean or it does not; you do not know which, and there is no randomness left in it. The 95% describes the reliability of the method, and your confidence in this particular interval is borrowed from that track record.

Analogy. A ring toss at a fair where the player hits the peg 95% of the time. Once a ring has been thrown, it is either on the peg or not, and if the peg is hidden behind a curtain, you are "95% confident" it landed, because of how reliable the thrower is. The peg (parameter) never moves; the ring (interval) does.

How intervals are built: pivots. Most intervals come from a pivotal quantity: a function of the data and the parameter whose distribution is known and does not depend on the parameter. For a normal mean with known $\sigma$,

$$Z = \frac{\bar X - \mu}{\sigma/\sqrt n} \sim N(0, 1), \quad \text{so} \quad P\left(-z_{\alpha/2} \le \frac{\bar X - \mu}{\sigma/\sqrt n} \le z_{\alpha/2}\right) = 1 - \alpha.$$

Rearranging the inequality to isolate $\mu$ gives the interval $\bar X \pm z_{\alpha/2},\sigma/\sqrt n$. The following topics apply this recipe with different pivots ($t$, $\chi^2$, approximately normal).

import numpy as np

rng = np.random.default_rng(8131)
mu, sigma, n = 50.0, 10.0, 25
half = 1.96 * sigma / np.sqrt(n)                       # half-width of a 95% interval (sigma known)

# 20 repetitions of the same study: each produces its own interval
print("each line is one study's 95% interval; | marks the true mean 50")
lo_axis, hi_axis, width = 40.0, 60.0, 60
to_col = lambda v: min(max(int(round((v - lo_axis) / (hi_axis - lo_axis) * width)), 0), width)  # clamp to the axis
for i in range(20):
    xbar = rng.normal(mu, sigma, n).mean()
    lo, hi = xbar - half, xbar + half
    line = [" "] * (width + 1)
    for c in range(to_col(lo), to_col(hi) + 1):
        line[c] = "-"
    line[to_col(mu)] = "|"
    miss = "" if lo <= mu <= hi else "   <- misses"
    print(f"{i + 1:3d} [{lo:5.2f}, {hi:5.2f}] " + "".join(line).rstrip() + miss)

# Long-run coverage: the 95% is a property of the procedure
xbars = rng.normal(mu, sigma, size=(100_000, n)).mean(axis=1)
print(f"coverage over 100,000 studies: {np.mean(np.abs(xbars - mu) <= half):.4f}")

each line is one study's 95% interval; | marks the true mean 50
  1 [46.79, 54.63]                     ----------|--------------
  2 [43.99, 51.83]             ------------------|-----
  3 [44.26, 52.10]              -----------------|------
  4 [43.55, 51.39]            -------------------|----
  5 [45.02, 52.86]                ---------------|---------
  6 [41.08, 48.92]    -------------------------  |   <- misses
  7 [44.26, 52.10]              -----------------|------
  8 [43.48, 51.32]           --------------------|----
  9 [47.15, 54.99]                      ---------|---------------
 10 [44.77, 52.61]               ----------------|--------
 11 [49.43, 57.27]                             --|----------------------
 12 [46.35, 54.19]                    -----------|-------------
 13 [45.54, 53.38]                  -------------|----------
 14 [44.38, 52.22]              -----------------|-------
 15 [44.79, 52.63]               ----------------|--------
 16 [45.77, 53.61]                  -------------|-----------
 17 [45.66, 53.50]                  -------------|----------
 18 [45.34, 53.18]                 --------------|----------
 19 [47.05, 54.89]                      ---------|---------------
 20 [48.83, 56.67]                           ----|--------------------
coverage over 100,000 studies: 0.9503

Twenty repetitions of the same study produce twenty different intervals; one of them (study 6) misses the true mean. In the long run, 95.03% of the intervals cover it. No single interval comes labeled as a hit or a miss, which is why the confidence level describes the procedure rather than any one interval.

Common misinterpretations.

Statement	Correct?
"95% of intervals constructed this way contain $\mu$."	yes
"There is a 95% probability that $\mu$ lies in [2.9, 3.5]."	no (in frequentist terms): $\mu$ is fixed, the probability is 0 or 1. A Bayesian credible interval (Section 11.2) makes this kind of statement, with a prior.
"95% of the data lie in the interval."	no: that would be a prediction or tolerance interval, much wider (Section 10.1)
"A replication has a 95% chance of an estimate inside this interval."	no: the probability is lower (about 83% on average for a mean)
"Values near the center of the interval are more plausible than values near the edges."	roughly yes, for intervals derived from a likelihood or a normal pivot
"If two intervals overlap, the difference is not significant."	no (CI for Differences of Means and Proportions, below)

The trade-offs.

• Level vs. width. Higher confidence requires a wider interval: 99% intervals are about 31% wider than 95% ones for a normal pivot (2.576 vs. 1.960). A 100% interval is useless (the whole parameter space).
• Sample size vs. width. The width shrinks like $1/\sqrt n$ (Margin of Error and Sample Size Determination, below).
• Nominal vs. actual coverage. The formulas rely on assumptions (normality, independence, large $n$). When those fail, the actual coverage can differ from the advertised level, sometimes badly. Checking coverage by simulation, as the code above does, is the honest test of an interval method.

Key takeaways

• The confidence level is the long-run proportion of intervals, from repeated samples, that contain the fixed parameter.
• A computed interval either contains $\theta$ or not; "95% probability that $\theta$ is in [a, b]" is a Bayesian statement, not a frequentist one.
• Intervals come from pivots: quantities with a known distribution that does not depend on $\theta$.
• Higher confidence means wider intervals; more data means narrower intervals; violated assumptions mean the real coverage may differ from the nominal one.

🧪 Practice

1. A report says: "The 95% CI for the mean salary is [52k, 58k], so 95% of employees earn between 52k and 58k." Explain the mistake.
2. In 200 independent studies, each reporting a 95% CI for its own parameter, how many intervals would you expect to miss? What is the probability that all 200 contain their parameter?
3. Modify the code to use 80% intervals. How many of the 20 intervals would you expect to miss, and how does the width change?
4. Interview: A stakeholder asks, "So is there a 95% chance that the true uplift is between 1% and 3%?" How do you answer in a way that is correct but still useful? Hint: distinguish the procedure's guarantee from a probability about this interval; mention that a Bayesian analysis with a weak prior gives a similar interval with exactly that interpretation.

CI for a Mean with Known Variance

Why it matters. The known-variance interval is the simplest case: the pivot is exactly standard normal, and every step of the construction is visible. It is used in practice when the measurement process has a long, stable track record (a calibrated instrument, an industrial process with years of quality-control data), and it is the template for the large-sample intervals used everywhere else.

The interval. For $X_1, \dots, X_n$ i.i.d. $N(\mu, \sigma^2)$ with $\sigma$ known:

$$\bar x \pm z_{\alpha/2}\frac{\sigma}{\sqrt n},$$

where $z_{\alpha/2}$ is the upper $\alpha/2$ quantile of the standard normal ($\Phi(z_{\alpha/2}) = 1 - \alpha/2$). The quantity $z_{\alpha/2},\sigma/\sqrt n$ is the margin of error.

Confidence level	80%	90%	95%	99%	99.9%
$z_{\alpha/2}$	1.282	1.645	1.960	2.576	3.291

One-sided bounds. When only one direction matters (a minimum fill volume, a maximum defect rate), put all of $\alpha$ in one tail: $\mu \ge \bar x - z_\alpha \sigma/\sqrt n$ is a $100(1 - \alpha)%$ lower confidence bound. For 95% this uses $z_{0.05} = 1.645$, not 1.96.

When the data are not normal. By the central limit theorem (Section 6.5), $\bar X$ is approximately normal for large $n$ whatever the population shape, so the same interval is approximately valid. "Large" depends on skewness: for symmetric data a handful of observations suffice; for strongly skewed data many more are needed.

import numpy as np
from statistics import NormalDist

rng = np.random.default_rng(814)

# Critical values z_{alpha/2} for common confidence levels
for level in [0.80, 0.90, 0.95, 0.99, 0.999]:
    print(f"{level:.1%} -> z = {NormalDist().inv_cdf(1 - (1 - level) / 2):.3f}")

# A filling machine: sigma = 2.0 ml is known from years of calibration; 16 bottles measured
sigma, n = 2.0, 16
x = rng.normal(501.2, sigma, n)
z = NormalDist().inv_cdf(0.975)
xbar = x.mean()
print(f"\nxbar = {xbar:.2f} ml, 95% CI = [{xbar - z * sigma / np.sqrt(n):.2f}, {xbar + z * sigma / np.sqrt(n):.2f}]")
print(f"one-sided 95% lower bound: mu >= {xbar - NormalDist().inv_cdf(0.95) * sigma / np.sqrt(n):.2f}")

# Non-normal data (exponential, sigma = 1 known): overall coverage is fine, but the misses are lopsided
print("\nnominal 95% z-interval: coverage, and how the 5% of misses split")
for n in [5, 30, 200]:
    for name, draw in [("normal", lambda s: rng.normal(1, 1, s)), ("exponential", lambda s: rng.exponential(1, s))]:
        xb = draw((100_000, n)).mean(axis=1)               # true mean 1 in both cases
        h = z / np.sqrt(n)
        below, above = np.mean(xb + h < 1), np.mean(xb - h > 1)   # interval entirely below / above mu
        print(f"  n = {n:3d}, {name:11s}: coverage {1 - below - above:.4f}   "
              f"interval too low {below:.4f}, too high {above:.4f}")

80.0% -> z = 1.282
90.0% -> z = 1.645
95.0% -> z = 1.960
99.0% -> z = 2.576
99.9% -> z = 3.291

xbar = 501.02 ml, 95% CI = [500.04, 502.00]
one-sided 95% lower bound: mu >= 500.20

nominal 95% z-interval: coverage, and how the 5% of misses split
  n =   5, normal     : coverage 0.9500   interval too low 0.0250, too high 0.0250
  n =   5, exponential: coverage 0.9555   interval too low 0.0004, too high 0.0441
  n =  30, normal     : coverage 0.9501   interval too low 0.0256, too high 0.0243
  n =  30, exponential: coverage 0.9526   interval too low 0.0140, too high 0.0335
  n = 200, normal     : coverage 0.9511   interval too low 0.0237, too high 0.0252
  n = 200, exponential: coverage 0.9496   interval too low 0.0216, too high 0.0287

For the filling machine, the 95% interval [500.04, 502.00] ml is entirely above 500 ml, and the one-sided bound says the mean is at least 500.20 ml. For exponential data, the overall coverage looks fine even with $n = 5$, but the misses are lopsided: the interval is almost never entirely below the true mean and too high 4.4% of the time, instead of 2.5% on each side. If you only care about one direction, as with a one-sided bound, the skewness matters a lot. The imbalance fades as $n$ grows and the CLT takes hold.

Key takeaways

• With known $\sigma$: $\bar x \pm z_{\alpha/2},\sigma/\sqrt n$; 1.645, 1.960 and 2.576 are the 90%, 95% and 99% critical values.
• One-sided bounds use $z_\alpha$ instead of $z_{\alpha/2}$.
• For non-normal data the interval is approximately valid by the CLT; skewness unbalances the two tails before it affects the overall coverage.
• Knowing $\sigma$ is rare; the next topic handles the usual case.

🧪 Practice

1. A scale has a known measurement SD of 0.4 g. Nine weighings of an object average 25.13 g. Compute the 90% and 99% confidence intervals.
2. Derive the interval from the pivot $Z$ by solving $-z_{\alpha/2} \le Z \le z_{\alpha/2}$ for $\mu$, step by step.
3. For the same data, why is the 95% one-sided lower bound closer to $\bar x$ than the lower end of the 95% two-sided interval?
4. Interview: A monitoring system computes an interval for the mean response time every minute, using $\sigma$ estimated from last year's data. What could go wrong? Hint: the "known" $\sigma$ may have changed (new hardware, traffic mix); response times are skewed and correlated in time, which breaks both the normal approximation and the $\sigma/\sqrt n$ standard error.

CI for a Mean with Unknown Variance

Why it matters. In almost every real problem the population standard deviation is unknown and must be estimated from the same data. Replacing $\sigma$ by $s$ adds uncertainty, and ignoring it gives intervals that are too narrow, especially for small samples. The $t$-interval is the standard tool for a mean, and it is also the model for the intervals of regression coefficients (Chapter 10).

The pivot. For i.i.d. normal data, $\frac{\bar X - \mu}{S/\sqrt n} \sim t_{n-1}$ (Section 5.4), which does not depend on $\mu$ or $\sigma$. The interval is

$$\bar x \pm t_{\alpha/2,,n-1}\frac{s}{\sqrt n},$$

where $s$ is the sample SD (divisor $n - 1$), $s/\sqrt n$ is the estimated standard error, and $t_{\alpha/2,,n-1}$ is the upper $\alpha/2$ quantile of the $t$ distribution with $n - 1$ degrees of freedom.

$n$	5	10	15	20	30	50	100	$\infty$
$t_{0.025,,n-1}$	2.776	2.262	2.145	2.093	2.045	2.010	1.984	1.960

With small samples the $t$ critical value is noticeably larger than 1.96, widening the interval to pay for not knowing $\sigma$.

Assumptions and robustness.

Assumption	How much it matters
Independent observations	critical: dependence (repeated measures, time series, clusters) makes $s/\sqrt n$ the wrong standard error
Normal population	important for small $n$, less so as $n$ grows (CLT); skewness matters most, heavy tails next
No gross outliers	a single outlier inflates both $\bar x$ and $s$; check a plot first

A practical rule: for $n < 15$, use the $t$-interval only if the data look roughly symmetric with no outliers; for $15 \le n < 40$, it is fine unless the data are strongly skewed; for $n \ge 40$, it usually works even with skewness, but very heavy tails (revenue, file sizes) can need thousands of observations, or the bootstrap (Section 8.4), or a transformation.

import numpy as np

rng = np.random.default_rng(815)
T975 = {9: 2.262, 29: 2.045, 99: 1.984}                 # t_{0.975, df} from a t table (Section 5.4)

# Battery life (hours) of 10 prototypes
x = np.array([18.2, 21.5, 19.8, 22.1, 17.6, 20.4, 23.0, 19.1, 20.9, 18.8])
n, xbar, s = len(x), x.mean(), x.std(ddof=1)            # ddof=1: the sample SD
half = T975[n - 1] * s / np.sqrt(n)
print(f"xbar = {xbar:.2f}, s = {s:.2f}, SE = {s / np.sqrt(n):.3f}")
print(f"95% t-interval: [{xbar - half:.2f}, {xbar + half:.2f}]   (z instead of t would give +-{1.96 * s / np.sqrt(n):.2f})")

# Coverage: t vs z critical values, normal vs skewed (lognormal) data
mu_ln = np.exp(0.5)                                     # mean of lognormal(0, 1)
print(f"\n{'n':>4s} {'normal, t':>10s} {'normal, z':>10s} {'lognormal, t':>13s}")
for n in [10, 30, 100]:
    xn = rng.normal(0, 1, (100_000, n))
    xl = rng.lognormal(0, 1, (100_000, n))
    cover = lambda d, c, m: np.mean(np.abs(d.mean(axis=1) - m) <= c * d.std(axis=1, ddof=1) / np.sqrt(n))
    print(f"{n:4d} {cover(xn, T975[n - 1], 0):10.4f} {cover(xn, 1.96, 0):10.4f} {cover(xl, T975[n - 1], mu_ln):13.4f}")

xbar = 20.14, s = 1.76, SE = 0.555
95% t-interval: [18.88, 21.40]   (z instead of t would give +-1.09)

   n  normal, t  normal, z  lognormal, t
  10     0.9484     0.9172        0.8403
  30     0.9510     0.9411        0.8824
 100     0.9499     0.9471        0.9162

The 10 prototypes give a 95% interval of [18.88, 21.40] hours, with a margin of error of 1.26 hours; using 1.96 instead of $t_{0.025, 9} = 2.262$ would understate it as 1.09. The simulation shows the cost of that shortcut: with $n = 10$ normal observations, the $z$-based interval covers only 91.7% of the time. The $t$-interval is exact for normal data at every $n$. For lognormal data, which are strongly right-skewed, even the $t$-interval falls short: 84% coverage at $n = 10$ and still only 92% at $n = 100$. For such data, a log transformation or a bootstrap interval is a better choice.

Paired data. When observations come in pairs (before/after on the same subject, two measurements on the same item), compute the differences $d_i$ and apply the one-sample $t$-interval to them. This uses the pairing to remove between-subject variability (Section 7.2), and it is the next topic's first special case.

Key takeaways

• With unknown $\sigma$, use $\bar x \pm t_{\alpha/2,,n-1},s/\sqrt n$; the $t$ quantile accounts for estimating $\sigma$.
• The interval is exact for normal data and approximately valid for moderately non-normal data with enough observations.
• Independence is the critical assumption; skewness and outliers are the main threats for small samples.
• For strongly skewed data, consider a transformation or the bootstrap rather than trusting the nominal level.

🧪 Practice

1. Twelve measurements have $\bar x = 4.83$ and $s = 0.62$. Compute a 95% confidence interval for $\mu$ ($t_{0.025, 11} = 2.201$).
2. Explain why, for the same data, the 95% $t$-interval is always wider than the 95% interval you would compute if $\sigma$ happened to equal $s$.
3. Run the code with $t_3$-distributed data instead of lognormal. Is the coverage too low, too high, or about right? Why is symmetry so helpful?
4. Interview: An analyst computes a 95% CI for average revenue per user from 50 users and gets [-3, 41] dollars. What does the negative lower limit tell you, and what would you do? Hint: the interval is symmetric but the data are highly skewed and bounded at 0; the normal approximation is poor; use more data, a bootstrap interval, or model the zero and positive parts separately.

CI for a Proportion

Why it matters. Conversion rates, defect rates, approval ratings, click-through rates and model accuracies are all proportions, and they are often estimated from samples where the event is rare or the sample is small. The textbook formula breaks down exactly in those situations, and the failure is not subtle: it can produce intervals with negative lower limits, intervals of zero width, and actual coverage far below the advertised 95%. Knowing which interval to use is a practical skill.

The Wald interval. From $\hat p \approx N\big(p, p(1-p)/n\big)$, plugging in $\hat p$ for $p$ in the standard error:

$$\hat p \pm z_{\alpha/2}\sqrt{\frac{\hat p(1 - \hat p)}{n}}.$$

It is simple, and it is what most introductory courses teach, but it has serious problems. When $\hat p = 0$ or $1$, the interval has width zero. Near 0 or 1 it can extend outside $[0, 1]$. And its coverage oscillates and is often far below nominal even for $n$ in the hundreds (Brown, Cai and DasGupta, 2001).

The Wilson score interval. Instead of plugging in $\hat p$, solve $|\hat p - p| \le z\sqrt{p(1-p)/n}$ for $p$ (a quadratic inequality). The result is

$$\frac{\hat p + \frac{z^2}{2n}}{1 + \frac{z^2}{n}} \pm \frac{z}{1 + \frac{z^2}{n}}\sqrt{\frac{\hat p(1 - \hat p)}{n} + \frac{z^2}{4n^2}}.$$

The center is pulled from $\hat p$ toward $1/2$, the interval always stays in $[0, 1]$, and the coverage is close to nominal across $p$ and $n$.

The Agresti-Coull interval. A simple approximation to Wilson: add $z^2/2 \approx 2$ successes and 2 failures, then apply the Wald formula to $\tilde p = (x + 2)/(n + 4)$ with $\tilde n = n + 4$. Easy to compute by hand, slightly conservative.

The Clopper-Pearson "exact" interval. Inverts the binomial distribution directly, using beta quantiles. Its coverage is guaranteed to be at least $1 - \alpha$ for every $p$, which makes it conservative (often 97% to 99% in practice). It is the choice when a guarantee is required, such as in regulatory or safety settings.

Zero events: the rule of three. If none of $n$ trials is a success, a 95% upper confidence bound for $p$ is approximately $3/n$ (from $(1 - p)^n = 0.05$, so $p \approx -\ln(0.05)/n \approx 3/n$). No adverse events in 300 patients means the rate could still plausibly be as high as 1%.

import numpy as np
from math import comb, sqrt

z = 1.96

def wald(k, n):
    p = k / n
    h = z * sqrt(p * (1 - p) / n)
    return p - h, p + h

def wilson(k, n):
    p = k / n
    center = (p + z**2 / (2 * n)) / (1 + z**2 / n)
    h = z / (1 + z**2 / n) * sqrt(p * (1 - p) / n + z**2 / (4 * n**2))
    return center - h, center + h

def agresti_coull(k, n):
    n_t = n + z**2                                    # add about 2 successes and 2 failures
    p_t = (k + z**2 / 2) / n_t
    h = z * sqrt(p_t * (1 - p_t) / n_t)
    return p_t - h, p_t + h

# A concrete case: 3 defective parts out of 40
for name, f in [("Wald", wald), ("Wilson", wilson), ("Agresti-Coull", agresti_coull)]:
    lo, hi = f(3, 40)
    print(f"{name:14s} [{lo:.3f}, {hi:.3f}]")
lo_w, hi_w = wald(0, 40)
print(f"0 of 40 defective: Wald [{lo_w:.3f}, {hi_w:.3f}], Wilson [0, {wilson(0, 40)[1]:.3f}], "
      f"rule of three [0, {3 / 40:.3f}]")

# EXACT coverage: sum the binomial probabilities of all outcomes k whose interval contains p
def coverage(f, n, p):
    return sum(comb(n, k) * p**k * (1 - p) ** (n - k) for k in range(n + 1) if f(k, n)[0] <= p <= f(k, n)[1])

n = 50
print(f"\nexact coverage of nominal 95% intervals, n = {n}")
print(f"{'p':>6s} {'Wald':>7s} {'Wilson':>7s} {'A-C':>7s}")
for p in [0.01, 0.03, 0.05, 0.10, 0.20, 0.35, 0.50]:
    print(f"{p:6.2f} {coverage(wald, n, p):7.3f} {coverage(wilson, n, p):7.3f} {coverage(agresti_coull, n, p):7.3f}")

Wald           [-0.007, 0.157]
Wilson         [0.026, 0.199]
Agresti-Coull  [0.019, 0.206]
0 of 40 defective: Wald [0.000, 0.000], Wilson [0, 0.088], rule of three [0, 0.075]

exact coverage of nominal 95% intervals, n = 50
     p    Wald  Wilson     A-C
  0.01   0.395   0.911   0.986
  0.03   0.781   0.937   0.983
  0.05   0.920   0.962   0.962
  0.10   0.879   0.970   0.970
  0.20   0.938   0.951   0.951
  0.35   0.945   0.963   0.963
  0.50   0.935   0.935   0.935

With 3 defects in 40 parts, the Wald interval has a negative lower limit, an impossible defect rate. With 0 defects it collapses to the single point [0, 0], claiming certainty from 40 parts. The coverage table is computed exactly by summing binomial probabilities, without simulation noise. At $p = 0.01$ the Wald interval's actual coverage is 40%, not 95%; it is still below 90% at $p = 0.10$ with 50 observations. Wilson stays near 95% throughout (with some oscillation, unavoidable for a discrete distribution), and Agresti-Coull is slightly conservative for small $p$.

Which interval to use.

Situation	Recommended
general use	Wilson
quick hand calculation	Agresti-Coull ("add two successes and two failures")
guaranteed minimum coverage required	Clopper-Pearson
zero successes, quick upper bound	rule of three, $3/n$
large $n$ and $p$ not near 0 or 1 (at least 10 successes and 10 failures)	Wald is acceptable, and all methods agree

Key takeaways

• The Wald interval $\hat p \pm z\sqrt{\hat p(1 - \hat p)/n}$ fails for small samples and rare events: zero width at $\hat p = 0$, limits outside $[0, 1]$, low coverage.
• The Wilson interval inverts the score test, stays in $[0, 1]$ and has near-nominal coverage; use it by default.
• Agresti-Coull approximates Wilson by adding about two successes and two failures; Clopper-Pearson guarantees at least nominal coverage.
• With zero events in $n$ trials, the 95% upper bound is about $3/n$.

🧪 Practice

1. In a survey, 412 of 1,000 respondents approve of a policy. Compute the Wald and Wilson 95% intervals. Do they differ meaningfully?
2. A new feature caused no crashes in 1,500 sessions. Give a 95% upper bound for the crash rate per session.
3. Derive the Wilson interval by solving $(\hat p - p)^2 = z^2 p(1 - p)/n$ for $p$.
4. Interview: A model's accuracy on a test set of 200 examples is 97%. Your colleague reports "97% ± 2.4%". What is wrong with that, and what would you report? Hint: the Wald interval reaches above 99% and is unreliable near 1; report the Wilson interval, and consider whether the test set is large enough to distinguish 97% from a competitor's 95%.

CI for a Variance

Why it matters. Sometimes the spread is the quantity of interest: the precision of a measuring instrument, the consistency of a manufacturing process, the volatility of an asset, the variability of a drug's effect across patients. A quality engineer who needs $\sigma$ to be below a specification limit needs an interval for $\sigma$, not just for the mean. The interval for a variance also teaches an important lesson: unlike the $t$-interval, it is not robust to non-normality.

The pivot. For i.i.d. $N(\mu, \sigma^2)$ data (Section 7.3),

$$\frac{(n - 1)S^2}{\sigma^2} \sim \chi^2_{n-1}.$$

With $\chi^2_{\alpha/2}$ and $\chi^2_{1-\alpha/2}$ the lower and upper quantiles of $\chi^2_{n-1}$, $P\left(\chi^2_{\alpha/2} \le \frac{(n-1)S^2}{\sigma^2} \le \chi^2_{1-\alpha/2}\right) = 1 - \alpha$. Solving for $\sigma^2$ flips the inequalities:

$$\left[\frac{(n - 1)s^2}{\chi^2_{1 - \alpha/2,,n-1}},; \frac{(n - 1)s^2}{\chi^2_{\alpha/2,,n-1}}\right].$$

The larger quantile produces the lower limit. Take square roots of both ends for an interval for $\sigma$.

Asymmetry. The chi-square distribution is right-skewed, so the interval is not symmetric around $s^2$: it extends further above than below. Small samples estimate variability very poorly, and the upper limit can be several times the lower limit.

$n$	10	20	30	100
95% CI for $\sigma$, as multiples of $s$	[0.69, 1.83]	[0.76, 1.46]	[0.80, 1.34]	[0.88, 1.16]

import numpy as np

rng = np.random.default_rng(817)
n = 20
chi2_lo, chi2_hi = 8.907, 32.852                       # 2.5% and 97.5% quantiles of chi-square, 19 df

# Thickness of 20 machined parts (mm): is the process precise enough?
x = rng.normal(5.00, 0.03, n)
s2 = x.var(ddof=1)
lo, hi = (n - 1) * s2 / chi2_hi, (n - 1) * s2 / chi2_lo   # note: the LARGE quantile gives the LOWER limit
print(f"s = {np.sqrt(s2):.4f} mm; 95% CI for sigma^2: [{lo:.2e}, {hi:.2e}]; for sigma: [{np.sqrt(lo):.4f}, {np.sqrt(hi):.4f}]")
print(f"asymmetric around s: {np.sqrt(s2) - np.sqrt(lo):.4f} below, {np.sqrt(hi) - np.sqrt(s2):.4f} above")

# Coverage depends heavily on normality (all three populations have variance 1)
pops = {
    "normal":            lambda size: rng.normal(0, 1, size),
    "uniform":           lambda size: rng.uniform(-np.sqrt(3), np.sqrt(3), size),        # light tails
    "t, 5 df (scaled)":  lambda size: rng.standard_t(5, size) / np.sqrt(5 / 3),          # heavy tails
    "exponential":       lambda size: rng.exponential(1, size),
}
for name, draw in pops.items():
    s2 = draw((100_000, n)).var(axis=1, ddof=1)
    cover = np.mean(((n - 1) * s2 / chi2_hi <= 1) & (1 <= (n - 1) * s2 / chi2_lo))
    print(f"coverage, {name:17s}: {cover:.3f}")

s = 0.0248 mm; 95% CI for sigma^2: [3.56e-04, 1.31e-03]; for sigma: [0.0189, 0.0362]
asymmetric around s: 0.0059 below, 0.0114 above
coverage, normal           : 0.950
coverage, uniform          : 0.996
coverage, t, 5 df (scaled) : 0.842
coverage, exponential      : 0.729

With 20 parts, the standard deviation is known only to within about $-24%$ and $+46%$ of the sample value. The coverage simulation is the important result: all four populations have variance 1, yet the interval's real coverage ranges from 99.6% (light-tailed uniform) to 73% (exponential). The reason is that $\operatorname{Var}(S^2)$ depends on the kurtosis of the population, $\operatorname{Var}(S^2) \approx \frac{\sigma^4}{n}\left(\kappa - 1 + \frac{2}{n-1}\right)$ with $\kappa$ the kurtosis (3 for the normal), and the chi-square pivot assumes $\kappa = 3$. Unlike the mean, the CLT does not rescue this interval: the error does not shrink as $n$ grows.

Alternatives for non-normal data. Bootstrap intervals for the variance or SD (Section 8.4), or intervals based on a kurtosis-adjusted standard error. To compare variances of two groups, Levene's or the Brown-Forsythe test is preferred over the normal-theory $F$-test for the same reason (Section 9.2).

Key takeaways

• For normal data, $(n-1)S^2/\sigma^2 \sim \chi^2_{n-1}$ gives $\left[(n-1)s^2/\chi^2_{\text{upper}},, (n-1)s^2/\chi^2_{\text{lower}}\right]$.
• The interval is asymmetric, and wide for small samples: variability is hard to estimate.
• For an interval for $\sigma$, take square roots of the endpoints.
• The interval is very sensitive to non-normality (kurtosis), even for large $n$; prefer the bootstrap when normality is doubtful.

🧪 Practice

1. A sample of 10 measurements has $s = 0.8$. Compute a 95% CI for $\sigma$ using $\chi^2_{0.025, 9} = 2.700$ and $\chi^2_{0.975, 9} = 19.023$.
2. A specification requires $\sigma \le 0.05$ mm. From 25 parts, $s = 0.042$. Using a one-sided 95% upper bound ($\chi^2_{0.05, 24} = 13.848$), can you claim the process meets the specification?
3. Explain why the variance interval is accurate for uniform data in the sense of covering too often, while it under-covers for the exponential.
4. Interview: A finance team reports a 95% CI for the daily volatility of a stock using the chi-square formula on 250 daily returns. What would you warn them about? Hint: daily returns have heavy tails (high kurtosis) and volatility clusters in time; both make the chi-square interval too narrow; a block bootstrap is more honest.

CI for Differences of Means and Proportions

Why it matters. Most practical questions are comparisons: is the new checkout flow better than the old one, does the drug lower blood pressure more than placebo, is plant A more productive than plant B? The parameter of interest is a difference, and its confidence interval answers two questions at once. It tells you whether there is evidence of a difference (does the interval exclude 0?), and it tells you how large that difference plausibly is, which is usually what matters for decisions.

The general form. For independent estimates, variances add:

$$(\hat\theta_1 - \hat\theta_2) \pm (\text{critical value}) \times \sqrt{\operatorname{SE}_1^2 + \operatorname{SE}_2^2}.$$

Two independent means.

Method	Standard error	Degrees of freedom	When
Welch (default)	$\sqrt{\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2}}$	Welch-Satterthwaite (below)	variances may differ; always safe
Pooled	$s_p\sqrt{\frac{1}{n_1} + \frac{1}{n_2}}$, $s_p^2 = \frac{(n_1 - 1)s_1^2 + (n_2 - 1)s_2^2}{n_1 + n_2 - 2}$	$n_1 + n_2 - 2$	variances known to be similar
Large samples	as Welch	use $z$	both $n$ large

The Welch-Satterthwaite degrees of freedom are

$$\nu = \frac{\left(\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2}\right)^2}{\frac{(s_1^2/n_1)^2}{n_1 - 1} + \frac{(s_2^2/n_2)^2}{n_2 - 1}},$$

usually not an integer; software uses it directly, and a $t$ table is read at the integer below it (slightly conservative). The pooled interval can be badly wrong when the variances differ and the group sizes differ, while Welch loses almost nothing when the variances are equal, which is why it is the default in modern practice.

Paired means. When each observation in one group is matched with one in the other, the two samples are not independent, and the formula above is wrong. Compute the differences $d_i$ and use the one-sample interval $\bar d \pm t_{\alpha/2,,n-1},s_d/\sqrt n$.

Two proportions.

$$(\hat p_1 - \hat p_2) \pm z_{\alpha/2}\sqrt{\frac{\hat p_1(1 - \hat p_1)}{n_1} + \frac{\hat p_2(1 - \hat p_2)}{n_2}}.$$

This Wald-type interval works when each group has at least about 10 successes and 10 failures. For small counts, the Newcombe interval (built from two Wilson intervals) or the Agresti-Caffo adjustment (add one success and one failure to each group) has better coverage.

import numpy as np

rng = np.random.default_rng(818)

# 1. Two independent groups, unequal variances: Welch interval
a = rng.normal(52, 6, 25)        # time to first byte (ms), old CDN
b = rng.normal(47, 12, 18)       # new CDN: faster on average, but more variable
va, vb = a.var(ddof=1) / len(a), b.var(ddof=1) / len(b)
se = np.sqrt(va + vb)
df = (va + vb) ** 2 / (va**2 / (len(a) - 1) + vb**2 / (len(b) - 1))   # Welch-Satterthwaite
t_crit = 2.069                                                         # t_{0.975, 23}: df rounded down (t table)
diff = a.mean() - b.mean()
print(f"difference {diff:.2f}, SE {se:.2f}, Welch df {df:.1f}")
print(f"95% CI for mu_A - mu_B: [{diff - t_crit * se:.2f}, {diff + t_crit * se:.2f}]")

# 2. Paired data: the same 12 users measured before and after; analyze the differences
before = rng.normal(100, 15, 12)
after = before - 4 + rng.normal(0, 3, 12)          # a real drop of 4, small within-user noise
d = after - before
half = 2.201 * d.std(ddof=1) / np.sqrt(len(d))     # t_{0.975, 11}
print(f"\npaired:  mean diff {d.mean():.2f}, 95% CI [{d.mean() - half:.2f}, {d.mean() + half:.2f}]")
se_wrong = np.sqrt(before.var(ddof=1) / 12 + after.var(ddof=1) / 12)
print(f"treated as independent (wrong): SE {se_wrong:.2f} instead of {d.std(ddof=1) / np.sqrt(12):.2f}")

# 3. Two proportions: A/B test conversion
x1, n1, x2, n2 = 230, 4_000, 274, 4_000
p1, p2 = x1 / n1, x2 / n2
se = np.sqrt(p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2)
print(f"\nconversion A {p1:.2%}, B {p2:.2%}; 95% CI for p_B - p_A: "
      f"[{p2 - p1 - 1.96 * se:+.2%}, {p2 - p1 + 1.96 * se:+.2%}]")

# 4. The overlap fallacy: individual CIs overlap, yet the CI for the difference excludes 0
m1, m2, se1, se2 = 10.0, 13.0, 1.0, 1.0
print(f"\nCI1 [{m1 - 1.96 * se1:.2f}, {m1 + 1.96 * se1:.2f}], CI2 [{m2 - 1.96 * se2:.2f}, {m2 + 1.96 * se2:.2f}] overlap; "
      f"CI for difference [{m2 - m1 - 1.96 * np.hypot(se1, se2):.2f}, {m2 - m1 + 1.96 * np.hypot(se1, se2):.2f}]")

difference 7.17, SE 2.99, Welch df 23.2
95% CI for mu_A - mu_B: [1.00, 13.35]

paired:  mean diff -4.11, 95% CI [-6.48, -1.74]
treated as independent (wrong): SE 7.47 instead of 1.08

conversion A 5.75%, B 6.85%; 95% CI for p_B - p_A: [+0.04%, +2.16%]

CI1 [8.04, 11.96], CI2 [11.04, 14.96] overlap; CI for difference [0.23, 5.77]

The Welch interval for the CDN comparison excludes 0, so the data support a real difference, but its width (from 1 to 13 ms) shows that the size of the improvement is poorly determined. In the paired example, treating the before/after measurements as independent samples gives a standard error of 7.47 instead of 1.08, because it leaves the large between-user variability in the noise. The paired analysis detects the 4-unit drop easily; the unpaired one would not. In the A/B test, B's conversion is higher by between 0.04 and 2.16 percentage points: the difference is probably real but could be negligible. The last line shows the overlap fallacy: two 95% intervals can overlap while the interval for the difference excludes 0, because the SE of a difference is $\sqrt{2},\text{SE}$, not $2,\text{SE}$, when the SEs are equal.

Key takeaways

• For independent groups, the SE of a difference is the square root of the sum of squared SEs.
• Use the Welch interval for two means by default; use the paired interval on differences when observations are matched.
• For two proportions, the Wald difference interval needs about 10 successes and 10 failures per group; otherwise use Newcombe or Agresti-Caffo.
• Read both whether the interval excludes 0 and how wide it is; never compare two separate intervals by eye.

🧪 Practice

1. Group A: $n = 30$, $\bar x = 72.1$, $s = 8.4$; group B: $n = 35$, $\bar x = 67.5$, $s = 9.9$. Compute the Welch SE and degrees of freedom, and a 95% interval using the $z$ approximation. Would the $t$ value change the conclusion?
2. Show that when $n_1 = n_2$, the pooled and Welch standard errors are equal. Where do the two methods still differ?
3. Two non-overlapping 95% CIs: what can you say about the difference? Two overlapping ones: what can you say? Use the equal-SE case to quantify how much overlap is compatible with a significant difference.
4. Interview: An A/B test shows conversion of 5.0% vs. 5.3%, with a 95% CI for the difference of [-0.1%, +0.7%]. The product manager says the test "failed". How would you explain the result and what would you recommend? Hint: the interval is compatible with no effect and with a meaningful positive effect; the test was underpowered for effects of that size; discuss the smallest effect worth detecting and the sample size needed (next topic and Section 9.1).

Margin of Error and Sample Size Determination

Why it matters. "How much data do we need?" is one of the first questions in any study, and it should be answered before the data are collected. Too few observations produce intervals too wide to support any decision; too many waste money, time and participants' effort. Inverting the confidence interval formula answers the question directly, and the same reasoning explains the familiar "±3 percentage points" attached to opinion polls.

Margin of error. The margin of error (ME) is the half-width of a symmetric confidence interval:

$$E = z_{\alpha/2} \times \operatorname{SE}.$$

For a mean, $E = z_{\alpha/2},\sigma/\sqrt n$; for a proportion, $E = z_{\alpha/2}\sqrt{p(1 - p)/n}$.

Solving for $n$. Fix the desired margin $E$ and confidence level, then solve:

$$n = \left(\frac{z_{\alpha/2},\sigma}{E}\right)^2 \quad \text{(mean)}, \qquad n = \frac{z_{\alpha/2}^2,p(1 - p)}{E^2} \quad \text{(proportion)}.$$

Always round up. The formulas need a planning value for the unknown $\sigma$ or $p$:

• For $\sigma$: a pilot study, previous studies, or the range rule $\sigma \approx \text{range}/4$. Use a conservatively large value: underestimating $\sigma$ undersizes the study.
• For $p$: $p(1 - p)$ is largest at $p = 0.5$ (value $0.25$), so $n = z^2/(4E^2)$ is a worst-case size that works for any $p$. If you know $p$ is near 0.1, using it reduces $n$ substantially.

Finite populations. When the sample is a sizable fraction of a population of size $N$, apply the finite population correction (Section 7.3): $n = \frac{n_0}{1 + (n_0 - 1)/N}$, where $n_0$ is the size from the formula above.

The square-root law. Because $E \propto 1/\sqrt n$, halving the margin of error requires four times the sample, and one more decimal place of precision requires a hundred times the sample.

import math
import numpy as np

rng = np.random.default_rng(819)
z = 1.96

def n_for_mean(sigma, E):
    return math.ceil((z * sigma / E) ** 2)                  # always round UP

def n_for_proportion(E, p=0.5, N=None):
    n0 = z**2 * p * (1 - p) / E**2
    if N is not None:
        n0 = n0 / (1 + (n0 - 1) / N)                        # finite population correction
    return math.ceil(n0)

print("mean, sigma = 15:", {E: n_for_mean(15, E) for E in [5, 2, 1, 0.5]})
print("proportion, worst case p = 0.5:", {E: n_for_proportion(E) for E in [0.05, 0.03, 0.02, 0.01]})
print("proportion near 0.1:", {E: n_for_proportion(E, p=0.1) for E in [0.05, 0.03, 0.02, 0.01]})
print("E = 0.03 in a population of N = 2,000:", n_for_proportion(0.03, N=2_000))

# Check by simulation: with n = 1,068 and p = 0.5, how often is the poll within 3 points?
n = n_for_proportion(0.03)
p_hat = rng.binomial(n, 0.5, size=100_000) / n
print(f"n = {n}: P(|p_hat - p| <= 0.03) = {np.mean(np.abs(p_hat - 0.5) <= 0.03):.3f}")

# Margin of error for a subgroup: 20% of the sample (e.g. one region)
print(f"subgroup of {n // 5}: margin of error = +-{z * math.sqrt(0.25 / (n // 5)):.3f}")

mean, sigma = 15: {5: 35, 2: 217, 1: 865, 0.5: 3458}
proportion, worst case p = 0.5: {0.05: 385, 0.03: 1068, 0.02: 2401, 0.01: 9604}
proportion near 0.1: {0.05: 139, 0.03: 385, 0.02: 865, 0.01: 3458}
E = 0.03 in a population of N = 2,000: 697
n = 1068: P(|p_hat - p| <= 0.03) = 0.953
subgroup of 213: margin of error = +-0.067

The familiar numbers appear: about 385 respondents for ±5 points, 1,068 for ±3 points and 9,604 for ±1 point, whatever the population size (a national poll and a city poll need about the same $n$). If the proportion is known to be near 10%, the sample needed falls to 36% of the worst case. For a population of only 2,000, the finite population correction reduces 1,068 to 697. The simulation confirms that a poll of 1,068 lands within 3 points 95% of the time. The last line is often forgotten: a subgroup of 213 respondents has a margin of error of almost ±7 points, so breakdowns by region or age in a national poll are much less precise than the headline number.

What the margin of error does not cover. It measures only sampling error, the random variation from drawing a random sample. Nonresponse, undercoverage, question wording and other nonsampling errors (Section 7.1) are not included, and in real surveys they are often larger. A poll with a ±3% margin of error can easily be off by more than 3 points.

Sample size for comparisons. For a difference between two groups, the SE involves two variances, so each group needs roughly twice the single-group size for the same margin on the difference. When the goal is to detect an effect with a test, the sample size also depends on the desired power (Section 9.1).

Key takeaways

• Margin of error = critical value × standard error; it is the half-width of the interval.
• $n = (z\sigma/E)^2$ for a mean and $n = z^2 p(1 - p)/E^2$ for a proportion; round up, and use conservative planning values ($p = 0.5$ in the worst case).
• Halving the margin of error requires four times the data; population size barely matters unless the sample is a large fraction of it.
• The margin of error covers only sampling error, and subgroups have much larger margins than the full sample.

🧪 Practice

1. How many measurements are needed to estimate a mean within ±0.5 units with 99% confidence if $\sigma \approx 3$?
2. A poll of 600 people finds 54% support. Compute the margin of error. Can the poll claim a majority supports the measure?
3. A company with 800 employees wants a ±5-point margin of error for an engagement survey proportion. How many responses are needed with the finite population correction? What if only 60% of invited employees respond?
4. Interview: A stakeholder asks why a survey of 1,000 people can represent a country of 60 million. How do you explain it? Hint: the standard error depends on $n$, not on $N$, when the sample is random and a small fraction of the population (like tasting soup: a spoonful from a well-stirred pot works for any pot size); the real risks are non-random selection and nonresponse.

<a id="84-resampling-methods"></a>

8.4 Resampling Methods

Resampling methods replace difficult mathematics with computation: by reusing the observed data many times, in carefully chosen ways, they estimate standard errors, bias, confidence intervals, null distributions and prediction error for statistics that have no convenient formula.

Bootstrap Principle

Why it matters. The formulas of Section 8.3 cover means, proportions and variances under specific assumptions. But what is the standard error of a median, a 90th percentile, a correlation coefficient, a ratio of two means, or the AUC of a classifier? Deriving each one is hard or impossible, and the asymptotic answers may be poor for your sample size. The bootstrap (Efron, 1979) gives a single, general recipe that works for almost any statistic, using nothing but the data and a computer. It is one of the most important ideas in modern statistics.

The problem. The standard error of a statistic $\hat\theta$ is the SD of its sampling distribution, the distribution we would see if we could draw many samples from the population (Section 7.3). We cannot: we have one sample.

The idea: the plug-in principle. We do not know the population distribution $F$, but we have a good estimate of it: the empirical distribution $\hat F_n$, which puts probability $1/n$ on each observed data point. So we pretend that the sample is the population, and simulate the sampling process from it. Drawing $n$ points from $\hat F_n$ means drawing from the data with replacement.

        REAL WORLD (unknowable)                   BOOTSTRAP WORLD (computable)

   population F  --sample-->  data x         data x = F_hat  --resample-->  x*_1, x*_2, ..., x*_B
                               |                                              |     |          |
                         theta_hat = s(x)                           theta*_b = s(x*_b) for each b
                               |                                              |
   sampling distribution of theta_hat       ~~~~     bootstrap distribution of theta*
   (what we want; needs many samples)         approximates  (spread of theta* around theta_hat)

The algorithm.

1. From the data $x_1, \dots, x_n$, draw a bootstrap resample $x_1^, \dots, x_n^$ of the same size, with replacement.
2. Compute the statistic on it: $\hat\theta^*_b$.
3. Repeat $B$ times (typically 1,000 to 10,000).
4. Use the distribution of $\hat\theta^_1, \dots, \hat\theta^_B$ as an approximation of the sampling distribution of $\hat\theta$:
   • bootstrap SE: $\widehat{\operatorname{SE}}_{\text{boot}} = \operatorname{SD}(\hat\theta^*_b)$;
   • bootstrap bias: $\widehat{\operatorname{Bias}}_{\text{boot}} = \overline{\hat\theta^*} - \hat\theta$;
   • confidence intervals: next topic.

Why "with replacement"? Sampling $n$ values without replacement from $n$ data points just returns the same dataset in a different order. Sampling with replacement creates new datasets in which some points appear several times and others not at all. On average a resample contains about $1 - (1 - 1/n)^n \approx 1 - 1/e \approx 63%$ of the distinct original points.

Two layers of approximation. The bootstrap answer differs from the truth for two reasons:

Error source	Cause	How to reduce it
Statistical error	$\hat F_n$ is not $F$: the sample is not the population	more data (only)
Monte Carlo error	finite number $B$ of resamples	more resamples (cheap)

Increasing $B$ removes only the second. The bootstrap cannot create information that the sample does not contain.

import numpy as np

rng = np.random.default_rng(820)

# The "real world": a skewed population we normally cannot see (e.g. session lengths in minutes)
population = rng.lognormal(mean=1.5, sigma=0.8, size=1_000_000)
n, B = 50, 5_000
stats = {"mean": np.mean, "median": np.median, "SD": lambda a, axis=None: np.std(a, ddof=1, axis=axis)}

def bootstrap(sample, f, B):
    idx = rng.integers(0, len(sample), size=(B, len(sample)))   # B resamples of size n, as index arrays
    return f(sample[idx], axis=1)                               # the statistic on every resample

sample = rng.choice(population, size=n)                         # the ONE sample we actually observe
print(f"{'statistic':10s} {'estimate':>9s} {'true SE':>8s} {'bootstrap SE':>13s} {'boot SE, 200 other samples':>28s}")
for name, f in stats.items():
    true_se = f(rng.choice(population, size=(B, n)), axis=1).std()     # needs the population: impossible in practice
    boot_se = bootstrap(sample, f, B).std()
    others = [bootstrap(rng.choice(population, size=n), f, 1_000).std() for _ in range(200)]
    q10, q50, q90 = np.quantile(others, [0.1, 0.5, 0.9])
    print(f"{name:10s} {f(sample):9.2f} {true_se:8.3f} {boot_se:13.3f}     median {q50:.3f} [{q10:.3f}, {q90:.3f}]")

# How many distinct observations does a typical resample contain?
idx = rng.integers(0, n, size=(2_000, n))
share = np.mean([len(np.unique(row)) for row in idx]) / n
print(f"\naverage share of distinct original points in a resample: {share:.3f}  (theory 1 - (1 - 1/n)^n = {1 - (1 - 1 / n) ** n:.3f})")

# A failure case: the bootstrap distribution of the sample maximum
boot_max = bootstrap(sample, np.max, B)
print(f"P(bootstrap max == sample max) = {np.mean(boot_max == sample.max()):.3f}: a point mass at the observed max")

statistic   estimate  true SE  bootstrap SE   boot SE, 200 other samples
mean            5.03    0.832         0.594     median 0.706 [0.505, 1.048]
median          4.33    0.638         0.727     median 0.646 [0.435, 0.936]
SD              4.21    1.725         1.218     median 0.876 [0.414, 1.895]

average share of distinct original points in a resample: 0.636  (theory 1 - (1 - 1/n)^n = 0.636)
P(bootstrap max == sample max) = 0.633: a point mass at the observed max

From a single sample, and without any formula, the bootstrap gives standard errors of the right order of magnitude for all three statistics, including the median and the SD, whose formulas are awkward for skewed data. The match is not perfect (0.59 against a true 0.83 for the mean), and the last column shows why: the bootstrap SE is itself an estimate, computed from one sample, and it varies from sample to sample. For skewed data it also tends to run low: the typical bootstrap SE of the mean is 0.71 against a true 0.83, and for the SD the estimates range from 0.41 to 1.90 (10th to 90th percentile) around a true 1.73. A sample of 50 rarely contains the rare large values that drive the population's variability. The resamples contain 63.6% distinct points, as the theory predicts. The final line is a classic failure: 63% of the bootstrap maxima equal the observed maximum exactly, so the bootstrap distribution of the maximum is a poor approximation of its true, continuous sampling distribution.

When the bootstrap works and when it does not.

Works well	Fails or needs modification
smooth statistics: means, variances, correlations, regression coefficients, ratios	extremes (max, min) and statistics at the boundary of the parameter space
quantiles such as the median (for SE; with moderate $n$)	very small samples ($n < 10$ or so): $\hat F_n$ is a poor stand-in for $F$
i.i.d. data	dependent data: use the block bootstrap (time series) or resample whole clusters
moderate skewness	infinite-variance populations (very heavy tails)

The resampling must mimic how the data were generated: resample subjects, not individual measurements, when there are repeated measurements per subject; resample within strata for a stratified design; resample blocks of consecutive observations for time series.

Key takeaways

• The bootstrap treats the sample as the population and simulates sampling by drawing $n$ points with replacement.
• The spread of the statistic over $B$ resamples estimates its standard error; the shift of their mean estimates its bias.
• It works for almost any smooth statistic, without formulas, but cannot add information: its accuracy is limited by the sample.
• It fails for extremes and tiny samples, and must respect the data's structure (clusters, strata, time dependence).

🧪 Practice

1. For the data ${2, 5, 9}$, list all $3^3 = 27$ equally likely ordered resamples' means and compute the exact bootstrap SE of the mean. Compare it with $s/\sqrt n$ and with $\sqrt{(n-1)/n};s/\sqrt n$.
2. Show that the probability that a given observation is not in a bootstrap resample is $(1 - 1/n)^n$, and find its limit.
3. Use the bootstrap to estimate the SE of the sample correlation for 30 points from a bivariate normal with $\rho = 0.6$, and compare with the approximation $(1 - \rho^2)/\sqrt{n}$.
4. Interview: You have 50 users with 20 sessions each and want the SE of the average session length. A colleague bootstraps the 1,000 sessions. What is wrong, and how would you fix it? Hint: sessions from the same user are correlated; resampling individual sessions treats them as independent and understates the SE; resample users (with all their sessions) instead, a cluster bootstrap.

Bootstrap Confidence Intervals

Why it matters. A standard error is half the story; decisions usually need an interval. The bootstrap distribution contains more information than its SD: its shape reveals skewness, which symmetric formula-based intervals ignore. Several ways to turn a bootstrap distribution into an interval exist, and they differ in accuracy, especially for skewed statistics and small samples. Knowing which to use, and that none is perfect, matters in practice, because bootstrap intervals are what many analysts use for medians, ratios, percentiles and machine learning metrics.

The main methods. Let $\hat\theta$ be the estimate, $\hat\theta^_b$ the bootstrap replicates, $\widehat{\operatorname{SE}}$ their SD, and $q^_\alpha$ their $\alpha$ quantile.

Method	95% interval	Idea	Pros and cons
Normal	$\hat\theta \pm 1.96,\widehat{\operatorname{SE}}$	assume $\hat\theta$ is normal, use the bootstrap only for the SE	simple; ignores skewness
Percentile	$[q^_{0.025},; q^_{0.975}]$	read the interval directly off the bootstrap distribution	intuitive, respects the range of $\theta$ (e.g. stays in $[0, 1]$), transformation-invariant; can be biased for skewed statistics
Basic (reverse percentile)	$[2\hat\theta - q^_{0.975},; 2\hat\theta - q^_{0.025}]$	the distribution of $\hat\theta^* - \hat\theta$ mimics that of $\hat\theta - \theta$	theoretically motivated; can leave the parameter space
BCa (bias-corrected and accelerated)	percentiles at adjusted levels	corrects the percentile method for bias ($z_0$) and for skewness, via how the SE changes with $\theta$ ($a$)	most accurate general-purpose choice; needs more computation (and a jackknife)
Studentized (bootstrap-$t$)	$[\hat\theta - t^_{0.975}\widehat{\operatorname{SE}},; \hat\theta - t^_{0.025}\widehat{\operatorname{SE}}]$	bootstrap the $t$-statistic $(\hat\theta^* - \hat\theta)/\widehat{\operatorname{SE}}^*$	very accurate when a good SE estimate exists for each resample; can be unstable

BCa in brief. The adjusted quantile levels are

$$\alpha_1 = \Phi\left(z_0 + \frac{z_0 + z_{\alpha/2}}{1 - a(z_0 + z_{\alpha/2})}\right), \qquad \alpha_2 = \Phi\left(z_0 + \frac{z_0 + z_{1-\alpha/2}}{1 - a(z_0 + z_{1-\alpha/2})}\right),$$

with $z_0 = \Phi^{-1}\big(\text{proportion of } \hat\theta^*_b < \hat\theta\big)$ measuring median bias, and the acceleration $a$ estimated from the skewness of jackknife values (next topic). When $z_0 = a = 0$, BCa reduces to the percentile interval.

How many resamples? About 1,000 suffice for a standard error, but interval endpoints depend on tail quantiles, so use at least 2,000 to 10,000 for 95% intervals, more for 99%.

import numpy as np
from statistics import NormalDist

rng = np.random.default_rng(821)
Phi, Phi_inv = NormalDist().cdf, NormalDist().inv_cdf

def boot_intervals(x, B=4_000, alpha=0.05):
    n, est = len(x), x.mean()
    reps = x[rng.integers(0, n, size=(B, n))].mean(axis=1)
    se = reps.std(ddof=1)
    z = Phi_inv(1 - alpha / 2)
    q_lo, q_hi = np.quantile(reps, [alpha / 2, 1 - alpha / 2])
    # BCa: bias correction z0 from the bootstrap distribution, acceleration a from the jackknife
    z0 = Phi_inv(np.clip(np.mean(reps < est), 1e-6, 1 - 1e-6))
    jack = (x.sum() - x) / (n - 1)                         # leave-one-out means
    d = jack.mean() - jack
    a = np.sum(d**3) / (6 * np.sum(d**2) ** 1.5)
    adj = lambda q: Phi(z0 + (z0 + Phi_inv(q)) / (1 - a * (z0 + Phi_inv(q))))
    b_lo, b_hi = np.quantile(reps, [adj(alpha / 2), adj(1 - alpha / 2)])
    return {
        "normal":     (est - z * se, est + z * se),
        "percentile": (q_lo, q_hi),
        "basic":      (2 * est - q_hi, 2 * est - q_lo),      # reflect the percentiles around the estimate
        "BCa":        (b_lo, b_hi),
    }

# One skewed sample: hospital stay lengths (days), lognormal
x = rng.lognormal(1.0, 1.0, size=30)
t_half = 2.045 * x.std(ddof=1) / np.sqrt(len(x))          # t_{0.975, 29}
print(f"sample mean {x.mean():.2f};  t-interval [{x.mean() - t_half:.2f}, {x.mean() + t_half:.2f}]")
for name, (lo, hi) in boot_intervals(x).items():
    print(f"{name:10s} [{lo:.2f}, {hi:.2f}]")

# Coverage over 2,000 repeated studies (true mean of lognormal(1, 1) = e^1.5)
mu = np.exp(1.5)
hits = {k: 0 for k in ["t", "normal", "percentile", "basic", "BCa"]}
reps = 2_000
for _ in range(reps):
    x = rng.lognormal(1.0, 1.0, size=30)
    h = 2.045 * x.std(ddof=1) / np.sqrt(30)
    hits["t"] += x.mean() - h <= mu <= x.mean() + h
    for name, (lo, hi) in boot_intervals(x, B=2_000).items():
        hits[name] += lo <= mu <= hi
print("\ncoverage of nominal 95% intervals, n = 30, lognormal data:")
print("  " + "   ".join(f"{k} {v / reps:.3f}" for k, v in hits.items()))

sample mean 5.95;  t-interval [3.10, 8.80]
normal     [3.28, 8.62]
percentile [3.57, 8.86]
basic      [3.04, 8.33]
BCa        [4.04, 9.89]

coverage of nominal 95% intervals, n = 30, lognormal data:
  t 0.892   normal 0.878   percentile 0.882   basic 0.860   BCa 0.901

On this sample, the intervals differ noticeably: BCa is shifted to the right compared to the others, because the sampling distribution of a mean of lognormal data is right-skewed and the true mean is more likely to be above the estimate than a symmetric interval suggests. The basic interval goes the wrong way, shifting left. The coverage simulation is sobering: with 30 strongly skewed observations, every method falls short of 95%. BCa comes closest at 90%, and the basic interval is worst at 86%. The bootstrap does not fix a small sample from a heavy-tailed population; it only works with the information the sample contains, and samples of 30 rarely contain the rare large values that pull the mean up.

Practical guidance.

• Default to BCa (for example scipy.stats.bootstrap(method="BCa")), or the percentile interval for a quick look.
• Check whether the bootstrap distribution is roughly symmetric and smooth. Strong discreteness (few distinct values) or extreme skewness is a warning sign.
• For a critical decision with an unusual statistic, check the method's coverage by simulation from a plausible model, as the code does.
• Bootstrap intervals are approximate even when everything is done right; their errors shrink as $n$ grows.

Key takeaways

• Percentile: read the 2.5% and 97.5% quantiles of the bootstrap replicates; normal: estimate $\pm 1.96$ bootstrap SE.
• BCa adjusts the percentile levels for bias and skewness and is the best general-purpose bootstrap interval.
• Use thousands of resamples for intervals, since tail quantiles need more Monte Carlo precision than SEs.
• With small samples from skewed populations, all bootstrap intervals under-cover; the method cannot exceed the information in the data.

🧪 Practice

1. A bootstrap distribution has quantiles $q^_{0.025} = 11.2$ and $q^_{0.975} = 15.8$, and the estimate is $\hat\theta = 12.9$. Compute the percentile and basic intervals by hand. Why is the basic interval shifted to the other side?
2. Show that the percentile interval is invariant under monotone transformations: if $[a, b]$ is the interval for $\theta$, then $[g(a), g(b)]$ is the interval for $g(\theta)$ for increasing $g$. Is the normal interval invariant?
3. In the code, change the sample size to $n = 100$ and then $n = 300$. How do the coverages change?
4. Interview: You need a 95% CI for the AUC of a classifier evaluated on a test set of 2,000 examples, of which 40 are positives. How would you bootstrap, and what would you watch out for? Hint: resample the test set (stratified by class, to keep 40 positives), recompute the AUC on each resample, use percentile or BCa; with only 40 positives the interval will be wide and discrete.

Jackknife

Why it matters. The jackknife (Quenouille, 1949, for bias; Tukey, 1958, for variance) is the bootstrap's older, deterministic cousin. Instead of random resamples, it uses the $n$ datasets obtained by leaving out one observation at a time. It is cheap, needs no random numbers, gives reproducible results, and is particularly good at estimating and removing bias. It is also a building block of other methods: the BCa interval uses it, and influence diagnostics in regression (Section 10.3) are leave-one-out ideas.

The leave-one-out replicates. For a statistic $\hat\theta = s(x_1, \dots, x_n)$, define

$$\hat\theta_{(i)} = s(x_1, \dots, x_{i-1}, x_{i+1}, \dots, x_n), \qquad \bar\theta_{(\cdot)} = \frac{1}{n}\sum_{i=1}^n \hat\theta_{(i)}.$$

Jackknife bias estimate and bias-corrected estimator.

$$\widehat{\operatorname{Bias}}{\text{jack}} = (n - 1)\left(\bar\theta{(\cdot)} - \hat\theta\right), \qquad \hat\theta_{\text{jack}} = \hat\theta - \widehat{\operatorname{Bias}}{\text{jack}} = n\hat\theta - (n - 1)\bar\theta{(\cdot)}.$$

Why the factor $n - 1$? If the bias has the common form $b/n$, then each leave-one-out estimate has bias $b/(n - 1)$, and the difference $\bar\theta_{(\cdot)} - \hat\theta$ has expectation $\frac{b}{n-1} - \frac{b}{n} = \frac{b}{n(n-1)}$. Multiplying by $n - 1$ recovers $b/n$. The correction removes bias of order $1/n$ exactly and leaves only bias of order $1/n^2$.

Jackknife standard error.

$$\widehat{\operatorname{SE}}{\text{jack}} = \sqrt{\frac{n - 1}{n}\sum{i=1}^n \left(\hat\theta_{(i)} - \bar\theta_{(\cdot)}\right)^2}.$$

The factor $(n - 1)/n$ (instead of the $1/(n-1)$ of an ordinary variance) compensates for the leave-one-out replicates being very similar to each other: each pair shares $n - 2$ observations.

Jackknife vs. bootstrap.

	Jackknife	Bootstrap
Resamples	$n$ deterministic leave-one-out sets	$B$ random with-replacement sets
Cost	$n$ evaluations	$B$ evaluations (usually 1,000+)
Reproducible	yes	only with a fixed seed
Bias correction	excellent for $O(1/n)$ bias	possible, noisier
Smooth statistics (means, ratios, correlations)	good SE	good SE
Non-smooth statistics (median, quantiles)	SE unreliable	works
Confidence intervals	only via normal approximation	full distribution: percentile, BCa

The jackknife is a linear approximation of the bootstrap; it is accurate when the statistic changes smoothly as observations are slightly reweighted.

import numpy as np

rng = np.random.default_rng(822)

def jackknife(x, stat):
    n = len(x)
    loo = np.array([stat(np.delete(x, i)) for i in range(n)])     # leave-one-out replicates
    est = stat(x)
    bias = (n - 1) * (loo.mean() - est)
    se = np.sqrt((n - 1) / n * np.sum((loo - loo.mean()) ** 2))
    return est, bias, est - bias, se

x = rng.normal(10, 3, size=15)

# 1. Plug-in variance (divisor n) is biased; the jackknife correction recovers S^2 exactly
est, bias, corrected, se = jackknife(x, lambda a: a.var())
print(f"plug-in variance {est:.4f}, jackknife bias {bias:+.4f}, corrected {corrected:.4f}, S^2 = {x.var(ddof=1):.4f}")

# 2. For the mean, the jackknife SE equals the textbook s / sqrt(n)
_, _, _, se = jackknife(x, np.mean)
print(f"jackknife SE of the mean {se:.4f}   s/sqrt(n) = {x.std(ddof=1) / np.sqrt(len(x)):.4f}")

# 3. A ratio estimator: revenue per visit = sum(revenue) / sum(visits) over 20 stores
visits = rng.poisson(400, 20).astype(float)
revenue = visits * rng.normal(2.5, 0.6, 20)
data = np.column_stack([revenue, visits])
ratio = lambda d: d[:, 0].sum() / d[:, 1].sum()
loo = np.array([ratio(np.delete(data, i, axis=0)) for i in range(20)])
se_jack = np.sqrt(19 / 20 * np.sum((loo - loo.mean()) ** 2))
boot = [ratio(data[rng.integers(0, 20, 20)]) for _ in range(5_000)]
print(f"ratio {ratio(data):.3f}: jackknife SE {se_jack:.4f}, bootstrap SE {np.std(boot):.4f}")

# 4. Failure: the median is not smooth. Leave-one-out medians take only a couple of distinct values,
#    so the jackknife SE is a very noisy estimate (compare with the bootstrap on the same samples)
x = rng.normal(10, 3, size=15)
print(f"distinct leave-one-out medians in one sample of 15: {np.unique([np.median(np.delete(x, i)) for i in range(15)]).size}")
true_se = np.median(rng.normal(10, 3, size=(50_000, 15)), axis=1).std()
jack_ses, boot_ses = [], []
for _ in range(2_000):
    x = rng.normal(10, 3, 15)
    jack_ses.append(jackknife(x, np.median)[3])
    boot_ses.append(np.median(x[rng.integers(0, 15, size=(500, 15))], axis=1).std())
print(f"median, true SE {true_se:.3f}:  jackknife SE mean {np.mean(jack_ses):.3f} (SD {np.std(jack_ses):.3f})   "
      f"bootstrap SE mean {np.mean(boot_ses):.3f} (SD {np.std(boot_ses):.3f})")

plug-in variance 4.1808, jackknife bias -0.2986, corrected 4.4794, S^2 = 4.4794
jackknife SE of the mean 0.5465   s/sqrt(n) = 0.5465
ratio 2.364: jackknife SE 0.1114, bootstrap SE 0.1081
distinct leave-one-out medians in one sample of 15: 3
median, true SE 0.961:  jackknife SE mean 0.923 (SD 0.585)   bootstrap SE mean 1.023 (SD 0.316)

The jackknife correction of the divisor-$n$ variance gives exactly $S^2$: the bias of the plug-in variance is $-\sigma^2/n$, precisely of the form the jackknife removes. For the mean, the jackknife SE reproduces $s/\sqrt n$ exactly. For a ratio estimator, where the SE formula would need the delta method, the jackknife and bootstrap agree closely. The median shows the limitation: with 15 observations, leaving one out produces only 3 distinct medians, so the jackknife SE rests on very little variation. Its average is about right, but it is almost twice as noisy from sample to sample as the bootstrap SE.

Key takeaways

• The jackknife recomputes the statistic $n$ times, leaving out one observation each time.
• Bias estimate: $(n - 1)(\bar\theta_{(\cdot)} - \hat\theta)$; it removes $O(1/n)$ bias exactly (it turns the divisor-$n$ variance into $S^2$).
• SE estimate: $\sqrt{\frac{n-1}{n}\sum(\hat\theta_{(i)} - \bar\theta_{(\cdot)})^2}$; it equals $s/\sqrt n$ for the mean.
• It is deterministic and cheap, good for smooth statistics, unreliable for non-smooth ones like the median.

🧪 Practice

1. Show algebraically that for $\hat\theta = \bar x$, the leave-one-out means are $\bar x_{(i)} = (n\bar x - x_i)/(n - 1)$, and that the jackknife SE equals $s/\sqrt n$.
2. For the data ${3, 7, 8, 12}$, compute the jackknife bias-corrected estimate of $\sigma^2$ starting from the divisor-$n$ variance.
3. Apply the jackknife to estimate the bias of $\hat\lambda = 1/\bar x$ for an exponential sample of size 10, and compare the corrected estimate's bias (by simulation) with the uncorrected one.
4. Interview: Training your model takes 10 minutes. You want an SE for its test-set metric that accounts for training randomness and data sampling. Would you use the bootstrap or a jackknife, and how? Hint: count model fits: $B = 1{,}000$ bootstrap refits take a week; a delete-$d$ or grouped jackknife over $k$ folds needs only $k$ refits; also consider bootstrapping only the test-set predictions if training variability is small.

Permutation Resampling

Why it matters. The bootstrap estimates the uncertainty of an effect. A different question is: could an effect this large have arisen by chance, if the groups were really no different? Permutation resampling answers it with the same "simulate instead of derive" philosophy. If the group labels are irrelevant, shuffling them should produce differences as large as the observed one fairly often; if the observed difference is extreme compared to the shuffled ones, the labels matter. Permutation methods need almost no assumptions (no normality, no large samples), and they are exact in a precise sense. Formal hypothesis testing is the subject of Chapter 9; this topic focuses on the resampling mechanism.

The setting. Two groups, $x_1, \dots, x_m$ and $y_1, \dots, y_n$, and a statistic such as $T = \bar x - \bar y$. Under the null hypothesis that the group label has no effect, all $N = m + n$ observations come from the same distribution, and the assignment of labels is arbitrary. More precisely, the observations are exchangeable: every rearrangement of them is equally likely.

The algorithm.

1. Compute the observed statistic $T_{\text{obs}}$.
2. Pool all $N$ observations. Randomly reassign $m$ of them to group 1 and the rest to group 2: a random permutation of the labels (resampling without replacement).
3. Compute $T^*$ for the relabeled data.
4. Repeat $B$ times. The values $T^_1, \dots, T^_B$ form the permutation distribution: the distribution of $T$ if the labels were meaningless.
5. The p-value is the proportion of permuted statistics at least as extreme as the observed one. For a Monte Carlo approximation use $p = \frac{1 + #{|T^*b| \ge |T{\text{obs}}|}}{1 + B}$, which counts the observed labeling as one of the permutations and never reports $p = 0$.

If $N$ is small, all $\binom{N}{m}$ relabelings can be enumerated, giving an exact p-value.

Why it is exact. If the labels were assigned at random (a randomized experiment) and the treatment has no effect, the observed labeling is just one random draw from all possible labelings. The permutation distribution is then exactly the distribution of $T$ under the null hypothesis, whatever the shape of the data's distribution and whatever the sample size. This is Fisher's randomization argument (1935), and it links the design of experiments (Section 7.2) directly to their analysis.

Bootstrap vs. permutation.

	Bootstrap	Permutation
Resampling	with replacement, within each group	without replacement, across groups (shuffle labels)
World simulated	the one that produced the data (effect included)	the one with no effect (null hypothesis)
Distribution centered at	the estimate $\hat\theta$	the null value (0 for a difference)
Main output	SE, bias, confidence interval	p-value
Key assumption	sample represents the population	exchangeability under the null

import numpy as np
from itertools import combinations

rng = np.random.default_rng(823)

# Task completion times (s) with the old and new UI design
old = np.array([41.2, 38.5, 45.1, 50.3, 39.8, 47.6, 44.0, 52.1, 36.9, 43.3])
new = np.array([35.4, 40.1, 33.8, 37.2, 42.5, 31.9, 36.6, 39.0])
observed = old.mean() - new.mean()
pooled = np.concatenate([old, new])
n_old = len(old)

# Monte Carlo permutation distribution: if the design has no effect, the labels are arbitrary
B = 20_000
perm_diffs = np.empty(B)
for b in range(B):
    shuffled = rng.permutation(pooled)                  # resample WITHOUT replacement: a relabeling
    perm_diffs[b] = shuffled[:n_old].mean() - shuffled[n_old:].mean()
p_mc = (np.sum(np.abs(perm_diffs) >= abs(observed)) + 1) / (B + 1)
print(f"observed difference {observed:.2f} s; permutation SD {perm_diffs.std():.2f}; Monte Carlo p = {p_mc:.4f}")

# Exact version: enumerate all C(18, 10) = 43,758 ways to choose the 'old' group
total = pooled.sum()
diffs = np.array([pooled[list(c)].sum() / n_old - (total - pooled[list(c)].sum()) / len(new)
                  for c in combinations(range(len(pooled)), n_old)])
print(f"exact: {len(diffs):,} relabelings, p = {np.mean(np.abs(diffs) >= abs(observed) - 1e-12):.4f}")

# Contrast: the bootstrap resamples each group WITH replacement to estimate the effect's uncertainty
boot = [rng.choice(old, n_old).mean() - rng.choice(new, len(new)).mean() for _ in range(B)]
print(f"bootstrap 95% percentile CI for the difference: [{np.quantile(boot, 0.025):.2f}, {np.quantile(boot, 0.975):.2f}]")
print(f"centers: permutation distribution {perm_diffs.mean():+.2f}, bootstrap distribution {np.mean(boot):+.2f}")

observed difference 6.82 s; permutation SD 2.61; Monte Carlo p = 0.0045
exact: 43,758 relabelings, p = 0.0052
bootstrap 95% percentile CI for the difference: [3.13, 10.50]
centers: permutation distribution +0.00, bootstrap distribution +6.80

The new design is 6.82 s faster on average. Among 20,000 random relabelings, a difference this large (in either direction) occurred only about 0.45% of the time, and full enumeration of all 43,758 possible labelings gives the exact two-sided p-value 0.0052: if the design had no effect, a difference this large would be very unusual. The bootstrap answers a different question with a different resampling scheme: its distribution is centered at the observed effect (+6.80), not at 0, and gives a 95% interval of about 3 to 10.5 seconds for the size of the improvement. In a report, the two are complementary: the permutation p-value says the effect is unlikely to be chance, and the bootstrap interval says how big it plausibly is.

Beyond two means. Any statistic can be used: a difference of medians, a ratio of variances, a correlation (permute one variable relative to the other), or a classifier's accuracy (permute the labels to see what accuracy chance alone gives). For paired data, the permutation is a random sign flip of each within-pair difference. The null hypothesis must justify the exchange: for comparing means of groups with very different variances, the groups are not exchangeable, and the permutation test can be misleading.

Key takeaways

• Permutation resampling shuffles group labels (without replacement) to simulate the world where the labels do not matter.
• The p-value is the share of permuted statistics at least as extreme as the observed one; use $(1 + \text{count})/(1 + B)$ for Monte Carlo approximations.
• It is exact under exchangeability, for any distribution and any sample size; in randomized experiments, exchangeability comes from the randomization itself.
• Bootstrap: uncertainty of an estimate (centered at the estimate). Permutation: null distribution of a statistic (centered at no effect).

🧪 Practice

1. Groups ${3, 5}$ and ${8, 9, 10}$: enumerate all $\binom{5}{2} = 10$ relabelings, compute the difference in means for each, and find the exact two-sided p-value.
2. Describe how to run a permutation test for the correlation between hours studied and exam score in a sample of 25 students. What exactly is permuted?
3. Explain why a permutation p-value can never be smaller than $1/\binom{N}{m}$, and compute this minimum for two groups of 4.
4. Interview: Your fraud model reaches 92% accuracy on a dataset where 90% of transactions are legitimate. How could a permutation approach tell you whether the model has learned anything? Hint: shuffle the labels, retrain and evaluate many times to get the accuracy distribution under "no relationship"; compare 92% with it; also note that 90% accuracy is reachable by always predicting "legitimate", so choose a better metric too.

Cross-Validation

Why it matters. A model's error on the data it was fitted to is optimistic, and the more flexible the model, the more optimistic it is: a model can memorize its training data and still predict new cases poorly. What we want is the prediction error on new data, a parameter of the data-generating process, and it must be estimated. Holding out a separate test set is the simplest estimator, but it wastes data and depends on one random split. Cross-validation reuses the data so that every observation serves for both fitting and evaluation, never at the same time. It is the standard tool for choosing among models and tuning hyperparameters (Sections 10.4 and 12.5).

The quantity being estimated. For a model fitted by some procedure on $n$ observations, the expected prediction error is $\operatorname{Err} = E\big[L(Y_{\text{new}}, \hat f(X_{\text{new}}))\big]$, for example the expected squared error on a new observation from the same population. The training error $\frac{1}{n}\sum L(y_i, \hat f(x_i))$ underestimates it, because $\hat f$ was chosen to fit exactly these points.

$k$-fold cross-validation.

1. Randomly split the data into $k$ roughly equal folds.
2. For each fold $j = 1, \dots, k$: fit the model on the other $k - 1$ folds, and compute its error on fold $j$.
3. Average the $k$ fold errors: $\operatorname{CV}k = \frac{1}{k}\sum{j=1}^k \operatorname{Err}_j$.

  fold:        1        2        3        4        5
  round 1:  [ test ][ train ][ train ][ train ][ train ]   -> Err_1
  round 2:  [ train][  test ][ train ][ train ][ train ]   -> Err_2
  round 3:  [ train][ train ][  test ][ train ][ train ]   -> Err_3
  round 4:  [ train][ train ][ train ][  test ][ train ]   -> Err_4
  round 5:  [ train][ train ][ train ][ train ][  test ]   -> Err_5
                                                 CV error = average of Err_1 ... Err_5

Choosing $k$. The choice is itself a bias-variance trade-off for the error estimate:

Scheme	Training size per fit	Bias of the error estimate	Variance	Cost
Holdout (single split, e.g. 70/30)	$0.7n$	pessimistic	high (one split)	1 fit
5-fold or 10-fold	$0.8n$ to $0.9n$	slightly pessimistic	moderate	5 to 10 fits
Leave-one-out (LOOCV, $k = n$)	$n - 1$	nearly unbiased	can be high	$n$ fits
Repeated $k$-fold	as $k$-fold	as $k$-fold	lower	$r \times k$ fits

Each fit uses fewer than $n$ observations, so CV estimates the error of a model trained on slightly less data, and it is slightly pessimistic. 5-fold and 10-fold CV are the usual compromise.

import numpy as np

rng = np.random.default_rng(824)

# True relationship: a smooth curve plus noise; we fit polynomials of increasing degree
f = lambda x: np.sin(2 * x) + 0.5 * x
n = 40
x = rng.uniform(0, 3, n)
y = f(x) + rng.normal(0, 0.4, n)
x_test = rng.uniform(0, 3, 100_000)                       # a huge fresh sample: the "true" error
y_test = f(x_test) + rng.normal(0, 0.4, x_test.size)

def mse(deg, x_tr, y_tr, x_ev, y_ev):
    coef = np.polyfit(x_tr, y_tr, deg)
    return np.mean((np.polyval(coef, x_ev) - y_ev) ** 2)

def kfold_cv(deg, k=5):
    folds = np.array_split(rng.permutation(n), k)          # shuffle, then split into k folds
    errs = []
    for i in range(k):
        val = folds[i]
        tr = np.concatenate([folds[j] for j in range(k) if j != i])
        errs.append(mse(deg, x[tr], y[tr], x[val], y[val]))  # fit on k-1 folds, evaluate on the held-out fold
    return np.mean(errs)

def loocv(deg):
    return np.mean([mse(deg, np.delete(x, i), np.delete(y, i), x[i:i + 1], y[i:i + 1]) for i in range(n)])

print(f"{'degree':>6s} {'training':>9s} {'5-fold CV':>10s} {'LOOCV':>8s} {'true test':>10s}")
for deg in [1, 2, 3, 4, 5, 7, 10]:
    print(f"{deg:6d} {mse(deg, x, y, x, y):9.3f} {kfold_cv(deg):10.3f} {loocv(deg):8.3f} "
          f"{mse(deg, x, y, x_test, y_test):10.3f}")
print("irreducible error (noise variance): 0.160")

degree  training  5-fold CV    LOOCV  true test
     1     0.318      0.378    0.358      0.341
     2     0.290      0.364    0.355      0.345
     3     0.165      0.221    0.223      0.190
     4     0.141      0.220    0.213      0.210
     5     0.137      0.221    0.228      0.209
     7     0.135      0.259    0.277      0.216
    10     0.127      1.653    0.424      0.221
irreducible error (noise variance): 0.160

The training error falls steadily as the degree increases, so it would always choose the most complex model. The true test error, computed on 100,000 fresh points, is lowest at degree 3 and rises slowly after that. Both cross-validation estimates track the true error and select degree 4, with degree 3 essentially tied: close to the best model. They are somewhat pessimistic, as expected, since each fit uses fewer points. For degree 10, 5-fold CV reports a huge error: fitted on 32 points, a degree-10 polynomial swings wildly near the edges of the data. This exaggerates the error of the model fitted on all 40 points, but it correctly flags the high-degree model as unstable.

Pitfalls.

• Leakage. Every step that learns from data (scaling, imputation, feature selection, oversampling) must be fitted inside each training fold only. Selecting features on the full dataset and then cross-validating gives wildly optimistic estimates.
• Dependent observations. If the same patient, user or device appears in several rows, put all of its rows in the same fold (grouped $k$-fold). Otherwise the model is evaluated on near-copies of its training data.
• Time series. Never use the future to predict the past: use forward-chaining (rolling-origin) validation, training on data up to time $t$ and testing after $t$ (Section 12.1).
• Selection optimism. The CV error of the model chosen as best among many is itself optimistic, because it was selected for being low. To report an honest error for the whole procedure, use a separate test set or nested cross-validation (an inner loop for tuning, an outer loop for evaluation).
• Class imbalance. Use stratified folds so that each fold has the same class proportions.

Key takeaways

• Training error underestimates prediction error, and more so for flexible models.
• $k$-fold CV fits on $k - 1$ folds and evaluates on the remaining one, rotating through all folds; 5 or 10 folds are standard.
• CV is slightly pessimistic (smaller training sets) and has its own variance; repeated CV reduces the variance.
• Avoid leakage: all data-dependent preprocessing goes inside the folds; respect groups and time order; use nested CV to evaluate a tuned model.

🧪 Practice

1. With 1,000 observations and 10-fold CV, how many observations are used to fit each model and how many to evaluate it? How many model fits does LOOCV need?
2. Modify the code to run 5-fold CV 20 times with different random splits for each degree, and report the mean and SD of the CV error. Which degrees are clearly distinguishable?
3. Explain why selecting the 20 features most correlated with the outcome on the full dataset, then running 10-fold CV, gives an optimistic estimate. How would you fix it?
4. Interview: Your model's 5-fold CV accuracy is 94%, but in production it is 81%. List possible causes. Hint: leakage (preprocessing or features using future information), duplicate or grouped records across folds, a shift between the training and production populations, selection optimism from heavy tuning, and time-ordered data validated with random folds.

<a id="9-hypothesis-testing"></a>