Chapter 8
Estimation

Lectures 19–22 · from data to numbers we can trust: estimators, bias and variance, the method of moments, maximum likelihood and confidence intervals

an estimate is a rule for turning data into a guess, and some rules are better than others

In the Second World War, the Allies badly wanted to know how many tanks Germany was building each month. Intelligence reports, based on spies, aerial photographs and interrogations, put the figure at well over a thousand. Statisticians tried something different. Captured German tanks carried serial numbers on their gearboxes and other parts, and the numbers ran in sequence. If you capture a handful of tanks at random and the largest serial number you see is 61, the enemy probably hasn't built thousands. After the war, German production records were compared with both sets of figures. The intelligence estimates had been wildly too high. The statistical estimates were remarkably close.

That is estimation: using data to put a number on something we can't see directly. Chapters 3 to 7 assumed that we knew the model's parameters: the rate $λ$ , the probability $p$ , the mean $μ$ . In practice we almost never do. This chapter turns the direction round. Given data, what values of the parameters are reasonable, and how sure can we be?

Estimation is the heart of inference. A point estimate is our best single guess at an unseen quantity: the true detection rate of a test, the true average speed of light, the true effect of a drug. A confidence interval adds an honest statement of how far off that guess might be. Estimation also serves prediction, because every prediction from a model uses estimated parameters, and the uncertainty in those estimates carries into the prediction. Section 8.4 shows the difference between an interval for a mean and an interval for a new observation, which is much wider.

Here is the route through the chapter.

Section 8.1 describes what makes an estimator good: bias, standard error and mean squared error.
Sections 8.2 and 8.3 give two general methods for building estimators: the method of moments and maximum likelihood.
Section 8.4 builds confidence intervals for a mean, and Section 8.5 for a proportion.
Section 8.6 explains what "95% confidence" means, and what it doesn't.

Think About It: Before You Start

Five captured tanks carry the serial numbers 19, 56, 24, 16 and 61. Tanks are numbered $1, 2, \dots, N$ . Write down your best guess for $N$ , and the rule you used to get it. Is your rule more likely to guess too high or too low?

Many people say 61, the largest number seen. But the largest serial number in a sample can never be bigger than $N$ , and is usually smaller, so this rule always guesses too low (or exactly right). Others double the average: the average of 1 to $N$ is about $N / 2$ , and the sample average is 35.2, giving about 70 (or 69.4 after a small correction, $2 \overline{x} - 1$ ). A third rule adds to the largest serial number the typical gap between the numbers seen: 61 plus the average gap, $(61$ $- 5) / 5$ $\approx 11$ , giving about 72. Section 8.1 compares the three rules, and this one turns out to be the best.

Lecture 19 — Estimators and the method of moments

8.1 Estimators and their errors

Parameters, estimators and estimates

A parameter is a number that describes the population or the model: $μ$ , $σ$ , $p$ , $λ$ , $N$ . We use the Greek letter $θ$ (theta) for a generic parameter.

Definition: Estimator and Estimate

An estimator $\hat{θ}$ of a parameter $θ$ is a statistic, a rule $\hat{θ}$ $= g (X_{1},$ $\dots,$ $X_{n})$ for turning a sample into a guess. The number it produces from a particular dataset is an estimate, or point estimate.

An estimator is a random variable: it has a sampling distribution (Chapter 7). That's what lets us say how good it is. A single estimate can't be judged on its own, since we don't know $θ$ ; but a rule can be judged by how it behaves over many possible samples.

Definition: Bias, Standard Error and Mean Squared Error

For an estimator $\hat{θ}$ of $θ$ :

its bias is $Bias [\hat{θ}]$ $= E [\hat{θ}]$ $- θ$ ; the estimator is unbiased if the bias is 0;
its standard error is $SE [\hat{θ}] = SD [\hat{θ}]$ ;
its mean squared error is $MSE [\hat{θ}]$ $= E [{(\hat{θ} - θ)}^{2}]$ , and its square root, the root mean squared error (RMSE), is the typical size of its error.

The picture to keep in mind is a dartboard with the true $θ$ at the bullseye and each estimate a dart. Bias is aiming at the wrong place; standard error is a shaky hand. Chapter 1's bathroom scale that always reads 2 kg light (Section 1.3) has large bias and no variance at all: precisely wrong.

Bias is aiming at the wrong place (the average of the darts is off-centre); variance is a shaky hand (the darts are spread out). Only the top-left estimator is both accurate and precise.

Throw some darts yourself, and watch how the two sources of error add up:

Theorem: MSE = Variance + Bias²

$MSE [\hat{θ}] = Var [\hat{θ}] + (Bias [\hat{θ}])^{2} .$

Write $m = E [\hat{θ}]$ and add and subtract it: $\begin{matrix} E [{(\hat{θ} - θ)}^{2}] & = E [{(\hat{θ} - m)}^{2}] + 2 (m - θ) E [\hat{θ} - m] + {(m - θ)}^{2} \\ = Var [\hat{θ}] + 0 + Bias {[\hat{θ}]}^{2} . ◼ \end{matrix}$ $E [{(\hat{θ} - θ)}^{2}]$ $= E [{(\hat{θ} - m)}^{2}]$ $+ 2 (m$ $- θ) E [\hat{θ}$ $- m]$ $+ {(m - θ)}^{2}$ $= Var [\hat{θ}]$ $+ 0$ $+ Bias {[\hat{θ}]}^{2} . ◼$ ∎

So the total error has two sources, and a little bias can be worth accepting if it buys a large cut in variance. This bias–variance trade-off runs through all of statistics and machine learning.

Example 8.1: How Many Tanks? ★ lecture example

Tanks are numbered $1, \dots, N$ , and we capture $k$ at random (without replacement). Let $M$ be the largest serial number among them. Compare three estimators of $N$ by simulation, with $N = 300$ and $k = 5$ :

${\hat{N}}_{1} = 2 \overline{X} - 1$ (twice the sample mean, minus 1);
${\hat{N}}_{2} = M$ (the largest serial number);
${\hat{N}}_{3}$ $= M (1 + \frac{1}{k})$ $- 1$ (the largest, plus the average gap between serial numbers).

Then apply all three to the five tanks in "before you start".

Why these rules? The serial numbers are equally likely to be any of $1, \dots, N$ , whose mean is $\frac{N + 1}{2}$ ; setting $\overline{X} = \frac{N + 1}{2}$ gives ${\hat{N}}_{1}$ . The largest number seen, $M$ , can't exceed $N$ . And the $k$ numbers seen split the range $1$ to $M$ into gaps averaging about $M / k$ , so we expect a gap of about that size above $M$ as well, which gives ${\hat{N}}_{3}$ .

Simulate 20,000 sets of 5 captured tanks:

set.seed(8)
N <- 300; k <- 5
sims <- replicate(20000, { s <- sample(1:N, k)
  c(twice_mean = 2 * mean(s) - 1, largest = max(s), largest_plus_gap = max(s) * (1 + 1/k) - 1) })
round(rbind(average = rowMeans(sims), SE = apply(sims, 1, sd),
            RMSE = sqrt(rowMeans((sims - N)^2))), 1)

output

        twice_mean largest largest_plus_gap
average      300.4   251.3            300.5
SE            76.7    41.7             50.1
RMSE          76.7    64.1             50.1

d <- data.frame(estimate = as.vector(t(sims)),
                rule = rep(c("twice the mean", "largest", "largest + gap"), each = 20000))
d$rule <- factor(d$rule, levels = c("twice the mean", "largest", "largest + gap"))
ggplot(d, aes(estimate)) + geom_histogram(binwidth = 15, boundary = 0, colour = "white") +
  geom_vline(xintercept = N, linetype = "dashed") + facet_wrap(~ rule) +
  scale_x_continuous(breaks = c(150, 300, 450)) +
  labs(x = "estimate of N (true N = 300)", y = "count")

rule	average	SE	RMSE
twice the mean	300.4	76.7	76.7
largest	251.3	41.7	64.1
largest + gap	300.5	50.1	50.1

Twice the mean is unbiased but very variable: one unusually large or small serial number moves it a lot.
The largest is much less variable, but biased: on average it's 49 tanks too low, and it can never be too high.
Largest plus gap is unbiased and less variable than twice the mean, so it has the smallest typical error. It is the best unbiased estimator possible here: no unbiased estimator has a smaller standard error (a fact we won't prove).

The five captured tanks (19, 56, 24, 16, 61) give ${\hat{N}}_{1}$ $= 2$ $\times 35.2$ $- 1$ $= 69.4$ , ${\hat{N}}_{2} = 61$ and ${\hat{N}}_{3}$ $= 61$ $\times 1.2$ $- 1$ $= 72.2$ . We'd estimate about 72 tanks.

The German tank problem: three ways to guess N

The enemy has numbered its tanks 1, 2, …, N. We capture k of them at random and read their serial numbers. Each press of capture simulates many such samples and applies three estimators: twice the mean (method of moments), the largest serial number (maximum likelihood), and the largest plus the average gap, m(1 + 1/k) − 1. The dashed line is the true N. Which estimator would you want your army to use?

true number of tanks N = tanks captured k =

History: Counting Tanks from Serial Numbers

Allied economists in the Second World War read the serial numbers on captured German tanks, including those stamped on gearboxes and on the moulds of road wheels, and estimated monthly production. After the war, the figures were compared with German records. In one widely quoted comparison from Richard Ruggles and Henry Brodie's 1947 account, for June 1940 the statistical estimate was 169 tanks a month, intelligence estimates said about 1,000, and German records showed 122; for August 1942 the figures were 327, about 1,550, and 342. Similar methods were later used to estimate the production of other equipment, and the idea is still used, for example to estimate how many units of a product a company has sold from the serial numbers of a few.

Source: R. Ruggles and H. Brodie (1947), "An empirical approach to economic intelligence in World War II", Journal of the American Statistical Association 42, 72–91.

Why the sample variance divides by n − 1

Chapter 2 promised to explain why the sample variance divides by $n - 1$ rather than $n$ . Here is the answer.

Theorem: The Sample Variance Is Unbiased

If $X_{1}, \dots, X_{n}$ is a random sample with variance $σ^{2}$ , then $S^{2}$ $= \frac{1}{n - 1} \sum_{i = 1}^{n} {(X_{i} - \overline{X})}^{2}$ satisfies $E [S^{2}] = σ^{2}$ .

Expand the square: $\sum {(X_{i} - \overline{X})}^{2}$ $= \sum X_{i}^{2}$ $- n {\overline{X}}^{2}$ . For any random variable, $E [Y^{2}]$ $= Var [Y]$ $+ E {[Y]}^{2}$ , so $E [X_{i}^{2}]$ $= σ^{2}$ $+ μ^{2}$ and $E [{\overline{X}}^{2}]$ $= \frac{σ^{2}}{n}$ $+ μ^{2}$ (Section 7.2). Therefore $\begin{matrix} E [\sum {(X_{i} - \overline{X})}^{2}] & = n (σ^{2} + μ^{2}) - n (\frac{σ^{2}}{n} + μ^{2}) \\ = (n - 1) σ^{2}, \end{matrix}$ $E [\sum {(X_{i} - \overline{X})}^{2}]$ $= n (σ^{2}$ $+ μ^{2})$ $- n (\frac{σ^{2}}{n} + μ^{2})$ $= (n$ $- 1) σ^{2},$ and dividing by $n - 1$ gives exactly $σ^{2}$ . $◼$ ∎

The deviations are measured from $\overline{X}$ , which sits closer to the data than the true $μ$ does, so they're slightly too small on average; dividing by $n - 1$ instead of $n$ compensates exactly.

Example 8.2: Is Unbiased Always Better?

For samples of size 5 from $N (0, 1)$ (so $σ^{2} = 1$ ), compare by simulation the unbiased $S^{2}$ with ${\hat{σ}}^{2} = \frac{n - 1}{n} S^{2}$ , which divides by $n$ .

set.seed(8)
v <- replicate(20000, { x <- rnorm(5); c(divide_by_4 = var(x), divide_by_5 = var(x) * 4 / 5) })
rbind(average = rowMeans(v), MSE = rowMeans((v - 1)^2))

output

        divide_by_4 divide_by_5
average   1.0024569   0.8019656
MSE       0.5003244   0.3594214

$S^{2}$ averages 1.00 (unbiased), and ${\hat{σ}}^{2}$ averages 0.80 (biased low by 20%). But the MSEs are 0.50 and 0.36: the biased estimator has the smaller typical error, because shrinking every estimate by $\frac{4}{5}$ cuts the variance more than it adds in squared bias. "Unbiased" is a good property, not the only one. We still use $S^{2}$ by convention, and because the $t$ intervals of Section 8.4 are built on it.

Comparing estimators: bias and spread

Each estimator is a rule for turning a sample into a guess. Draw many samples from a population whose true value we know (the vertical line), apply every rule to each sample, and compare. A good estimator is centred on the truth (unbiased) and tightly clustered (small standard error), like darts grouped on the bullseye.

estimate sample size n =

One more property is worth naming. An estimator is consistent if it gets arbitrarily close to $θ$ with probability approaching 1 as $n \to \infty$ . The Law of Large Numbers (Section 7.3) says that $\overline{X}$ is consistent for $μ$ . Consistency is the least we should ask: with enough data, the estimator should find the truth.

Check Yourself

An estimator has bias 2 and standard error 3. What is its mean squared error?

MSE =

Estimator A is unbiased with SE 10. Estimator B has bias 3 and SE 5. Which has the smaller typical error? A, because it is unbiased.B: its MSE is 34, against A's 100.They are equally good.

Think About It: A Biased Scale

A lab's scale reads, on average, 0.3 mg too heavy, with a measurement SD of 0.1 mg. A second scale is correctly calibrated but has SD 0.5 mg. A forensic chemist must weigh a seized sample once. Which scale would you use, and why? What would you do if you could weigh it 100 times?

For one weighing, the first scale's RMSE is $\sqrt{{0.1}^{2} + {0.3}^{2}}$ $\approx 0.32$ mg and the second's is 0.5 mg, so the biased scale is better. But averaging 100 weighings shrinks the standard error by a factor of 10 and leaves bias untouched: the first scale's average would still be about 0.3 mg too heavy, with RMSE about 0.30, while the second's would have RMSE 0.05. Repetition fixes variance, not bias. In court, a known bias should be corrected (by recalibrating or subtracting it) and disclosed, since a systematic error affects every case the lab handles in the same direction.

Key Takeaways

An estimator is a rule; its bias, standard error and MSE describe how it behaves over repeated samples.
$MSE$ $= Var$ $+ {Bias}^{2}$ : total error has two sources, and they can be traded off.
$S^{2}$ (dividing by $n - 1$ ) is unbiased for $σ^{2}$ . Unbiased isn't always best.
More data reduces variance, not bias.

8.2 The method of moments

The oldest general recipe for estimators is to make the model match the data's averages. Karl Pearson developed it in the 1890s.

Definition: Method of Moments

The $k$ th population moment is $E [X^{k}]$ and the $k$ th sample moment is $\frac{1}{n} \sum_{i = 1}^{n} X_{i}^{k}$ . For a model with $m$ unknown parameters, the method of moments sets the first $m$ population moments (which depend on the parameters) equal to the first $m$ sample moments, and solves for the parameters.

Equivalently, for two parameters, match the mean and the variance, using the variance with divisor $n$ .

Example 8.3: Suspicious Logins ★ lecture example

A security analyst models the time between suspicious login attempts as $Exp (λ)$ , and records five gaps, in minutes: 10.2, 12.6, 7.5, 8.3, 10.5. Estimate $λ$ by the method of moments.

One parameter, so match the first moment. The population mean of $Exp (λ)$ is $\frac{1}{λ}$ (Section 5.4), and the sample mean is $\overline{x} = 9.82$ minutes. Setting $\frac{1}{\hat{λ}} = 9.82$ gives $\hat{λ}$ $= \frac{1}{\overline{x}}$ $= \frac{1}{9.82}$ $\approx 0.102 per minute,$ about one suspicious attempt every 10 minutes.

Critical thinking. Is an exponential model plausible? For the exponential, the SD equals the mean. These five gaps have SD only 2.0 minutes, a fifth of their mean, so they're far more regular than exponential gaps would be. That suggests automated attempts at a near-steady rate, not random human ones. Estimating a parameter doesn't check the model; Chapter 9 gives tools that do.

Example 8.4: A Model for Rainfall ★ lecture example

R's precip gives the average annual rainfall, in inches, of 70 US cities. Rainfall is positive and right-skewed, so fit a gamma model, $Gamma (α, λ)$ , by the method of moments. Does it fit?

Two parameters, two moments. A $Gamma (α, λ)$ has mean $\frac{α}{λ}$ and variance $\frac{α}{λ^{2}}$ (Section 5.4). Matching them to the sample mean $\overline{x}$ and the divide-by- $n$ variance $v$ : $\frac{α}{λ} = \overline{x}, \frac{α}{λ^{2}} = v .$ Dividing the first equation by the second gives $λ = \overline{x} / v$ , and then $α = \overline{x} λ = {\overline{x}}^{2} / v$ .

m <- mean(precip); v <- mean((precip - m)^2)
c(n = length(precip), mean = m, sd = sqrt(v))
c(alpha = m^2 / v, lambda = m / v)

output

       n     mean       sd 
70.00000 34.88571 13.60839 
    alpha    lambda 
6.5717576 0.1883796

With $\overline{x} = 34.9$ inches and $\sqrt{v} = 13.6$ inches, the estimates are $\hat{α} = 6.57$ and $\hat{λ} = 0.188$ per inch.

ggplot(data.frame(precip), aes(precip)) +
  geom_histogram(aes(y = after_stat(density)), binwidth = 5, boundary = 0, colour = "white") +
  stat_function(fun = dgamma, args = list(shape = m^2 / v, rate = m / v), linewidth = 1) +
  labs(x = "average annual rainfall (inches)", y = "density")

The fit is fair but not great. The data have a sharper peak around 35 to 45 inches and a handful of very dry desert cities in the left tail, which the smooth gamma curve can't capture: the cities come from two different climates. As in Chapter 5, a model that matches the mean and variance can still miss the shape.

The method of moments is quick and needs no calculus, but it doesn't always give sensible answers (in Example 8.1 it can even produce an estimate smaller than the largest serial number seen, which is impossible), and it usually isn't the most precise method. The next section's method usually is.

Check Yourself

A server receives 4, 3, 2, 5, 4, 6, 3, 2, 5, 3, 4, 3 packets in 12 consecutive one-second intervals. Fit a Poisson model by the method of moments.

λ̂ =

Key Takeaways

Method of moments: match the first few population moments to the sample moments and solve.
Easy and general, but not always sensible or efficient; check the fitted model against the data.

Lecture 20 — Maximum likelihood

8.3 Maximum likelihood

The likelihood

Here is a more principled idea. Different parameter values make the observed data more or less probable. Choose the value that makes the data we actually saw most probable.

Definition: Likelihood and Maximum Likelihood Estimator

Let $x_{1}, \dots, x_{n}$ be observed data from a model with PMF or density $f (x; θ)$ . The likelihood function is the probability (or density) of the observed data, viewed as a function of $θ$ : $L (θ) = f (x_{1}; θ) f (x_{2}; θ) \dots f (x_{n}; θ)$ for an IID sample. The maximum likelihood estimate (MLE) $\hat{θ}$ is the value of $θ$ that maximizes $L (θ)$ . It is usually easier to maximize the log-likelihood $ℓ (θ)$ $= ln L (θ)$ $= \sum_{i} ln f (x_{i}; θ)$ , which has its maximum at the same place because $ln$ is increasing.

The likelihood is not a probability distribution for $θ$ (it needn't integrate to 1 over $θ$ ). It's a measure of how well each parameter value explains the data. The recipe is the calculus of Section 0.4: differentiate $ℓ$ , set the derivative to zero, solve, and check that it's a maximum.

Example 8.5: Who Uses Two-Factor Authentication? ★ lecture example

In a random sample of 10 people, 3 use two-factor authentication (2FA) on their email: the 1st, 3rd and 10th. Model each person as $Bernoulli (p)$ . Find the MLE of $p$ .

Likelihood. Each user contributes a factor $p$ and each non-user $1 - p$ : $L (p)$ $= p$ $\cdot (1$ $- p)$ $\cdot p$ $\cdot {(1 - p)}^{6}$ $\cdot p$ $= p^{3} {(1 - p)}^{7} .$

Log-likelihood and its derivative: $ℓ (p)$ $= 3 ln p$ $+ 7 ln (1$ $- p),$ $ℓ' (p)$ $= \frac{3}{p}$ $- \frac{7}{1 - p} .$

Solve $ℓ' (\hat{p}) = 0$ : $3 (1 - \hat{p}) = 7 \hat{p}$ , so $\hat{p} = \frac{3}{10}$ . The second derivative, $ℓ ″ (p)$ $= - \frac{3}{p^{2}}$ $- \frac{7}{{(1 - p)}^{2}}$ , is negative, so this is a maximum.

p <- seq(0.02, 0.9, by = 0.005)
ggplot(data.frame(p, ll = 3 * log(p) + 7 * log(1 - p)), aes(p, ll)) +
  geom_line(linewidth = 1) + geom_vline(xintercept = 0.3, linetype = "dashed") +
  labs(x = "p", y = "log-likelihood")

In general, $x$ successes in $n$ Bernoulli trials give $\hat{p} = \frac{x}{n}$ : the sample proportion, as in Section 0.4's preview. Only the number of successes matters, not which people they were; the same answer comes from the binomial likelihood $(\binom{10}{3}) p^{3} {(1 - p)}^{7}$ , since the constant doesn't move the peak.

Example 8.6: Time to a Toxic Effect ★ lecture example

Researchers record the time, in minutes, until a compound produces a toxic effect in five test subjects: 36.2, 42.8, 47.2, 50.2, 51.7. Assuming an exponential model, find the MLE of $λ$ .

Likelihood. Each time $x_{i}$ contributes $λ e^{- λ x_{i}}$ : $L (λ)$ $= \prod_{i = 1}^{5} λ e^{- λ x_{i}}$ $= λ^{5} e^{- λ \sum x_{i}}$ $= λ^{5} e^{- 228.1 λ} .$

Log-likelihood: $ℓ (λ)$ $= 5 ln λ$ $- 228.1 λ$ , so $ℓ' (λ)$ $= \frac{5}{λ}$ $- 228.1$ , which is zero at $\hat{λ} = \frac{5}{228.1} \approx 0.0219 per minute,$ and $ℓ ″ (λ)$ $= - 5 / λ^{2}$ $< 0$ confirms a maximum. In general, $\hat{λ}$ $= n / \sum x_{i}$ $= 1 / \overline{x}$ , the same as the method of moments here. The estimated mean time to effect is $1 / \hat{λ} = 45.6$ minutes.

Critical thinking. As in Example 8.3, these times are much more regular (SD 6.3 minutes) than an exponential model allows (SD equal to the mean, about 46 minutes). A toxic effect that requires a dose to build up is better described by a gamma or Weibull model, which can have a rising risk over time (Sections 5.4 and 5.6).

Example 8.7: Weekly Burglaries

A precinct records 3, 3, 5 and 8 burglaries in four weeks. Assuming the weekly counts are IID $Pois (λ)$ , find the MLE of $λ$ .

$L (λ) = \prod_{i = 1}^{4} \frac{e^{- λ} λ^{x_{i}}}{x_{i}!} = \frac{e^{- 4 λ} λ^{19}}{3! 3! 5! 8!},$ so $ℓ (λ)$ $= - 4 λ$ $+ 19 ln λ$ $- ln (3! 3! 5! 8!)$ and $ℓ' (λ)$ $= - 4$ $+ \frac{19}{λ}$ . Setting this to zero gives $\hat{λ} = \frac{19}{4} = 4.75$ burglaries a week, the sample mean. The constant $ln (3! 3! 5! 8!)$ doesn't depend on $λ$ , so it never affects the answer.

Find the maximum of the likelihood

The likelihood asks, for each possible parameter value, how probable the data we actually observed would be. Choose one of the chapter's datasets and slide the parameter. The top panel shows the data with the model for your chosen value; the bottom panel shows the log-likelihood. The best value is at the peak, the maximum likelihood estimate.

data parameter =

Example 8.8: The Tanks, by Maximum Likelihood ★ lecture example

Find the MLE of $N$ in the tank problem of Example 8.1, from $k$ distinct serial numbers with largest value $m$ .

Every set of $k$ distinct serial numbers from $1, \dots, N$ is equally likely, with probability $1 / (\binom{N}{k})$ , provided all the numbers are at most $N$ . So $L (N) = {\begin{matrix} 1 / (\binom{N}{k}) & if N \geq m, \\ 0 & if N < m . \end{matrix}$ For $N \geq m$ , $(\binom{N}{k})$ increases with $N$ , so $L (N)$ decreases. The likelihood is largest at the smallest possible value, $\hat{N} = m$ , the largest serial number seen. Calculus doesn't help here, since the maximum is at the edge of the allowed range.

So the MLE is the estimator ${\hat{N}}_{2}$ of Example 8.1, which is biased: it's never too high, and it's too low on average. Maximum likelihood is a superb general method, but its estimates are not always unbiased, especially in small samples. Here the fix ( $M (1$ $+ 1 / k)$ $- 1$ ) is simple.

The likelihood for the five captured tanks (largest serial number 61) is zero below 61 and falls steadily above it, so its peak is at the edge:

N <- 40:150
lik <- ifelse(N >= 61, 1 / choose(N, 5), 0)
ggplot(data.frame(N, lik = lik / max(lik)), aes(N, lik)) +
  geom_step(linewidth = 1) + geom_vline(xintercept = 61, linetype = "dashed") +
  labs(x = "number of tanks N", y = "likelihood (scaled)")

Properties of maximum likelihood

Why is maximum likelihood the standard method? Because for large samples, under mild conditions on the model, MLEs have three excellent properties.

Theorem: Large-Sample Properties of the MLE

For an IID sample from a model satisfying mild smoothness conditions, as $n \to \infty$ :

Consistent: $\hat{θ}$ converges to the true $θ$ .
Approximately Gaussian: $\hat{θ} \approx N (θ, {SE}^{2})$ , with the standard error estimated by the curvature of the log-likelihood at its peak: $SE [\hat{θ}] \approx \frac{1}{\sqrt{- ℓ ″ (\hat{θ})}} .$
Efficient: in the large-sample limit its variance reaches the Cramér–Rao lower bound (Section 8.7), so no well-behaved estimator does better.

Also, MLEs are invariant: the MLE of any function $g (θ)$ is $g (\hat{θ})$ .

The curvature formula makes intuitive sense: a sharply peaked log-likelihood (large $- ℓ ″$ ) means the data rule out all but a narrow range of $θ$ , so the estimate is precise. For the 2FA data, $- ℓ ″ (\hat{p})$ $= \frac{3}{{0.3}^{2}}$ $+ \frac{7}{{0.7}^{2}}$ $= \frac{10}{0.3 \times 0.7}$ $\approx 47.6,$ $SE$ $\approx \frac{1}{\sqrt{47.6}}$ $\approx 0.145,$ which is the familiar $\sqrt{\hat{p} (1 - \hat{p}) / n}$ of Section 7.5. For the toxic-effect times, $- ℓ ″ (\hat{λ})$ $= 5 / {\hat{λ}}^{2}$ , giving $SE$ $\approx \hat{λ} / \sqrt{5}$ $\approx 0.0098$ : with five observations, the rate is known only to within about 45%.

Example 8.9: Gaussian Maximum Likelihood

For an IID sample from $N (μ, σ^{2})$ with both parameters unknown, find the MLEs.

The log-likelihood is $ℓ (μ,$ $σ^{2})$ $= - \frac{n}{2} ln (2 π σ^{2})$ $- \frac{1}{2 σ^{2}} \sum_{i = 1}^{n} {(x_{i} - μ)}^{2} .$ Maximize over $μ$ first. Only the sum depends on $μ$ , and it's smallest at $μ = \overline{x}$ (Section 2.2: the mean minimizes the sum of squared deviations). So $\hat{μ} = \overline{x}$ , whatever $σ^{2}$ is.

Then over $σ^{2}$ . Writing $v = σ^{2}$ and $Q$ $= \sum {(x_{i} - \overline{x})}^{2}$ , the derivative of $- \frac{n}{2} ln v - \frac{Q}{2 v}$ with respect to $v$ is $- \frac{n}{2 v} + \frac{Q}{2 v^{2}}$ , which is zero at ${\hat{σ}}^{2} = \frac{1}{n} \sum_{i = 1}^{n} {(x_{i} - \overline{x})}^{2} .$ This divides by $n$ , so it is biased (Example 8.2), though the bias vanishes as $n$ grows. By invariance, the MLE of $σ$ is $\sqrt{{\hat{σ}}^{2}}$ . Least squares and maximum likelihood give the same $\hat{μ}$ for Gaussian data: that is Gauss's 1809 argument (Section 5.3), in modern dress.

Check Yourself

In 20 biometric database queries, 7 produce a match. Find the MLE of the match probability $p$ .

p̂ =

The log-likelihood for one parameter is very flat near its peak. What does that tell you? The MLE is wrong.The model doesn't fit.The data are consistent with a wide range of parameter values, so the standard error is large.

Think About It: The Likelihood of What, Exactly?

A defence lawyer says: "The likelihood of the evidence under my client's innocence is 0.01, so the probability that my client is innocent is only 1%." What's wrong with this argument? How is it related to the likelihood in this section?

It confuses $Pr (evidence ∣ innocent)$ with $Pr (innocent ∣ evidence)$ : the prosecutor's fallacy of Section 3.6, here used by the defence. A likelihood measures how well a hypothesis explains the data. To turn it into a probability of the hypothesis, Bayes' theorem needs the prior probability and the likelihood under the alternatives as well. In this section the likelihood compares parameter values by how well they explain the data, and the MLE picks the best explainer, without claiming that it has any particular probability of being true. Bayesian statistics (Section 8.6's box) is what adds the prior and turns likelihoods into probabilities about parameters.

Key Takeaways

The likelihood $L (θ)$ is the probability of the observed data as a function of $θ$ ; the MLE maximizes it (usually via the log-likelihood and calculus).
Common MLEs: $\hat{p} = x / n$ , $\hat{λ} = \overline{x}$ (Poisson), $\hat{λ} = 1 / \overline{x}$ (exponential), $\hat{μ} = \overline{x}$ and ${\hat{σ}}^{2}$ $= \frac{1}{n} \sum {(x_{i} - \overline{x})}^{2}$ (Gaussian).
For large samples MLEs are consistent, approximately Gaussian and efficient, with $SE$ $\approx 1 / \sqrt{- ℓ ″ (\hat{θ})}$ . In small samples they can be biased.

Lecture 21 — Confidence intervals for a mean

8.4 Confidence intervals for a mean

The idea

A point estimate is almost never exactly right. A confidence interval gives a range of plausible values, built so that the method captures the true value a known fraction of the time. Here's the recipe, using the Gaussian sampling distribution of Chapter 7.

If $\overline{X}$ $\approx N (μ,$ $σ^{2} / n)$ , then $Z = \frac{\overline{X} - μ}{σ / \sqrt{n}}$ is approximately standard Gaussian, so $Pr (- 1.96 \leq \frac{\overline{X} - μ}{σ / \sqrt{n}} \leq 1.96)$ $\approx 0.95 .$ Rearranging the inequalities to isolate $μ$ : $Pr (\overline{X} - 1.96 \frac{σ}{\sqrt{n}} \leq μ \leq \overline{X} + 1.96 \frac{σ}{\sqrt{n}})$ $\approx 0.95 .$ The interval $\overline{X}$ $\pm 1.96 σ / \sqrt{n}$ is random (it moves with $\overline{X}$ ), and it contains the fixed number $μ$ with probability 0.95.

Definition: Confidence Interval

A $100 (1 - α) %$ confidence interval for $θ$ is an interval computed from the sample by a rule that, over repeated samples, contains the true $θ$ with probability $1 - α$ , called the confidence level or coverage probability. Many intervals have the form $estimate$ $\pm (critical value)$ $\times (standard error),$ and the second term is called the margin of error.

The critical value $z_{α / 2}$ is the point with area $α / 2$ above it under the standard Gaussian: qnorm(1 - alpha/2). For 90%, 95% and 99% confidence it is 1.645, 1.960 and 2.576.

Case 1: Gaussian population, σ known

Theorem: Z Interval for a Mean

If $X_{1}, \dots, X_{n}$ is a random sample from $N (μ, σ^{2})$ with $σ$ known, then $\overline{x} \pm z_{α / 2} \frac{σ}{\sqrt{n}}$ is an exact $100 (1 - α) %$ confidence interval for $μ$ .

Example 8.10: A Well-Calibrated Instrument

An instrument measures the concentration of a heavy metal in water with Gaussian error of known SD $σ = 0.1$ mg/L, from long calibration studies. Forty measurements of one sample average $\overline{x} = 5.426$ mg/L. (a) Find 90%, 95% and 99% confidence intervals for the true concentration $μ$ . (b) How many measurements would give a 95% interval of total width 0.02 mg/L?

(a) The standard error is $0.1 / \sqrt{40}$ $\approx 0.0158$ .

confidence	zα/₂	margin of error	interval
90%	1.645	0.026	(5.400, 5.452)
95%	1.960	0.031	(5.395, 5.457)
99%	2.576	0.041	(5.385, 5.467)

se <- 0.1 / sqrt(40)
sapply(c(0.90, 0.95, 0.99), function(cl) 5.426 + c(-1, 1) * qnorm(1 - (1 - cl) / 2) * se)

output

         [,1]    [,2]     [,3]
[1,] 5.399993 5.39501 5.385273
[2,] 5.452007 5.45699 5.466727

More confidence requires a wider interval: to be more sure of catching $μ$ , you cast a wider net.

(b) Sample size. The width is $2$ $\times 1.96$ $\times 0.1 / \sqrt{n}$ . Setting it to 0.02 gives $\sqrt{n}$ $= 2$ $\times 1.96$ $\times 0.1 / 0.02$ $= 19.6$ , so $n = 384.2$ : 385 measurements. Halving the width of an interval needs four times the data, as always.

Case 2: Gaussian population, σ unknown: the t interval

Knowing $σ$ but not $μ$ is rare. Usually we must estimate $σ$ by the sample SD $s$ , and that adds uncertainty: in a small sample, $s$ may happen to be too small, making the interval too narrow. William Gosset ("Student", Section 5.6) found the exact correction in 1908; Fisher later gave the rigorous proof.

Theorem: Student's T

If $X_{1}, \dots, X_{n}$ is a random sample from $N (μ, σ^{2})$ , then $T = \frac{\overline{X} - μ}{S / \sqrt{n}}$ has a $t$ distribution with $n - 1$ degrees of freedom. Hence $\overline{x} \pm t_{n - 1, α / 2} \frac{s}{\sqrt{n}}$ is an exact $100 (1 - α) %$ confidence interval for $μ$ , where $t_{n - 1, α / 2}$ is qt(1 - alpha/2, n - 1).

The $t$ critical values are larger than the Gaussian ones, especially for small $n$ (with 9 degrees of freedom, $t_{9, 0.025} = 2.26$ instead of 1.96), and they approach the Gaussian values as $n$ grows (with 100 degrees of freedom, 1.98).

x <- seq(-4, 4, by = 0.02)
d <- rbind(data.frame(x, f = dnorm(x), dist = "Gaussian"),
           data.frame(x, f = dt(x, 9), dist = "t, 9 df"),
           data.frame(x, f = dt(x, 2), dist = "t, 2 df"))
d$dist <- factor(d$dist, levels = c("Gaussian", "t, 9 df", "t, 2 df"))
ggplot(d, aes(x, f, linetype = dist)) + geom_line(linewidth = 0.9) +
  scale_linetype_manual(values = c("solid", "dashed", "dotted")) +
  labs(x = NULL, y = "density", linetype = NULL)

Example 8.11: Student's Own Data ★ lecture example

In his 1908 paper, Student illustrated the $t$ distribution with data from a 1905 drug trial by Arthur Cushny and Alvin Peebles: 10 patients each took two sleeping drugs, and the extra hours of sleep were recorded. R's sleep data contain them. For each patient, take the difference between the two drugs. (a) Find a 95% confidence interval for the mean difference. (b) Find a 95% prediction interval for the difference in a new patient. Explain the difference between the two intervals.

(a) Paired differences. Each patient took both drugs, so we analyse one difference per patient, a single sample of 10:

d <- with(sleep, extra[group == 2] - extra[group == 1])
d
c(mean = mean(d), sd = sd(d), t_crit = qt(0.975, df = 9))
mean(d) + c(-1, 1) * qt(0.975, 9) * sd(d) / sqrt(10)
t.test(d)$conf.int                                  # the same, built in

output

 [1] 1.2 2.4 1.3 1.3 0.0 1.0 1.8 0.8 4.6 1.4
    mean       sd   t_crit 
1.580000 1.229995 2.262157 
[1] 0.7001142 2.4598858
[1] 0.7001142 2.4598858
attr(,"conf.level")
[1] 0.95

The differences have mean $\overline{d} = 1.58$ hours and SD $s = 1.23$ . With $t_{9, 0.025} = 2.262$ : $1.58$ $\pm 2.262$ $\times \frac{1.23}{\sqrt{10}}$ $= 1.58$ $\pm 0.88$ $= (0.70,$ $2.46) hours .$ The second drug gives between about 0.7 and 2.5 hours more sleep than the first, on average. The whole interval is above 0, so the data favour the second drug.

(b) A prediction interval for one new patient's difference $D_{new}$ must allow for two sources of uncertainty: how far $\overline{d}$ is from $μ$ , and how far one patient varies around $μ$ . Since $D_{new} - \overline{D}$ has variance $σ^{2} + σ^{2} / n$ , and assuming the individual differences are roughly Gaussian (a prediction for one patient can't lean on the CLT), $\overline{d}$ $\pm t_{9, 0.025} s \sqrt{1 + \frac{1}{n}}$ $= 1.58$ $\pm 2.262$ $\times 1.23$ $\times 1.049$ $= (- 1.34,$ $4.50) hours .$

mean(d) + c(-1, 1) * qt(0.975, 9) * sd(d) * sqrt(1 + 1/10)

output

[1] -1.338251  4.498251

The two goals. The confidence interval is about inference: where is the average effect? It shrinks to nothing as $n$ grows. The prediction interval is about prediction: what will happen to the next patient? It can never be narrower than the spread of individual patients, however much data we collect. A drug can help on average and still do nothing, or harm, for a particular person.

ci <- mean(d) + c(-1, 1) * qt(0.975, 9) * sd(d) / sqrt(10)
pi <- mean(d) + c(-1, 1) * qt(0.975, 9) * sd(d) * sqrt(1 + 1/10)
ggplot() +
  geom_point(aes(x = d, y = 3), size = 3, alpha = 0.7) +
  geom_segment(aes(x = ci[1], xend = ci[2], y = 2, yend = 2), linewidth = 2, colour = "steelblue") +
  geom_segment(aes(x = pi[1], xend = pi[2], y = 1, yend = 1), linewidth = 2, colour = "orange") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_y_continuous(breaks = 1:3, labels = c("prediction interval", "confidence interval", "data"),
                     limits = c(0.5, 3.5)) +
  labs(x = "extra sleep (hours): drug 2 minus drug 1", y = NULL)

Case 3: any population, large n

If the population isn't Gaussian but $n$ is large, the Central Limit Theorem makes $\overline{X}$ approximately Gaussian, and $s$ is close to $σ$ . Then $\overline{x} \pm z_{α / 2} \frac{s}{\sqrt{n}}$ is an approximate $100 (1 - α) %$ confidence interval for $μ$ . In practice most software uses the $t$ critical value anyway (it makes almost no difference for large $n$ , and it's exact when the population is Gaussian), which is what R's t.test does. How large "large" must be depends on the population's skewness (Section 7.4).

Example 8.12: Michelson and the Speed of Light ★ lecture example

In 1879, Albert Michelson made 100 measurements of the speed of light (R's morley, Section 2.3), recorded in km/s minus 299,000. (a) Find a 95% confidence interval for the mean of his measurements. (b) Michelson measured in air, where light travels at about 299,734.5 km/s (Section 2.3; in a vacuum it is exactly 299,792.458 km/s). Does the interval contain the true value? What went wrong?

(a) With $n = 100$ , $\overline{x} = 852.4$ and $s = 79.0$ : $852.4$ $\pm 1.984$ $\times \frac{79.0}{\sqrt{100}}$ $= 852.4$ $\pm 15.7$ $= (836.7,$ $868.1),$ that is, between 299,836.7 and 299,868.1 km/s.

t.test(morley$Speed)$conf.int

output

[1] 836.7226 868.0774
attr(,"conf.level")
[1] 0.95

(b) No. The true value, 734.5 on this scale, is about 102 km/s below the bottom of the interval, and $(852.4$ $- 734.5) / 7.9$ $\approx 15$ standard errors from the mean. That's far too far to be bad luck.

What went wrong is systematic error: something in the apparatus, its calibration or the corrections applied to the readings pushed every measurement upwards; historians of science have debated the exact cause. A confidence interval accounts only for random error, the scatter of measurements around their own average. It assumes that the average of the measurements is centred on the truth. If the instrument is biased, the interval is centred in the wrong place, and more measurements just give a narrower interval around the wrong value. Michelson's later measurements, with improved methods, came much closer. Every interval in this chapter carries the same unspoken assumption: no bias.

Check Yourself

Find the critical value for a 95% $t$ interval with a sample of 10. Use R: qt(0.975, 9).

t =

A 95% confidence interval for a mean is (12.1, 15.3). Which change would make it narrower, other things being equal? A larger sample.A higher confidence level, such as 99%.A population with a larger SD.

Think About It: A Precise Wrong Answer

A breath-alcohol device is calibrated once a year. Between calibrations it drifts, reading on average up to 0.005 g/dL too high. A driver's three readings give a 95% confidence interval for their true level of (0.081, 0.085). Should the driver be convicted of being over 0.08? What would you want to know first?

The interval is built on the assumption that readings scatter around the true level. If the device reads up to 0.005 too high, the true level could be anywhere from about 0.076 to 0.085, which includes values under the limit. The interval's narrowness reflects the consistency of the three readings, not their accuracy. You'd want the device's calibration record, the date of its last calibration, and ideally a blood test analysed by a different method. Many jurisdictions require exactly these safeguards. Like Michelson's interval, a confidence interval is only as good as the assumption of no bias.

Key Takeaways

A confidence interval is estimate ± critical value × standard error, built so that the method covers the truth $100 (1 - α) %$ of the time.
Gaussian population, σ known: $\overline{x} \pm z_{α / 2} σ / \sqrt{n}$ . σ unknown: $\overline{x}$ $\pm t_{n - 1, α / 2} s / \sqrt{n}$ . Large $n$ , any population: approximately the same.
A prediction interval for one new observation, $\overline{x}$ $\pm t s \sqrt{1 + 1 / n}$ , is much wider and never shrinks to zero. It needs the individual observations to be roughly Gaussian; the CLT doesn't rescue it.
Intervals capture random error only; they don't detect bias.

Lecture 22 — Intervals for proportions and what confidence means

8.5 Confidence intervals for a proportion

For a proportion, Section 7.5 gave $\hat{p}$ $\approx N (p,$ $p (1$ $- p) / n)$ . The difficulty is that the standard error depends on the unknown $p$ . There are two common ways to handle it.

Theorem: Wald and Wilson Intervals for a Proportion

With $x$ successes in $n$ trials and $\hat{p} = x / n$ :

The Wald (or traditional) interval replaces $p$ by $\hat{p}$ in the standard error: $\hat{p} \pm z_{α / 2} \sqrt{\frac{\hat{p} (1 - \hat{p})}{n}} .$
The Wilson (or score) interval keeps $p$ in the standard error and solves $| \hat{p} - p |$ $\leq z_{α / 2} \sqrt{p (1 - p) / n}$ for $p$ . Writing $z = z_{α / 2}$ , the solution is $\frac{\hat{p} + \frac{z^{2}}{2 n}}{1 + \frac{z^{2}}{n}} \pm \frac{z}{1 + \frac{z^{2}}{n}} \sqrt{\frac{\hat{p} (1 - \hat{p})}{n} + \frac{z^{2}}{4 n^{2}}} .$

For Wilson's interval, square both sides of $| \hat{p} - p |$ $\leq z \sqrt{p (1 - p) / n}$ to get the quadratic inequality ${(\hat{p} - p)}^{2}$ $\leq \frac{z^{2}}{n} p (1$ $- p)$ , that is, $(1 + \frac{z^{2}}{n}) p^{2}$ $- (2 \hat{p} + \frac{z^{2}}{n}) p$ $+ {\hat{p}}^{2}$ $\leq 0 .$ The set of $p$ satisfying it lies between the two roots of the quadratic, and the quadratic formula gives the expression above. $◼$ ∎

The Wilson interval's centre is pulled from $\hat{p}$ slightly towards $\frac{1}{2}$ , and it never extends below 0 or above 1. For large $n$ and $\hat{p}$ not near 0 or 1, the two intervals nearly coincide. For small $n$ or extreme $\hat{p}$ , Wilson's is much better: the Wald interval's actual coverage can fall well below its nominal 95%. In R, prop.test(x, n, correct = FALSE)$conf.int gives the Wilson interval.

Example 8.13: Pesticide Residues ★ lecture example

A lab tests a random sample of 150 agricultural products from a shipment, and 42 exceed the legal limit for pesticide residues. (a) Find the 95% Wald and Wilson intervals for the proportion $p$ over the limit in the whole shipment. (b) Compare them.

(a) $\hat{p}$ $= 42 / 150$ $= 0.28$ , and the Wald standard error is $\sqrt{0.28 \times 0.72 / 150}$ $\approx 0.0367$ : $Wald: 0.28$ $\pm 1.96$ $\times 0.0367$ $= (0.208,$ $0.352) .$ For Wilson, with $z^{2} / n$ $= 3.8416 / 150$ $\approx 0.0256$ : the centre is $\frac{0.28 + 0.0128}{1.0256}$ $\approx 0.2855$ , and the half-width is $\frac{1.96}{1.0256} \sqrt{0.001344 + 0.000043}$ $\approx 0.0712$ , giving $Wilson: (0.214, 0.357) .$

x <- 42; n <- 150; p_hat <- x / n
p_hat + c(-1, 1) * qnorm(0.975) * sqrt(p_hat * (1 - p_hat) / n)    # Wald
prop.test(x, n, correct = FALSE)$conf.int                            # Wilson

output

[1] 0.2081465 0.3518535
[1] 0.2143304 0.3566565
attr(,"conf.level")
[1] 0.95

(b) With $n = 150$ and $\hat{p} = 0.28$ , the intervals differ by only about 0.005 to 0.006 at each end; Wilson's is shifted slightly towards 0.5. Either way, between about a fifth and a third of the shipment is over the limit, which is far too many to release it.

Example 8.14: Zero Out of Twenty ★ lecture example

A new drug is given to 20 patients, and none has a serious side effect. (a) What does the Wald interval give? (b) What does the Wilson interval give? (c) The "rule of three" says that, after 0 events in $n$ trials, an approximate 95% upper bound for $p$ is $3 / n$ . Where does it come from?

(a) $\hat{p} = 0$ , so the Wald standard error is 0, and the interval is $[0, 0]$ , which claims certainty that the side effect never happens. That's absurd: 20 patients can't rule out a 1-in-50 risk.

(b) prop.test(0, 20, correct = FALSE)$conf.int gives $[0, 0.161]$ : the risk could be as high as about 16%.

(c) If the true risk is $p$ , the chance of 0 events in $n$ patients is ${(1 - p)}^{n}$ . The largest $p$ that still gives this outcome at least a 5% chance solves ${(1 - p)}^{n}$ $= 0.05$ , so $p$ $= 1$ $- {0.05}^{1 / n}$ . For $n = 20$ that's 0.139. Since $ln 0.05 \approx - 3$ and ${(1 - p)}^{n}$ $\approx e^{- n p}$ for small $p$ , the bound is about $3 / n = 0.15$ .

prop.test(0, 20, correct = FALSE)$conf.int
1 - 0.05^(1/20)          # exact one-sided 95% upper bound
3 / 20                   # rule of three

output

[1] 0.0000000 0.1611252
attr(,"conf.level")
[1] 0.95
[1] 0.1391083
[1] 0.15

This is why drug regulators require large trials and post-marketing surveillance: to detect a side effect that occurs in 1 patient in 1,000, you need to watch thousands of patients.

How badly can the Wald interval fail? A simulation for $n = 20$ and $p = 0.1$ shows the Wald "95%" interval covering the truth only about 87% of the time, while Wilson's covers it about 96% of the time. For $n = 1000$ and $p = 0.5$ , both cover about 95%. Try it in the widget below.

The coverage of each method can be computed exactly, for every possible true $p$ , by adding the binomial probabilities of all the outcomes whose interval contains $p$ :

cover <- function(p, n = 20, z = 1.96) {
  x <- 0:n; ph <- x / n
  wald <- abs(ph - p) <= z * sqrt(ph * (1 - ph) / n)
  ctr <- (ph + z^2 / (2 * n)) / (1 + z^2 / n)
  hw <- z / (1 + z^2 / n) * sqrt(ph * (1 - ph) / n + z^2 / (4 * n^2))
  wil <- abs(ctr - p) <= hw
  c(Wald = sum(dbinom(x, n, p)[wald]), Wilson = sum(dbinom(x, n, p)[wil])) }
p <- seq(0.01, 0.99, by = 0.002)
cv <- t(sapply(p, cover))
d <- rbind(data.frame(p, coverage = cv[, "Wald"], method = "Wald"),
           data.frame(p, coverage = cv[, "Wilson"], method = "Wilson"))
ggplot(d, aes(p, coverage, colour = method)) + geom_line() +
  geom_hline(yintercept = 0.95, linetype = "dashed") + coord_cartesian(ylim = c(0.5, 1)) +
  scale_colour_manual(values = c(Wald = "steelblue", Wilson = "orange")) +
  labs(x = "true proportion p (n = 20)", y = "actual coverage", colour = NULL)

Check Yourself

How many people must a poll survey for a 95% Wald margin of error of at most 3 percentage points, in the worst case $p = 0.5$ ? (Round up.)

n =

Think About It: Whose Rate Is It?

A hospital reports that 0 of its last 25 hip replacements became infected, and advertises "a 0% infection rate". A competitor reports 3 infections in 400 operations (0.75%). Which hospital has better evidence of a low infection rate? What would you want to know besides the numbers?

The rule of three gives the first hospital an upper bound of about $3 / 25 = 12 %$ : its data are consistent with a rate far higher than the competitor's. The competitor's Wilson interval is roughly 0.3% to 2.2%, much tighter. A "0% rate" from a small number of operations is weak evidence. Beyond the numbers, you'd want to know whether the hospitals treat similar patients (case mix), how infections are counted and over what follow-up period, and whether infections diagnosed elsewhere are included. Rates from different denominators and definitions aren't comparable, the same lesson as Section 7.5's funnel plots.

Key Takeaways

Wald: $\hat{p}$ $\pm z \sqrt{\hat{p} (1 - \hat{p}) / n}$ ; simple, but poor for small $n$ or $\hat{p}$ near 0 or 1.
Wilson: solves for $p$ inside the standard error; better coverage, stays within $[0, 1]$ . In R, prop.test(x, n, correct = FALSE).
After 0 events in $n$ trials, the 95% upper bound is about $3 / n$ .

8.6 What "95% confidence" means

After computing a 95% interval such as (0.70, 2.46) hours, it's tempting to say "there's a 95% probability that the true mean is between 0.70 and 2.46". In classical statistics, that statement isn't correct. The true mean $μ$ is a fixed number; it's either in that interval or it isn't. The 95% describes the method: if we repeated the study many times, computing an interval each time, 95% of those intervals would contain $μ$ .

What "95% confidence" means

Here we know the truth, which never happens in real life. Each horizontal line is a confidence interval computed from a fresh random sample. Intervals that capture the true value are blue; those that miss are red. Any single interval either covers the truth or it doesn't. The confidence level is the long-run proportion that do. Try the proportion intervals with a true p near 0 or 1, and compare the traditional and score methods.

interval confidence level = % sample size n = true p (proportions) =

Definition: Interpreting a Confidence Interval

Correct: "This interval was produced by a method that captures the true value in 95% of samples." Or: "Before the data were collected, there was a 95% probability that the interval we were about to compute would contain the true value."

Incorrect: "There is a 95% probability that $μ$ lies in (0.70, 2.46)"; "95% of the data lie in the interval"; "a repeat study would have a 95% chance of a mean in this interval".

The distinction matters most when we have information beyond the data. A 95% interval for the speed of light that misses the known value (Example 8.12) is not "95% likely" to contain it; we know it doesn't. The coverage guarantee is about the long run of intervals, and Michelson's was one of the 5%, or (more likely) its assumptions failed.

Bayesian Inference: Credible Intervals: The Probability Statement People Want

A Bayesian treats the unknown parameter as uncertain, gives it a prior distribution, and updates it with the data (Section 5.5). The middle 95% of the posterior distribution is a 95% credible interval, and for it the natural statement is correct: given the model, the prior and the data, there's a 95% probability that the parameter lies in the interval.

For the pesticide data (42 of 150 over the limit), a uniform $Beta (1, 1)$ prior gives the posterior $Beta (43, 109)$ :

qbeta(c(0.025, 0.975), 43, 109)          # 95% credible interval: 0.214 to 0.357
prop.test(42, 150, correct = FALSE)$conf.int   # 95% Wilson interval: 0.214 to 0.357

output

[1] 0.2143697 0.3568306
[1] 0.2143304 0.3566565
attr(,"conf.level")
[1] 0.95

The two intervals agree to three decimal places, though they mean different things. With plenty of data and a flat prior, Bayesian and classical intervals usually agree numerically (Section 7.5's box). They part company when data are scarce, when there is genuine prior information (the speed of light, a well-studied drug), or when the prior is strong. Which interpretation you prefer is one of the oldest debates in statistics.

p <- seq(0.1, 0.5, by = 0.001); d <- data.frame(p, f = dbeta(p, 43, 109))
ci <- qbeta(c(0.025, 0.975), 43, 109)
ggplot(d, aes(p, f)) +
  geom_area(data = subset(d, p >= ci[1] & p <= ci[2]), fill = "steelblue", alpha = 0.35) +
  geom_line(linewidth = 1) +
  annotate("text", x = 0.285, y = 3, label = "95% of the\nposterior") +
  labs(x = "proportion p over the limit", y = "posterior density")

History: Neyman, Fisher and the Meaning of Confidence

The theory of confidence intervals was set out in 1937 by the Polish-American statistician Jerzy Neyman, who insisted that probability statements should describe procedures, not fixed parameters. Ronald Fisher (Section 4.6) had proposed a rival idea, "fiducial probability", which tried to make probability statements about parameters without a prior; the two men disagreed bitterly for decades, and fiducial inference never became standard. The Wilson interval is older still: the American mathematician Edwin Bidwell Wilson published it in 1927, and it was largely ignored for the rest of the century, while textbooks taught the Wald interval, until a series of papers around 2000 showed how badly the Wald interval can perform.

Check Yourself

A study reports a 95% confidence interval of (2.1, 3.4) for a mean. Which interpretation is correct? There is a 95% probability that the true mean is between 2.1 and 3.4.95% of individuals have values between 2.1 and 3.4.The interval came from a method that captures the true mean in 95% of repeated samples.

Think About It: Reading a News Report

A news story says: "A new study found that people who drink coffee live 1.8 years longer, with a 95% confidence interval of 0.2 to 3.4 years." List three questions you'd ask before believing that coffee extends life.

(1) Was this an experiment or an observational study? If people chose whether to drink coffee, coffee drinkers may differ in income, health or habits (confounding, Chapter 1), and the interval measures only random error, not this bias. (2) How many outcomes did the researchers examine? If they looked at dozens of foods and outcomes, one interval excluding zero is expected by chance (Section 4.4's Paul the octopus, again). (3) Does the interval describe an average or an individual? Even if the average effect is real, an individual coffee drinker's lifespan varies by decades. A confidence interval answers one narrow question, about random sampling error, under the study's assumptions.

Key Takeaways

"95% confidence" describes the long-run success rate of the method, not the probability that a particular interval contains $θ$ .
A Bayesian credible interval does carry that probability, given a prior; with lots of data the two often agree numerically.
Coverage is guaranteed only if the model's assumptions (random sampling, no bias, the right distribution) hold.

Prediction and inference in this chapter

Inference. This chapter is inference at its core. An estimator turns data into a guess at an unseen parameter, and its bias, standard error and MSE say how good the rule is. The method of moments and maximum likelihood are general recipes for building estimators, and confidence intervals attach an honest margin of error. We estimated tank production, a speed, a drug effect, a rate of residues and a side-effect risk.
Prediction. Every prediction from a model uses estimated parameters. The prediction interval of Example 8.11 shows that predicting one new case carries two uncertainties, the estimate's and the individual's, and that the second never goes away.
The limits. Intervals cover random error, not bias (Michelson, the drifting breathalyser), and their coverage depends on the model's assumptions. "95% confidence" is a statement about a procedure, which is easy to misreport.
What comes next. Chapter 9 asks a related question: instead of "what values of the parameter are plausible?", it asks "is one particular value plausible?". That is a hypothesis test, and there is a close link: a 95% confidence interval contains exactly the values that the corresponding two-sided test at the 5% level would not reject.

8.7 Interludes

Interlude (Optional): Sufficient Statistics: Data Compression Without Loss

In Example 8.5, the MLE of $p$ depended only on the number of 2FA users, not on which people they were. In Example 8.7, only the total number of burglaries mattered. A statistic $T$ is sufficient for $θ$ if, once you know $T$ , the rest of the data carry no more information about $θ$ . Fisher introduced the idea in 1922, and the factorization theorem gives a test: $T$ is sufficient if the likelihood can be written as $g (T,$ $θ)$ $\times h (data)$ , a part that depends on $θ$ only through $T$ , times a part that doesn't involve $θ$ at all. For the Poisson, $L (λ)$ $= e^{- n λ} λ^{\sum x_{i}}$ $\times \frac{1}{\prod x_{i}!}$ , so $\sum x_{i}$ is sufficient. The "exponential families" of Section 5.8 are the main family of models with simple sufficient statistics of fixed size, however large the sample. (When the range of possible values depends on the parameter, other models can have them too: in the tank problem, the largest serial number $M$ is sufficient for $N$ .)

Why it matters: sufficiency tells you what to record. A lab that stores only the total count and the number of trials has lost nothing for estimating a Poisson or binomial rate; but it has lost the ability to check whether the model fits, which needs the individual values.

Interlude (Optional): The Cramér–Rao Bound: A Speed Limit for Estimators

How precise can an unbiased estimator possibly be? In the 1940s Harald Cramér and C. R. Rao independently proved a lower bound: for an unbiased estimator from an IID sample of size $n$ , $Var [\hat{θ}]$ $\geq \frac{1}{n I (θ)},$ $I (θ)$ $= - E [\frac{d^{2}}{d θ^{2}} ln f (X; θ)],$ where $I (θ)$ is the Fisher information in one observation. For a Bernoulli trial, $I (p) = \frac{1}{p (1 - p)}$ , so no unbiased estimator of $p$ can have variance below $\frac{p (1 - p)}{n}$ , and $\hat{p}$ achieves it exactly. The large-sample SE formula of Section 8.3, $1 / \sqrt{- ℓ ″ (\hat{θ})}$ , is an estimate of the square root of this bound from the data (using the observed curvature instead of its expected value), which is the precise sense in which MLEs are "efficient".

Why it matters: the bound tells you when to stop looking for a better estimator, and how much data an experiment needs before it starts.

Interlude (Optional): Maximum Likelihood for Classification: Logistic Regression

Chapter 1 asked whether a mother's smoking is associated with low birth weight (MASS::birthwt). Logistic regression models the probability of a yes/no outcome through its log-odds: $ln \frac{p}{1 - p}$ $= β_{0}$ $+ β_{1} x$ , where $x$ is a predictor. The coefficients are fitted by maximum likelihood, since each baby contributes a Bernoulli factor $p_{i}^{y_{i}} {(1 - p_{i})}^{1 - y_{i}}$ to the likelihood.

fit <- glm(low ~ smoke, data = MASS::birthwt, family = binomial)
coef(fit)
exp(coef(fit)["smoke"])           # odds ratio for smokers

output

(Intercept)       smoke 
 -1.0870515   0.7040592 
   smoke 
2.021944

The fitted odds ratio is about 2.0: the odds of a low-birth-weight baby are about twice as high for smokers. (With one yes/no predictor, this is exactly the odds ratio in the 2 × 2 table: $(30 / 44) / (29 / 86)$ $\approx 2.02$ .) Add more predictors and the same MLE machinery adjusts for confounders, and with enough predictors it becomes the classifier at the core of many spam filters, credit models and medical risk scores.

Why it matters: logistic regression is maximum likelihood in action, and it is often the first method a data scientist reaches for. Its coefficients are estimates with standard errors, so everything in this chapter applies, including the warning that they measure association, not cause.

8.8 Exercises

Each exercise has a hint and an answer. For interval questions, state which case applies (σ known, σ unknown, large sample, proportion) before calculating.

Exercises

Time to a toxic reaction. The times in minutes for a chemical to produce a toxic reaction in 7 mice are 4.2, 5.1, 6.7, 4.9, 5.5, 6.2, 5.0. Assuming an exponential model, find the MLE of $λ$ . Is the exponential model plausible?

$\hat{λ} = n / \sum x_{i}$ . Compare the sample SD with the sample mean.

$\sum x_{i} = 37.6$ , so $\hat{λ}$ $= 7 / 37.6$ $\approx 0.186$ per minute. Not very plausible: the SD (0.84) is far smaller than the mean (5.37), whereas an exponential's SD equals its mean.
Intrusion flags. Of 6 intrusion attempts, a system flagged the 1st, 4th, 5th and 6th. Find the MLE of the detection probability, and its approximate standard error.

$\hat{p} = x / n$ and $SE$ $\approx \sqrt{\hat{p} (1 - \hat{p}) / n}$ .

$\hat{p}$ $= 4 / 6$ $\approx 0.667$ , with SE $\approx \sqrt{0.667 \times 0.333 / 6}$ $\approx 0.19$ : with 6 attempts, the estimate is very imprecise.
A geometric MLE. A phishing tester sends emails until someone clicks. In four campaigns the first click came on email 3, 1, 4 and 2. Assuming each email is clicked independently with probability $p$ , find the MLE of $p$ .

The probability that the first click is on email $x$ is ${(1 - p)}^{x - 1} p$ . Write the likelihood for all four campaigns.

$L (p)$ $= p^{4} {(1 - p)}^{\sum (x_{i} - 1)}$ $= p^{4} {(1 - p)}^{6}$ , so $ℓ' (p)$ $= \frac{4}{p}$ $- \frac{6}{1 - p}$ $= 0$ gives $\hat{p} = 0.4$ , which is $n / \sum x_{i}$ $= 4 / 10$ .
Counting tanks. Four captured tanks have serial numbers 37, 82, 15 and 104. Estimate $N$ with each of the three rules of Example 8.1. Which would you report?

$\overline{x} = 59.5$ , $m = 104$ , $k = 4$ .

Twice the mean: 118. Largest: 104. Largest plus gap: $104$ $\times 1.25$ $- 1$ $= 129$ . Report 129: it's unbiased and has the smallest typical error.
Chemical concentration (σ known). Ten samples have concentrations 13.2, 12.8, 14.0, 13.5, 13.1, 13.3, 12.9, 13.6, 13.7, 13.0 ppm. The population is Gaussian with known SD 0.4 ppm. Find a 95% confidence interval for the mean.

$\overline{x} = 13.31$ ; SE $= 0.4 / \sqrt{10}$ .

$13.31$ $\pm 1.96$ $\times 0.126$ $= (13.06,$ $13.56)$ ppm.
Time to a verdict (σ unknown). Eight court cases took 41, 39, 43, 38, 37, 42, 40 and 44 days to reach a verdict. Assuming a Gaussian population, find a 90% confidence interval for the mean.

$\overline{x} = 40.5$ , $s = 2.45$ , $t_{7,,0.05} = $ qt(0.95, 7) $= 1.895$ .

$40.5$ $\pm 1.895$ $\times 2.45 / \sqrt{8}$ $= 40.5$ $\pm 1.64$ $= (38.86,$ $42.14)$ days.
Load times (large sample). Fifty measurements of a website's load time give $\overline{x} = 2.18$ seconds and $s = 0.24$ . Find a 95% confidence interval for the mean load time. Why is a Gaussian population not needed?

With $n = 50$ , the CLT makes $\overline{X}$ approximately Gaussian.

$2.18$ $\pm 1.96$ $\times 0.24 / \sqrt{50}$ $= (2.113,$ $2.247)$ seconds. The CLT covers the sampling distribution of the mean, whatever the shape of the population (unless it is extremely skewed).
Phishing detection. A filter correctly flags 72 of 100 phishing emails. Find 95% Wald and Wilson intervals for its detection rate, and compare.

Use the formulas of Section 8.5, or prop.test(72, 100, correct = FALSE).

Wald $(0.632, 0.808)$ ; Wilson $(0.625, 0.799)$ . Similar, with Wilson's shifted slightly towards 0.5.
Blood alcohol. Six forensic cases had blood alcohol levels 0.082, 0.089, 0.075, 0.081, 0.085, 0.080 g/dL. Assume a Gaussian population. (a) Estimate $μ$ and $σ$ by the method of moments. (b) Are these estimates unbiased? (c) Find a 95% confidence interval for $μ$ .

The method of moments divides by $n$ for the variance. For (c), use $s$ (dividing by $n - 1$ ) and $t_{5, 0.025} = 2.571$ .

(a) $\hat{μ} = 0.0820$ and $\hat{σ} = 0.00432$ . (b) $\hat{μ}$ is unbiased; ${\hat{σ}}^{2}$ is biased low (Example 8.2), and so is $\hat{σ}$ . (Even $s$ , though $s^{2}$ is unbiased for $σ^{2}$ , is slightly biased low for $σ$ .) (c) With $s = 0.00473$ : $0.082$ $\pm 2.571$ $\times 0.00473 / \sqrt{6}$ $= (0.0770,$ $0.0870)$ .
Two estimators of a Poisson rate. For a Poisson population, both $\overline{X}$ and $S^{2}$ are unbiased estimators of $λ$ (why?). Simulate 20,000 samples of size 20 from $Pois (4)$ and compare their standard errors. Which would you use?

A Poisson's mean and variance both equal $λ$ . Use replicate with c(mean(x), var(x)).

Both average about 4, but the SE of $\overline{X}$ is about 0.45 and of $S^{2}$ about 1.37. Use $\overline{X}$ , which is also the MLE.
Coverage by simulation. Simulate 10,000 samples of size 20 from a population with $p = 0.1$ , compute the 95% Wald and Wilson intervals for each, and find the fraction of each that contains 0.1.

x <- rbinom(10000, 20, 0.1), then compute both intervals for each x with vector arithmetic.

About 0.87 for Wald and 0.96 for Wilson. The Wald interval misses far more often than its 5% promise.
Galton's sons. For the 481 sons, $\overline{x} = 69.23$ and $s = 2.62$ inches. (a) Find a 95% confidence interval for the mean height of the population they represent. (b) Find a 95% prediction interval for one more son's height. (c) Why might the interval in (a) be too narrow? (Hint: Section 7.1's warning box.)

For (b), use $s \sqrt{1 + 1 / n}$ .

(a) $69.23$ $\pm 1.965$ $\times 0.12$ $= (69.0,$ $69.5)$ . (b) $69.23$ $\pm 1.965$ $\times 2.62$ $\times 1.001$ $\approx (64.1,$ $74.4)$ . (c) Brothers share families, so the 481 heights aren't independent; the effective sample size is smaller than 481, and the true standard error is larger.
Planning a study. A toxicologist wants a 95% interval of total width 1 µg/dL for the mean blood lead level of workers, whose SD is believed to be about 5. How many workers should be sampled? What if the SD is actually 7?

Width $= 2$ $\times 1.96$ $\times σ / \sqrt{n}$ .

$n$ $= {(2 \times 1.96 \times 5 / 1)}^{2}$ $= 384.2$ , so 385 workers. If $σ = 7$ , the interval from 385 workers would be $1 \times 7 / 5 = 1.4$ µg/dL wide; achieving width 1 would need ${(2 \times 1.96 \times 7)}^{2}$ $\approx 753$ workers.
Interpreting intervals. For each statement about a 95% confidence interval of (12, 18) for a city's mean commute time in minutes, say whether it's correct, and fix it if not. (a) "95% of commuters take between 12 and 18 minutes." (b) "If we repeated the survey many times, about 95% of the intervals would contain the true mean." (c) "There's a 5% chance the true mean is above 18."

Is each statement about individuals, about the method or about the parameter?

(a) Wrong: that's a statement about individuals, which needs a prediction interval (much wider). (b) Correct. (c) Wrong in classical statistics: the true mean is fixed. A Bayesian credible interval could support a statement like it, given a prior, but even then a central 95% interval leaves 2.5% above 18, not 5%.

Chapter 8Estimation

Think About It: Before You Start

Lecture 19 — Estimators and the method of moments

8.1 Estimators and their errors

Parameters, estimators and estimates

Definition: Estimator and Estimate

Definition: Bias, Standard Error and Mean Squared Error

Theorem: MSE = Variance + Bias²

Example 8.1: How Many Tanks? ★ lecture example

History: Counting Tanks from Serial Numbers

Why the sample variance divides by n − 1

Theorem: The Sample Variance Is Unbiased

Example 8.2: Is Unbiased Always Better?

Check Yourself

Think About It: A Biased Scale

Key Takeaways

8.2 The method of moments

Definition: Method of Moments

Example 8.3: Suspicious Logins ★ lecture example

Example 8.4: A Model for Rainfall ★ lecture example

Check Yourself

Key Takeaways

Lecture 20 — Maximum likelihood

8.3 Maximum likelihood

The likelihood

Definition: Likelihood and Maximum Likelihood Estimator

Example 8.5: Who Uses Two-Factor Authentication? ★ lecture example

Example 8.6: Time to a Toxic Effect ★ lecture example

Example 8.7: Weekly Burglaries

Example 8.8: The Tanks, by Maximum Likelihood ★ lecture example

Properties of maximum likelihood

Theorem: Large-Sample Properties of the MLE

Example 8.9: Gaussian Maximum Likelihood

Check Yourself

Think About It: The Likelihood of What, Exactly?

Key Takeaways

Lecture 21 — Confidence intervals for a mean

8.4 Confidence intervals for a mean

The idea

Definition: Confidence Interval

Case 1: Gaussian population, σ known

Theorem: Z Interval for a Mean

Example 8.10: A Well-Calibrated Instrument

Case 2: Gaussian population, σ unknown: the t interval

Theorem: Student's T

Example 8.11: Student's Own Data ★ lecture example

Case 3: any population, large n

Example 8.12: Michelson and the Speed of Light ★ lecture example

Check Yourself

Think About It: A Precise Wrong Answer

Key Takeaways

Lecture 22 — Intervals for proportions and what confidence means

8.5 Confidence intervals for a proportion

Theorem: Wald and Wilson Intervals for a Proportion

Example 8.13: Pesticide Residues ★ lecture example

Example 8.14: Zero Out of Twenty ★ lecture example

Check Yourself

Think About It: Whose Rate Is It?

Key Takeaways

8.6 What "95% confidence" means

Definition: Interpreting a Confidence Interval

Bayesian Inference: Credible Intervals: The Probability Statement People Want

History: Neyman, Fisher and the Meaning of Confidence

Check Yourself

Think About It: Reading a News Report

Key Takeaways

Prediction and inference in this chapter

8.7 Interludes

Interlude (Optional): Sufficient Statistics: Data Compression Without Loss

Interlude (Optional): The Cramér–Rao Bound: A Speed Limit for Estimators

Interlude (Optional): Maximum Likelihood for Classification: Logistic Regression

8.8 Exercises

Exercises

Chapter 8
Estimation