Probability and Statistics · Chapter 7

Chapter 7
Sampling distributions

Lectures 17–18 · why averages are predictable: statistics as random variables, the Law of Large Numbers, the Central Limit Theorem and the bootstrap

Cartoon: the best little schoolA stick figure holds up a trophy and says: Our tiny school has the best scores in the state! A second figure holding a newspaper replies: It also had the worst, two years ago. Our tiny school hasthe best scores in the state!It also had theworst, two years ago.
small samples produce the most extreme averages, in both directions

In the 1990s, researchers looking at US maps of kidney cancer noticed something striking. The counties with the lowest death rates were mostly small, rural counties in the Midwest, the South and the West. Clean air, fresh food, an outdoor life? Then they looked at the counties with the highest death rates. They were also mostly small, rural counties in the Midwest, the South and the West. Nothing about rural life can both protect people and endanger them. The explanation is a fact about averages: a rate computed from a few hundred people swings much further from the truth than a rate computed from a few million. This chapter is about exactly how far.

Every number we compute from data (a mean, a proportion, a median, a correlation) is a statistic, and a statistic computed from a random sample is a random variable: another sample would give a different value. Its distribution, the sampling distribution, is the bridge between probability and statistics. It serves both of the course's goals.

Here is the route through the chapter.

Think About It: Before You Start

A town has two hospitals. In the larger one about 45 babies are born each day, and in the smaller one about 15. About half of all babies are boys, but on any particular day the fraction varies. For a year, each hospital records the days on which more than 60% of the babies born were boys. Which hospital records more such days: the larger, the smaller, or about the same? Decide before reading on.

The smaller hospital, by a wide margin: a day with more than 60% boys happens about 15% of the time with 15 births, but only about 7% of the time with 45. When the psychologists Amos Tversky and Daniel Kahneman put this question to university students in 1972, most answered "about the same", because 60% sounds equally unusual for any hospital. But a proportion computed from a small sample is more variable. Section 7.5 explains why, and Exercise 7 asks you to compute the two probabilities.

Lecture 17 — Statistics as random variables and the Law of Large Numbers

7.1 Statistics are random variables

Random samples

Before we collect data, the values we'll observe are unknown, so we model them as random variables X1,…,Xn. The simplest and most important model says that the observations are drawn independently from one fixed population.

Definition: Random Sample

Random variables X1,…,Xn form a random sample of size n from a distribution (the population distribution) if they are

  1. independent, and
  2. identically distributed: each Xi has that same distribution.

We abbreviate this as IID (independent and identically distributed). The population mean and variance are written μ=E[Xi] and σ2=Var[Xi].

Sampling at random with replacement from a finite population gives exactly an IID sample. Sampling without replacement from a population much larger than the sample (say, 1,000 voters from millions) is so close to IID that the difference doesn't matter (Section 4.6 showed the same thing for the hypergeometric and binomial).

Statistics and their sampling distributions

Definition: Statistic and Sampling Distribution

A statistic is any quantity computed from the sample: T=g(X1,…,Xn). Because the Xi are random, T is a random variable, and its distribution is called the sampling distribution of T.

The most important statistics are the sample total and the sample mean, T=X1+⋯+Xn,X‾=X1+⋯+Xnn, along with the sample proportion, median and variance S2. We keep the convention of Chapter 3: capital letters (X‾, S) for the random variables, lower-case letters (x‾, s) for the numbers computed from one particular dataset.

The sampling distribution depends on three things: the population distribution, the sample size n, and the sampling method (with or without replacement, independent or clustered). The quickest way to see one is to simulate it.

Idea: Algorithm: Simulating a Sampling Distribution

  1. Draw a sample of size n from the population.
  2. Compute the statistic.
  3. Repeat steps 1 and 2 many times (10,000 is usually plenty).
  4. The histogram of the results approximates the sampling distribution.

In R, replicate(10000, statistic(sample_from_population(n))) does steps 1 to 3 in one line.

Example 7.1: Ten Sons at a Time ★ lecture example

Treat Galton's 481 sons (Example 5.1) as a population, with mean μ=69.23 inches and SD σ=2.62 inches. (a) Draw three random samples of 10 sons (with replacement) and compute their means. (b) Simulate the sampling distribution of X‾ for samples of 10. What are its centre and spread? (c) Compare the spread with σ/10.

(a) Three samples, three answers.

R
galton <- read.csv("https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/master/csv/HistData/GaltonFamilies.csv")
sons <- galton$childHeight[galton$gender == "male"]
set.seed(1)
replicate(3, mean(sample(sons, 10, replace = TRUE)))
output
[1] 69.09 67.71 69.32

The three sample means are 69.1, 67.7 and 69.3 inches. None equals the population mean exactly, and one is off by 1.5 inches: one unlucky sample of 10 can mislead.

(b) Ten thousand samples.

R
set.seed(7)
means <- replicate(10000, mean(sample(sons, 10, replace = TRUE)))
c(mean(means), sd(means))
sigma <- sqrt(mean((sons - mean(sons))^2))      # population SD of the 481 sons
ggplot(data.frame(means), aes(means)) +
  geom_histogram(aes(y = after_stat(density)), binwidth = 0.2, colour = "white") +
  stat_function(fun = dnorm, args = list(mean = mean(sons), sd = sigma / sqrt(10)), linewidth = 1) +
  stat_function(fun = dnorm, args = list(mean = mean(sons), sd = sigma), colour = "grey60", linewidth = 1) +
  labs(x = "mean height of a sample of 10 sons (inches)", y = "density")
output
[1] 69.2290930  0.8303354
Density histogram of 10,000 sample means of 10 sons' heights drawn from Galton's 481…Density histogram of 10,000 sample means of 10 sons' heights drawn from Galton's 481 sons, centred at 69.2 inches and ranging from about 66 to 72.5. A narrow Gaussian curve with SD 0.83 fits it closely. Behind it, a much wider grey Gaussian curve with SD 2.62 shows the spread of individual heights. 0.0 0.1 0.2 0.3 0.4 0.5 67 69 71 mean height of a sample of 10 sons (inches) density

The sample means are centred at 69.23 inches, exactly the population mean, with SD 0.83 inches.

(c) σ/10=2.62/3.16≈0.83. That's no coincidence: Section 7.2 proves that the SD of a sample mean is always σ/n. The grey curve shows the spread of individual heights; averaging 10 of them shrinks the spread by a factor of 10≈3.2.

Example 7.2: Mean or Median?

In a toxicology study, the time a drug takes to be absorbed is exponential with mean 5 hours, X∼Exp(0.2). A study measures 50 patients. Simulate the sampling distributions of the sample mean and the sample median. What does each one estimate?

R
set.seed(7)
sims <- replicate(10000, { x <- rexp(50, rate = 0.2); c(mean = mean(x), median = median(x)) })
apply(sims, 1, function(v) c(centre = mean(v), spread = sd(v)))
c(population_mean = 5, population_median = 5 * log(2))
output
            mean    median
centre 5.0067945 3.5213001
spread 0.7080361 0.7036034
  population_mean population_median 
         5.000000          3.465736 
R
d <- data.frame(value = c(sims["mean", ], sims["median", ]),
                statistic = rep(c("sample mean", "sample median"), each = 10000))
ggplot(d, aes(value, fill = statistic)) +
  geom_histogram(aes(y = after_stat(density)), binwidth = 0.1, alpha = 0.6, position = "identity") +
  scale_fill_manual(values = c("steelblue", "orange")) +
  labs(x = "statistic from a sample of 50 (hours)", y = "density", fill = NULL)
Two overlapping density histograms of 10,000 simulated statistics from samples of 50…Two overlapping density histograms of 10,000 simulated statistics from samples of 50 exponential absorption times. The sample medians (orange) centre near 3.5 hours; the sample means (blue) centre near 5 hours. Both have a spread of about 0.7 hours. 0.0 0.2 0.4 0.6 2 4 6 8 statistic from a sample of 50 (hours) density sample mean sample median

Both sampling distributions have about the same spread (0.71 and 0.70 hours), but they are centred in different places. The sample mean centres on the population mean, 5 hours. The sample median centres near the population median, 5ln2≈3.47 hours (Section 5.4), slightly above it at 3.52. For a skewed population these are different targets, and neither is "wrong": they answer different questions. The mean predicts the total drug exposure across many patients; the median describes a typical patient.

Watch Out: When the IID Model Fails

Two common situations break the random-sample model, and the formulas of this chapter then give wrong answers.

  • Time series. A patient's weight measured monthly during a weight-loss trial isn't IID: each month depends on the last (dependence), and the mean drifts down if the drug works (not identically distributed). Daily ozone readings (Chapter 1) behave the same way.
  • Clusters. Galton's sons came in families: brothers share genes and dinners, so their heights are correlated. A survey that interviews everyone in 50 households has far less information than one that interviews 1 person in each of 250 households, though both have 250 respondents.

Example 6.13 and Exercise 13 show what positive correlation does to an average: it erodes the benefit of a larger sample.

Check Yourself

A statistic is computed from a random sample. Which statement is true?

Think About It: Who Is in the Sample?

A university wants the average commuting time of its students. Option A: stand at the main entrance at 9 a.m. and ask the first 200 students who walk in. Option B: email 200 students chosen at random from the registrar's list. Which is closer to a random sample from the population of students? What could go wrong with each?

Option B is much closer, if most of the chosen students reply. Option A samples only students on campus at 9 a.m. on that day and at that entrance: part-time students, evening students and those with long commutes who arrive later are under-represented, and friends arriving together aren't independent. Option B's danger is non-response: if students with long commutes are too busy to reply, the answers are biased again. Neither problem is fixed by a bigger sample. The formulas in this chapter describe random sampling error only; Section 7.5's warning box returns to this.

Key Takeaways

  • A random sample is IID: independent draws from one population with mean μ and variance σ2.
  • A statistic computed from a random sample is a random variable; its distribution is the sampling distribution.
  • Simulation draws the sampling distribution: sample, compute, repeat, histogram.
  • Different statistics have different sampling distributions and may estimate different things.

7.2 The sample mean and total

For the sample mean and total, Chapter 6's rules give the centre and spread of the sampling distribution exactly, whatever the population.

Theorem: Mean and Variance of the Sample Mean and Total

If X1,…,Xn is a random sample from a population with mean μ and variance σ2, then E[T]=nμ,Var[T]=nσ2,SD[T]=σn, E[X‾]=μ,Var[X‾]=σ2n,SD[X‾]=σn.

By linearity of expectation (Section 6.4), E[T]=E[X1]+⋯+E[Xn]=nμ. The Xi are independent, so all covariances are zero and variances add (Section 6.6): Var[T]=nσ2. Then X‾=T/n, so E[X‾]=nμ/n=μ and Var[X‾]=nσ2/n2=σ2/n. ◼∎

Two features of this result run through the rest of the course.

The picture below shows the square-root law for Galton's sons, modelled as N(69.2,2.62). Each curve is the sampling distribution of the mean of n sons. Every fourfold increase in n halves the width.

R
x <- seq(61, 77.5, by = 0.02)
d <- do.call(rbind, lapply(c(1, 4, 16, 64), function(n)
  data.frame(n = paste("n =", n), x, f = dnorm(x, 69.2, 2.6 / sqrt(n)))))
d$n <- factor(d$n, levels = paste("n =", c(1, 4, 16, 64)))
peaks <- data.frame(n = levels(d$n), x = 69.2, f = dnorm(0, 0, 2.6 / sqrt(c(1, 4, 16, 64))))
ggplot(d, aes(x, f, group = n)) + geom_line(linewidth = 0.9) +
  geom_text(data = peaks, aes(label = n), x = 69.2 + c(4.3, 2.4, 1.4, 0.9), hjust = 0, vjust = 0.2, size = 3.8) +
  labs(x = "mean height of n sons (inches)", y = "density")
Four Gaussian density curves, all centred at 69.2 inches, for the mean height of samples…Four Gaussian density curves, all centred at 69.2 inches, for the mean height of samples of 1, 4, 16 and 64 sons. The curve for n = 1 is wide and low, with SD 2.6 inches; each fourfold increase in n halves the SD (1.3, 0.65, 0.325), so the curves become successively narrower and taller. n = 1 n = 4 n = 16 n = 64 0.00 0.25 0.50 0.75 1.00 1.25 65 70 75 mean height of n sons (inches) density

Definition: Standard Error

The SD of a statistic's sampling distribution is called its standard error (SE). For the sample mean, SE[X‾]=σ/n. When σ is unknown we estimate it by the sample SD s, giving the estimated standard error s/n.

The name distinguishes the two kinds of spread: the standard deviation σ describes how much individuals vary; the standard error describes how much a statistic varies from sample to sample.

If the population is Gaussian, we know more than the mean and variance: the whole sampling distribution.

Theorem: Sampling from a Gaussian Population

If X1,…,Xn is a random sample from N(μ,σ2), then exactly T∼N(nμ,nσ2),X‾∼N(μ,σ2n).

Sums of independent Gaussian random variables are Gaussian (Section 6.6), with means and variances adding; dividing by n rescales a Gaussian to another Gaussian. ◼∎

Example 7.3: How Many People Fit in the Elevator? ★ lecture example

An elevator is rated for 2,500 lb. Suppose adult passengers' weights have mean 185 lb (the figure the US Coast Guard has used for passenger vessels since 2011) and SD 40 lb, and model them as roughly Gaussian and independent. (a) If 12 people get in, what is the probability that they exceed 2,500 lb? (b) What is the largest number of passengers for which this probability is below 1 in 1,000?

(a) The total of 12 weights has mean 12×185=2,220 lb and SD 4012≈138.6 lb. So Pr(T>2500)≈Pr(Z>2500−2220138.6)=Pr(Z>2.02)≈0.022. About 1 load in 45. Notice the SD of the total: 12 people vary by 40 lb each, but the total varies by only 139 lb, not 12×40=480, because heavy and light passengers partly cancel.

(b) Try neighbouring values of n:

R
n <- 10:14
data.frame(n, mean = 185 * n, sd = round(40 * sqrt(n), 1),
           p_over = signif(1 - pnorm(2500, mean = 185 * n, sd = 40 * sqrt(n)), 2))
output
   n mean    sd  p_over
1 10 1850 126.5 1.4e-07
2 11 2035 132.7 2.3e-04
3 12 2220 138.6 2.2e-02
4 13 2405 144.2 2.6e-01
5 14 2590 149.7 7.3e-01

With 11 people the probability is about 0.0002; with 12 it jumps to 0.022, and with 13 to 0.26. The safe limit is 11 passengers. Most of the risk comes from the mean, which grows by 185 lb with each extra person, while the SD grows only slowly.

Why the average matters. In October 2005 the tour boat Ethan Allen capsized on Lake George, New York, and 20 passengers died. Investigators found that the boat was carrying too much weight, though not too many people: its capacity had been set decades earlier, using an assumption of 140 lb per passenger, and was never rechecked after a heavy canopy was added to the boat. In 2011 the Coast Guard raised its standard assumption to 185 lb per passenger. A model of a total is only as good as its mean.

Example 7.4: One Worker or Nine?

At a battery plant, workers' blood lead levels are N(25,52) µg/dL. (a) What is the probability that one randomly chosen worker is above 30? (b) For a random sample of 9 workers, what is the probability that their average is above 30? (c) What is the probability that their total is above 270?

(a) Pr(X>30)=Pr(Z>1)≈0.159.

(b) X‾∼N(25,52/9), with SE 5/3≈1.67: Pr(X‾>30)=Pr(Z>30−255/3)=Pr(Z>3)≈0.0013.

(c) The total exceeds 270 exactly when the average exceeds 270/9=30, so the answer is the same, 0.0013. Directly: T∼N(225,9×25), with SD 15, and Pr(Z>45/15=3).

R
1 - pnorm(30, mean = 25, sd = 5)
1 - pnorm(30, mean = 25, sd = 5 / 3)
1 - pnorm(270, mean = 225, sd = 15)
output
[1] 0.1586553
[1] 0.001349898
[1] 0.001349898

One worker in six is above 30, but a random group of nine averaging above 30 would be extraordinary. If it happened, we'd suspect the group wasn't random: perhaps they all work in the same high-exposure area. That's the logic of the tests in Chapter 9.

Check Yourself

A population has SD σ=4. What is the standard error of the mean of a random sample of 64?

A lab wants to make its average measurement 3 times as precise (SE one-third as big). How many times as many measurements does it need?

Think About It: Averages and Individuals

A newspaper reports that "the average American woman is 5 ft 4 in and the average man 5 ft 9 in, so a woman is shorter than a man." Is that a statement about individuals or about averages? If you pick one woman and one man at random, how likely is the woman to be taller? (Assume SDs of about 3 inches for each and Gaussian heights.)

It's a statement about averages, and it doesn't hold for every pair. The difference W−M between a random woman's and man's heights has mean 64−69=−5 inches and, by Section 6.6, SD 32+32≈4.24. So Pr(W>M)=Pr(Z>54.24)=Pr(Z>1.18)≈0.12: about one random pair in eight has a taller woman. Comparing averages is very different from comparing individuals, a point worth remembering whenever group averages are used to judge a person.

Key Takeaways

  • E[X‾]=μ and SE[X‾]=σ/n; E[T]=nμ and SD[T]=σn.
  • The standard error measures sample-to-sample variation of a statistic; it shrinks like 1/n.
  • From a Gaussian population, X‾ and T are exactly Gaussian.

7.3 The Law of Large Numbers

Chapter 3 promised a proof that a long-run proportion settles down to the probability. We now prove it for any average, using two simple inequalities that need nothing but a mean and a variance.

Theorem: Markov's and Chebyshev's Inequalities

  1. Markov. If X≥0 and a>0, then Pr(X≥a)≤E[X]a.
  2. Chebyshev. For any random variable with mean μ and SD σ, and any k>0, Pr(|X−μ|≥kσ)≤1k2.

(1) Let I=1 if X≥a and 0 otherwise. Then aI≤X always: if I=1 then X≥a, and if I=0 then X≥0. Taking expectations, aPr(X≥a)≤E[X].

(2) Apply Markov to the non-negative variable (X−μ)2 with a=k2σ2: Pr(|X−μ|≥kσ)=Pr((X−μ)2≥k2σ2)≤E[(X−μ)2]k2σ2=1k2.◼∎

Chebyshev's inequality is the probability version of the rule for data in Section 2.3: at least 75% of any distribution lies within 2 SDs of the mean, and at least 89% within 3. It's crude (for a Gaussian the true figures are 95% and 99.7%), but it holds for every distribution with a finite variance.

Theorem: The (Weak) Law of Large Numbers

Let X‾n be the mean of a random sample of size n from a population with mean μ and finite variance σ2. For any margin ε>0, however small, Pr(|X‾n−μ|≥ε)≤σ2nε2→0as n→∞.

X‾n has mean μ and variance σ2/n (Section 7.2). Chebyshev with kσ replaced by ε gives Pr(|X‾n−μ|≥ε)≤σ2/nε2, which goes to 0 as n grows. ◼∎

Applied to indicators, it says that the proportion of times an event happens converges to its probability: the long-run frequency interpretation of Section 3.1, now proved rather than assumed.

Example 7.5: A Coin in a Prison Camp ★ lecture example

When Germany invaded Denmark in April 1940, the South African mathematician John Kerrich was visiting Copenhagen, and he spent the rest of the war in an internment camp. With a fellow internee, Eric Christensen, he tossed a coin 10,000 times and got 5,067 heads (published in 1946). (a) What is the standard error of the proportion of heads in 10,000 tosses of a fair coin? (b) How many standard errors from 0.5 is Kerrich's 0.5067? (c) What does Chebyshev's inequality say about a deviation this big, and what does the Gaussian approximation say?

(a) Each toss is an indicator with σ2=p(1−p)=0.25, so the proportion p̂ has SE[p̂]=0.2510,000=0.5100=0.005.

(b) 0.5067−0.50.005=1.34 standard errors.

(c) Chebyshev with k=1.34: Pr(|p̂−0.5|≥0.0067)≤11.342≈0.56. That's a true but weak bound. The Gaussian approximation of Section 7.4 gives Pr(|Z|≥1.34)≈0.18: a deviation of this size or more would happen in about 1 in 5 such experiments, so Kerrich's coin looks fair.

R
set.seed(1940)
n <- 1:10000
heads <- cumsum(rbinom(10000, 1, 0.5))
d <- data.frame(n, prop = heads / n, lo = 0.5 - 2 * sqrt(0.25 / n), hi = 0.5 + 2 * sqrt(0.25 / n))
d <- d[unique(round(10^seq(0, 4, length.out = 600))), ]   # thin the points for plotting
ggplot(d, aes(n, prop)) +
  geom_line(aes(y = lo), linetype = "dashed") + geom_line(aes(y = hi), linetype = "dashed") +
  geom_line(colour = "steelblue") +
  scale_x_log10() + coord_cartesian(ylim = c(0, 1)) +
  labs(x = "number of tosses (log scale)", y = "proportion of heads")
Line plot of the running proportion of heads in 10,000 simulated tosses of a fair coin…Line plot of the running proportion of heads in 10,000 simulated tosses of a fair coin, against the number of tosses on a logarithmic axis. The proportion swings widely in the first few dozen tosses and then settles close to 0.5. Dashed curves at 0.5 plus and minus two standard errors form a funnel that narrows as the number of tosses grows, and the running proportion stays mostly inside it. 0.00 0.25 0.50 0.75 1.00 1 10 100 1000 10000 number of tosses (log scale) proportion of heads

The funnel narrows like 1/n: with 100 tosses, the proportion is usually within 0.10 of a half; with 10,000, within 0.01.

Watch Out: The Law of Averages Doesn't Balance the Books

The Law of Large Numbers says the proportion of heads approaches 12. It does not say that the number of heads catches up with the number of tails. After 10,000 tosses Kerrich had 134 more heads than tails, and the expected gap between heads and tails actually grows (like n). The proportion settles because each early surplus is swamped by an ever larger number of tosses, not because the coin compensates. Believing that a run of tails makes heads "due" is the gambler's fallacy: a coin has no memory.

Check Yourself

Use Chebyshev's inequality to bound Pr(|X−μ|≥2σ) for any distribution.

Think About It: Black, Black, Black...

On 18 August 1913, at the casino in Monte Carlo, the ball on a roulette wheel landed on black 26 times in a row. Gamblers lost millions betting on red, reasoning that red was "due". (a) Were they right? (b) Was 26 blacks in a row too unlikely to happen anywhere, ever?

(a) No. Each spin is independent; after 25 blacks, the chance of red on the next spin is still 1837 (a European wheel has 18 red, 18 black and one green zero). The Law of Large Numbers works by swamping, not by compensation. (b) The chance that a particular run of 26 spins is all black is (1837)26≈1 in 137 million. That's tiny for one sequence of spins, but casinos around the world spin roulette wheels billions of times, and long runs are bound to happen somewhere (the same logic as Paul the octopus in Section 4.4 and the cancer clusters in Section 4.7). The story is famous because it was selected for being extreme.

Key Takeaways

  • Chebyshev: Pr(|X−μ|≥kσ)≤1/k2 for any distribution with finite variance.
  • Law of Large Numbers: X‾n gets arbitrarily close to μ with probability approaching 1. Proportions converge to probabilities.
  • It works by swamping, not compensation: the gambler's fallacy is a fallacy.

Lecture 18 — The Central Limit Theorem and the bootstrap

7.4 The Central Limit Theorem

The Law of Large Numbers says X‾ settles near μ. The Central Limit Theorem says how it's distributed around μ on the way: approximately Gaussian, whatever the shape of the population.

Theorem: The Central Limit Theorem

Let X1,X2,… be IID with mean μ and finite variance σ2>0. Then for every z, Pr(X‾n−μσ/n≤z)→Φ(z)as n→∞. In practice: for large n, X‾≈N(μ,σ2n),T≈N(nμ,nσ2).

A sketch, using moment generating functions, for populations that have one. Standardize each observation, Yi=(Xi−μ)/σ, so that E[Yi]=0 and E[Yi2]=1, and let Zn=(Y1+⋯+Yn)/n, which is the standardized X‾n. By independence (Section 6.6), the MGF of Zn is a product: MZn(t)=[MY(tn)]n. A Taylor expansion of MY near 0 gives MY(s)=1+sE[Y]+s22E[Y2]+⋯=1+s22+⋯, so MZn(t)=[1+t22n+⋯]n→et2/2, using (1+a/n)n→ea. And et2/2 is the MGF of N(0,1). A theorem from advanced probability says that convergence of MGFs implies convergence of distributions. ◼∎

The proof shows why the population's shape washes out: only its mean and variance survive the limit. Skewness and every other feature are divided by higher powers of n and fade away.

The sampling distribution, built one sample at a time

Pick a population and a sample size n. Each time you draw a sample, its n values appear as dots in the middle panel, and the statistic computed from them drops into the histogram at the bottom. Draw thousands and the histogram becomes the sampling distribution. For the mean and the total, the orange curve is the Gaussian approximation from the Central Limit Theorem. Try the skewed and two-humped populations with n = 2, 5, 30.

TopTop: the population distribution. Middle: the latest sample as dots with its statistic marked. Bottom: histogram of the statistic over all samples drawn, with the Gaussian approximation.

Example 7.6: Averaging River Lengths ★ lecture example

The lengths of the 141 longest North American rivers (rivers, Example 5.12) are strongly right-skewed, with mean 591 miles and SD 494 miles. Treat them as a population, and simulate the sampling distribution of X‾ for samples of size 1, 5 and 30, drawn with replacement. Compare each with the Central Limit Theorem's Gaussian approximation.

R
set.seed(7)
sizes <- c(1, 5, 30)
d <- do.call(rbind, lapply(sizes, function(n)
  data.frame(n = paste("n =", n), xbar = replicate(10000, mean(sample(rivers, n, replace = TRUE))))))
d$n <- factor(d$n, levels = paste("n =", sizes))
mu <- mean(rivers); sigma <- sqrt(mean((rivers - mu)^2))
curves <- do.call(rbind, lapply(sizes, function(n) {
  x <- seq(max(0, mu - 4 * sigma / sqrt(n)), mu + 4 * sigma / sqrt(n), length.out = 300)
  data.frame(n = paste("n =", n), x, y = dnorm(x, mu, sigma / sqrt(n))) }))
curves$n <- factor(curves$n, levels = paste("n =", sizes))
ggplot(d, aes(xbar)) +
  geom_histogram(aes(y = after_stat(density)), bins = 40, colour = "white") +
  geom_line(data = curves, aes(x, y), linewidth = 0.9) +
  facet_wrap(~ n, scales = "free") + scale_x_continuous(n.breaks = 3) +
  labs(x = "sample mean of river lengths (miles)", y = "density")
Three density histograms of simulated sample means of river lengths, for samples of size…Three density histograms of simulated sample means of river lengths, for samples of size 1, 5 and 30, each with the Gaussian curve predicted by the Central Limit Theorem. For n = 1, the histogram is the strongly right-skewed population, peaking near 300 miles with a tail beyond 3,000, and the Gaussian curve fits badly. For n = 5, the histogram is narrower and still right-skewed. For n = 30, it is nearly symmetric, centred at 591 with SD about 90, and the Gaussian curve fits well. n = 1 n = 5 n = 30 0 2000 0 1000 500 1000 0.000 0.001 0.002 0.003 0.004 0.000 0.001 0.002 0.0000 0.0005 0.0010 0.0015 0.0020 0.0025 sample mean of river lengths (miles) density

For n=1 the "sample mean" is just one river, so its distribution is the skewed population itself, and a Gaussian curve is a poor description. By n=5 the distribution is narrower but still visibly skewed. By n=30 it's nearly symmetric, and the Gaussian curve with mean 591 and SE 494/30≈90 fits well. A strongly skewed population needs a larger sample before the approximation is good; a symmetric one needs very little.

How large is "large n"? There's no universal answer, because it depends on the population.

For proportions, a common check is that both np and n(1−p) should be at least 10.

Example 7.7: Why the Casino Always Wins ★ lecture example

An American roulette wheel has 18 red, 18 black and 2 green pockets. A $1 bet on red wins $1 with probability 1838 and loses $1 otherwise. (a) Find the mean and SD of the gambler's winnings on one bet. (b) A gambler makes 1,000 bets. What is the probability that they finish ahead? (c) Over a busy month the casino takes 1,000,000 such bets. Predict its profit.

(a) One bet. Let W=+1 with probability 1838 and −1 with probability 2038: E[W]=18−2038=−119≈−0.0526,Var[W]=E[W2]−E[W]2=1−1361≈0.997, so SD[W]≈0.999. The gambler loses about 5.3 cents per dollar on average, but each bet is mostly noise.

(b) 1,000 bets. The total has mean 1000×(−0.0526)=−52.6 dollars and SD 0.9991000≈31.6 dollars, so by the CLT Pr(T>0)≈Pr(Z>0−(−52.6)31.6)=Pr(Z>1.66)≈0.048. About 1 gambler in 20 is ahead after 1,000 bets. (The exact binomial answer, Pr(more than 500 wins), is 0.045.)

(c) 1,000,000 bets. The casino's profit −T has mean $52,632 and SD 0.999×1000≈$999. With 95% probability it lies within about 2 SDs: between roughly $50,600 and $54,600.

R
p <- 18 / 38; m <- 2 * p - 1; s <- sqrt(1 - m^2)
c(mean = m, sd = s)
1 - pnorm(0, mean = 1000 * m, sd = s * sqrt(1000))      # (b) CLT
1 - pbinom(500, 1000, p)                                 # (b) exact
-1e6 * m + c(-2, 2) * s * 1000                           # (c) casino's 95% range
output
       mean          sd 
-0.05263158  0.99861400 
[1] 0.04779035
[1] 0.04479591
[1] 50634.35 54628.81

The pattern. The mean of the total grows like n, its SD only like n. For one gambler, the small negative mean is buried in noise. For the casino, the noise is tiny compared with the mean, and the profit is as predictable as a salary. Insurance companies, lotteries and blood banks run on the same arithmetic.

History: From De Moivre to the 'Central' Limit Theorem

Abraham de Moivre, a French Protestant who fled to London to escape religious persecution and made a living tutoring mathematics and advising gamblers, proved the first version of the theorem in 1733: the binomial distribution, for large n, is approximately the bell curve (Section 5.3). Pierre-Simon Laplace extended it to sums of many kinds of random errors around 1810, and the Russian mathematician Aleksandr Lyapunov gave the first fully rigorous general proof in 1901. The name came later still: in 1920 George Pólya called it the zentraler Grenzwertsatz, the limit theorem that is central to probability.

Francis Galton, who built a board of pins (the "quincunx", now called a Galton board) to show balls piling up into a bell curve, wrote in 1889 of the bell-shaped "law of frequency of error" that "the law would have been personified by the Greeks and deified, if they had known of it." As Section 5.3 describes, Galton used the bell curve to rank people for eugenic purposes, a use the mathematics never justified. The theorem itself is about averages of independent random quantities, and it says nothing about the worth of the individuals being averaged.

Check Yourself

A population has mean 50 and SD 12 (shape unknown). Use the CLT to find the probability that the mean of a sample of 36 exceeds 54.

Which of these would need the largest sample before X‾ is close to Gaussian?

Think About It: When Can't You Average Your Way Out?

A city's fire department plans its budget around the average annual cost of fires, estimated from 30 years of data. In most years the cost is close to average, but once every few decades a huge fire costs as much as 50 normal years. Why might the CLT give the department false confidence? What would you advise?

The CLT's approximation is good when no single observation can dominate the total. With rare, enormous events, 30 years may include none of them (so the sample mean underestimates the true mean) or one of them (so the sample mean swings wildly). The "typical error" s/30 is then itself unreliable, because s depends heavily on whether a disaster happened to occur. The department should model the rare disasters separately, perhaps with reinsurance or a reserve fund, rather than treating the annual cost as approximately Gaussian. Catastrophe insurers, earthquake engineers and pandemic planners all face this problem.

Key Takeaways

  • Central Limit Theorem: for large n, X‾≈N(μ,σ2/n) and T≈N(nμ,nσ2), whatever the population (if its variance is finite).
  • How large n must be depends on the population: small for symmetric, larger for skewed, never for infinite variance.
  • The mean of a total grows like n, its SD like n: totals of many random quantities are predictable.

7.5 Gaussian approximations at work

Sample proportions

A proportion is a mean of indicators: if each Xi is 1 with probability p and 0 otherwise, the sample proportion p̂=X‾. Each indicator has mean p and variance p(1−p), so Section 7.2 and the CLT give E[p̂]=p,SE[p̂]=p(1−p)n,p̂≈N(p,p(1−p)n). This is the Gaussian approximation to the binomial (Section 5.3), divided by n.

Example 7.8: The Margin of Error in a Poll ★ lecture example

A poll asks 1,000 randomly chosen voters whether they support a proposal. (a) If the true level of support is 50%, what is the standard error of the poll's proportion? (b) Within what range around 50% will the poll fall 95% of the time? (c) How many voters would be needed to halve that range?

(a) SE=0.5×0.51000≈0.0158, about 1.6 percentage points.

(b) 95% of a Gaussian lies within 1.96 SDs of its mean: 0.5±1.96×0.0158=0.5±0.031, so between about 46.9% and 53.1%. This "±3 points" is the margin of error that news reports quote for polls of about 1,000 people. (It's largest at p=0.5, so pollsters quote that worst case.)

(c) Halving the margin needs four times as many voters: 4,000.

R
se <- sqrt(0.5 * 0.5 / 1000)
c(se = se, margin = qnorm(0.975) * se)
set.seed(7)
polls <- rbinom(10000, 1000, 0.5) / 1000            # 10,000 simulated polls
mean(abs(polls - 0.5) <= 0.031)                      # close to 0.95
output
        se     margin 
0.01581139 0.03098975 
[1] 0.9447

The surprising part. The margin of error depends on the sample size, not on the size of the population. A random sample of 1,000 is as precise for the whole United States as it is for a city of 100,000. What matters is how many people you ask and whether they were chosen at random, which is where real polls run into trouble (see the warning box below).

The most dangerous equation

The statistician Howard Wainer called SE=σ/n "the most dangerous equation", not because it's hard but because ignoring it has cost so much.

Example 7.9: The Kidney Cancer Map ★ lecture example

Wainer described maps of US kidney cancer death rates by county in the 1980s: the counties with the lowest rates and those with the highest rates were both mostly small, rural counties. Simulate 3,000 counties, with populations varying widely (median about 30,000), in which everyone has the same true death rate of 5 per 100,000 per year. Count deaths over 10 years. Which counties have the most extreme observed rates?

R
set.seed(7)
pop <- round(exp(rnorm(3000, log(30000), 1.2)))      # county populations
deaths <- rpois(3000, pop * 5e-5 * 10)                 # 10 years at 5 per 100,000 per year
rate <- deaths / (pop * 10) * 1e5
o <- order(rate)
group <- rep("middle", 3000); group[o[1:300]] <- "lowest 10%"; group[o[2701:3000]] <- "highest 10%"
c(all = median(pop), lowest = median(pop[group == "lowest 10%"]), highest = median(pop[group == "highest 10%"]))
ggplot(data.frame(pop, rate, group), aes(pop, rate, colour = group)) +
  geom_point(alpha = 0.6, size = 1.2) + geom_hline(yintercept = 5) +
  scale_x_log10(labels = scales::label_comma()) +
  scale_colour_manual(values = c("orange", "steelblue", "grey60")) +
  coord_cartesian(ylim = c(0, 25)) +
  labs(x = "county population (log scale)", y = "deaths per 100,000", colour = NULL)
output
    all  lowest highest 
  30099   10632    9956 
Scatter plot of simulated county death rates per 100,000 per year against county…Scatter plot of simulated county death rates per 100,000 per year against county population on a logarithmic axis, for 3,000 counties with the same true rate of 5. The points form a funnel: small counties scatter from 0 to over 20, large counties cluster tightly around 5. The 10% lowest rates (blue) and 10% highest rates (orange) are almost all in the small counties on the left. 0 5 10 15 20 25 1,000 10,000 100,000 1,000,000 county population (log scale) deaths per 100,000 highest 10% lowest 10% middle

The median county has about 30,000 people, but the median of the 300 lowest-rate counties is about 10,600, and of the 300 highest-rate counties about 10,000. Both extremes are dominated by small counties, although every county has the same true rate. A county of 2,000 people expects only 1 death in 10 years: zero deaths gives the "safest" rate in the country, and three deaths one of the "most dangerous". The SE of a rate is proportional to 1/population, so small places produce extremes in both directions.

The real cost. Around 2000, the Gates Foundation and others spent well over a billion dollars converting large high schools into small ones, partly because small schools were overrepresented among the best-performing schools. Wainer pointed out that they were also overrepresented among the worst: exactly what the equation predicts.

Source: H. Wainer (2007), "The most dangerous equation", American Scientist 95, 249–256.

The most dangerous equation: small places, extreme rates

Each dot is a simulated county. Every county has exactly the same true death rate, 5 per 100,000 people per year, and we count deaths over 10 years. Nothing distinguishes the counties except their size. The 10% of counties with the lowest observed rates and the 10% with the highest are coloured. The dashed curves are the true rate ± 2 standard errors, which shrink like 1/√(population). Press new counties a few times.

Scatter plot of observed death rate against county population on a logarithmic axis, with…Scatter plot of observed death rate against county population on a logarithmic axis, with a funnel of plus or minus two standard errors around the true rate. The lowest and highest rates are coloured and mostly belong to small counties.

History: The Trial of the Pyx

Since the Middle Ages, the Royal Mint in London has checked its coins in a ceremony called the Trial of the Pyx (the pyx is the box in which sample coins are kept). Coins are taken from production, and their total weight is compared with the standard. The Mint was allowed a tolerance for each coin, and for a sample of n coins the allowed tolerance on the total was n times as much. But totals vary like n, not like n (Section 7.2). With a sample of 100 coins, the tolerance was 10 times wider than honest random variation required. For centuries, a master of the Mint could make coins slightly light on purpose, keep the difference in gold, and still pass the Trial, and some masters were suspected of doing just that. The error persisted long after de Moivre's work could have corrected it. Wainer tells this story in the article cited in Example 7.9.

More Gaussian approximations

Example 7.10: Arsenic at a Contaminated Site

Arsenic concentrations in soil samples from a contaminated site are right-skewed, with mean 18 ppm and SD 6 ppm. (a) Can we use a Gaussian model for the concentration of a single sample? (b) For 36 independent samples, approximate the probability that their average exceeds 19.5 ppm. (c) What about their total exceeding 702 ppm?

(a) No: the population is skewed, so a Gaussian model for one sample would misstate the tail probabilities. We'd need the actual distribution.

(b) With n=36, the CLT gives X‾≈N(18,6236), with SE 1: Pr(X‾>19.5)≈Pr(Z>1.5)≈0.067.

(c) 702/36=19.5, so it's the same event and the same probability, 0.067.

Example 7.11: Is This Die Loaded? ★ lecture example

A die is rolled 600 times and shows 120 sixes. If the die is fair, how unusual is a count this high? Use the Gaussian approximation with continuity correction, and compare with the exact binomial.

If the die is fair, the number of sixes X∼Bin(600,16) has mean 100 and SD 600×16×56≈9.13. With the continuity correction (Section 5.3): Pr(X≥120)≈Pr(Z≥119.5−1009.13)=Pr(Z≥2.14)≈0.016.

R
1 - pnorm(119.5, mean = 100, sd = sqrt(600 * 5 / 36))
1 - pbinom(119, 600, 1 / 6)            # exact: 0.018
output
[1] 0.0163349
[1] 0.01801076

A fair die would show 120 or more sixes only about 2% of the time, so we have some evidence that this die favours six. Chapter 9 turns this reasoning into a formal test, and asks how surprised we should be when we've tested many dice.

Watch Out: A Bigger Sample Doesn't Fix a Biased One

The standard error measures only random sampling error: the scatter that comes from which units happened to be chosen. It says nothing about bias, the error that comes from how they were chosen. If people with a particular opinion are more likely to answer a poll, the poll is biased, and asking more people just gives a more precise wrong answer. The Literary Digest collected 2.4 million ballots in 1936 and called the election for the loser (Section 1.3). In 2018 the statistician Xiao-Li Meng showed that a self-selected sample of 2.3 million US voters (about 1% of the electorate), with the tiny correlation between candidate choice and willingness to respond that he estimated from 2016 survey data, would carry about as much information as a simple random sample of only 400 people. With bias, the error doesn't shrink like 1/n; it doesn't shrink at all.

Bayesian Inference: Posteriors Become Gaussian Too

The CLT has a Bayesian cousin. In Section 5.5, a free-throw shooter who made 7 of 10 shots had posterior Beta(8,4), a lopsided curve. After 70 makes in 100 shots, starting from the same uniform prior, the posterior is Beta(71,31), and it is almost exactly Gaussian, with mean 71102≈0.696 and SD 0.696×0.304/103≈0.045.

R
qbeta(c(0.025, 0.975), 71, 31)                          # exact posterior: 0.604 to 0.781
71/102 + c(-1, 1) * 1.96 * sqrt(71/102 * 31/102 / 103)  # Gaussian approximation: 0.607 to 0.785
output
[1] 0.6038529 0.7810213
[1] 0.6072510 0.7849059

With plenty of data, the posterior concentrates around the observed proportion, its spread shrinks like 1/n, and its shape becomes Gaussian, whatever the (reasonable) prior. That's why Bayesian and classical answers usually agree for large samples, and differ most when data are scarce.

Check Yourself

A poll of 1,000 voters finds 50% support. Find the 95% margin of error, 1.960.25/1000, to 4 decimal places.

A poll of 1,000 people from a city of 100,000 and a poll of 1,000 people from a country of 300 million, both random samples. Which has the smaller margin of error?

Think About It: Ranking Hospitals

A newspaper ranks 200 hospitals by their death rate after heart surgery. The "worst" hospital performed 12 operations last year, with 3 deaths (25%); the national average is 4%. Should patients avoid it? What would you want to know, and how would you present the rankings more fairly?

Not on this evidence alone. With 12 operations and a true rate of 4%, the number of deaths is roughly Bin(12,0.04), and 3 or more deaths happens about 1% of the time: unusual for one hospital, but among 200 hospitals, many of them small, we'd expect a couple to look this bad by chance. You'd want several years of data, the hospital's caseload (does it take the riskiest patients?), and the rates adjusted for risk. A fairer display is a funnel plot: each hospital's rate against its number of operations, with control limits that narrow as volume grows, like the widget above. Hospitals outside the funnel deserve a closer look. The statistician David Spiegelhalter, who led the statistical work for the inquiry into children's heart surgery at the Bristol Royal Infirmary in the 1990s, championed funnel plots for exactly this purpose, and they are now widely used to publish hospital and surgeon outcomes in the UK.

Key Takeaways

  • A sample proportion has SE=p(1−p)/n; a poll of 1,000 has a margin of error of about ±3 points, whatever the population size.
  • Small samples produce extreme averages in both directions: beware of rankings and maps of rates from units of very different sizes.
  • The SE measures random error only; a biased sample stays biased however large it is.

7.6 Simulation as a microscope: the bootstrap

For the mean, the total and the proportion, formulas give the standard error. For most other statistics (the median, a trimmed mean, a correlation, a ratio) there's no simple formula. Simulation still works, if we know the population: draw samples, compute the statistic, look at the spread. But in real problems we don't know the population. All we have is one sample.

In 1979 the Stanford statistician Bradley Efron proposed a beautifully simple idea: use the sample as a stand-in for the population. If the sample is random, it's a miniature version of the population, so drawing from it (with replacement, so that each draw is IID) mimics drawing new samples from the population. He called it the bootstrap, after the impossible feat of pulling yourself up by your own bootstraps.

Idea: Algorithm: The Bootstrap Standard Error

Given data x1,…,xn and a statistic T:

  1. Draw a resample of size n from x1,…,xn with replacement.
  2. Compute the statistic on the resample, T*.
  3. Repeat 1 and 2 many times (say 2,000).
  4. The SD of the T* values estimates the standard error of T.

Example 7.12: How Precise Is the Median River? ★ lecture example

The median of the 141 river lengths is 425 miles. (a) Use the bootstrap to estimate its standard error. (b) Check the method on the sample mean, where we know the formula s/n.

(a) Resample, recompute, repeat.

R
set.seed(7)
boot_medians <- replicate(2000, median(sample(rivers, replace = TRUE)))
c(median = median(rivers), boot_se = sd(boot_medians))
quantile(boot_medians, c(0.025, 0.975))
output
   median   boot_se 
425.00000  26.66176 
 2.5% 97.5% 
  380   500 
R
ggplot(data.frame(m = boot_medians), aes(m)) +
  geom_histogram(binwidth = 5, colour = "white") +
  geom_vline(xintercept = median(rivers)) +
  labs(x = "median of a bootstrap resample (miles)", y = "count")
Histogram of 2,000 bootstrap medians of the river lengths, mostly between about 375 and…Histogram of 2,000 bootstrap medians of the river lengths, mostly between about 375 and 505 miles. The values pile up on a few particular numbers, such as 420, 425 to 435 and 444, because a median is always one of the data values or the average of two. The original median, 425, is marked with a vertical line. 0 100 200 350 400 450 500 550 median of a bootstrap resample (miles) count

The bootstrap SE of the median is about 27 miles, and the middle 95% of the bootstrap medians runs from about 380 to 500 miles. Notice the lumpy histogram: a median can only take one of the data values (or the average of two), so the resampled medians pile up on a few numbers.

(b) Check on the mean.

R
set.seed(7)
boot_means <- replicate(2000, mean(sample(rivers, replace = TRUE)))
c(bootstrap = sd(boot_means), formula = sd(rivers) / sqrt(length(rivers)))
output
bootstrap   formula 
 41.23629  41.59143 

The bootstrap gives 41.2 miles and the formula 41.6: they agree, as they should. The bootstrap needs no formula, so it works for the median too. R's boot package automates all of this: boot::boot(rivers, function(x, i) median(x[i]), R = 2000).

Prediction and inference. The median river is 425 miles long in this sample, but if the almanac had listed a different random set of rivers, the median would typically differ by about 27 miles. That's the inference question every estimate must answer, and Chapter 8 turns standard errors like this one into confidence intervals.

The bootstrap: resampling the river data

The top strip shows the 141 river lengths (R's rivers), with their median, 425 miles, marked. We can't go back and measure more rivers, so we treat the sample as a stand-in for the population: draw 141 values from it with replacement, so some rivers appear twice and others not at all, and compute the median of this resample. Each resample's median drops into the histogram. Its spread estimates how much the median would vary from sample to sample: the standard error.

TopTop: strip of the 141 river lengths with the latest resample highlighted. Bottom: histogram of the bootstrap medians.

Watch Out: What the Bootstrap Can't Do

The bootstrap treats the sample as if it were the population, so it inherits every flaw of the sample. If the sample is biased, so are the resamples. If the sample is small (say 10 values), it's a poor picture of the population, and bootstrap SEs can be badly off. And for statistics that depend on the extremes, such as the maximum, a resample can never exceed the largest value already seen, so the bootstrap fails. It's a powerful microscope, but it can only magnify what's already on the slide.

Check Yourself

Why does the bootstrap resample with replacement?

Think About It: A Statistic for a Jury

A forensic analyst measures the refractive index of 8 glass fragments from a suspect's clothing and reports their median. The defence asks: "How precise is that median?" Could the analyst answer with a bootstrap? What are the risks with only 8 fragments, and what else would you want to know?

A bootstrap is possible, and it would give a rough SE, but with 8 values it's fragile: the resampled medians can take only a handful of values, and the 8 fragments may not represent the population of fragments well (they might all come from one part of the window). You'd also want to know whether the fragments are independent (fragments from one impact may share properties), how precise the measuring instrument is (Chapter 5's measurement error), and how variable the refractive index is within a single window, which forensic databases record. With small samples, assumptions carry most of the weight, and they should be stated openly.

Key Takeaways

  • When no formula exists, simulation gives the sampling distribution if the population is known.
  • The bootstrap resamples the data with replacement to estimate a statistic's standard error from a single sample.
  • The bootstrap inherits the sample's flaws: bias, small size and dependence on extremes.

Prediction and inference in this chapter

7.7 Interludes

Interlude (Optional): Four Ways to Converge

"X‾n approaches μ" can mean several different things, and probabilists distinguish them carefully.

  • In probability (the weak Law of Large Numbers): for each ε>0, Pr(|X‾n−μ|≥ε)→0. For each large n, a big error is unlikely.
  • Almost surely (the strong Law of Large Numbers, proved by Kolmogorov in the 1930s): with probability 1, the sequence X‾1,X‾2,… converges to μ. For almost every infinite sequence of coin tosses, the running proportion of heads tends to exactly 12.
  • In mean square: E[(X‾n−μ)2]=σ2/n→0. This is what Chebyshev used.
  • In distribution (the Central Limit Theorem): the CDF of the standardized X‾n converges to Φ. This says nothing about individual sequences, only about the shape of the distribution.

The differences matter in advanced probability, where one kind of convergence can hold without another. The Berry–Esseen theorem even quantifies the CLT: the error in the Gaussian approximation is at most a constant times E|X−μ|3/(σ3n), which is why skewed populations (large third moments) need bigger samples.

Why it matters: "converges" is not one idea but several, and the right one depends on whether you care about one long run (almost sure), a typical sample (in probability) or the shape of the error (in distribution).

Interlude (Optional): Random Walks, from Paris to Wall Street

Add up independent ±1 steps and you get a random walk: the position of a drunkard on a line, or a gambler's fortune. By the CLT, after n steps the position is approximately N(0,n), so the walker typically wanders about n from the start. Speed up time and shrink the steps in the right way, and the walk becomes a continuous, jagged path called Brownian motion, after the botanist Robert Brown, who in 1827 saw tiny particles from pollen grains jiggling in water.

In 1900, five years before Einstein used the same mathematics to explain Brown's jiggling (and to give convincing evidence that atoms exist), a French doctoral student, Louis Bachelier, used it to model prices on the Paris stock exchange. His thesis was largely ignored for half a century. Today its descendants, stochastic calculus and the Black–Scholes formula for pricing options, run trillions of dollars of financial markets, and they also describe the diffusion of drugs through tissue and of pollutants through groundwater.

Why it matters: the n law is everywhere. It's why a diffusing smell takes four times as long to travel twice as far, and why a stock's typical yearly move is about 252 times its typical daily move.

Interlude (Optional): Why Machine Learning Works at All

A spam filter is trained on 100,000 emails and makes 3% errors on them. Will it make about 3% errors on next year's emails? Each email it classifies is like a draw from a population, and its error rate on new emails is a proportion, so the Law of Large Numbers suggests yes. But there's a catch: the filter was chosen because it did well on these particular emails, out of millions of possible filters. Choosing the best of many is selection, and selection makes the training error look better than the truth (the cover jinx of Section 6.4, again).

Statistical learning theory (the "probably approximately correct", or PAC, framework of Leslie Valiant, and the earlier work of Vladimir Vapnik and Alexey Chervonenkis) handles this with inequalities like Chebyshev's but much sharper, such as Hoeffding's inequality: for a proportion of n independent indicators, Pr(|p̂−p|≥ε)≤2e−2nε2. Combined over all the filters that might have been chosen, these bounds say how much data is needed before the training error can be trusted to predict the error on new data.

Why it matters: every claim that an AI system "achieves 97% accuracy" is a statement about a sampling distribution, and it's only as trustworthy as the independence of the test data from the data used to build the system.

7.8 Exercises

Each exercise has a hint and an answer. Before calculating, decide whether the question is about one observation, an average or a total, and whether the population is Gaussian or the CLT is needed.

Exercises

  1. Simulating Poisson means. (a) In R, draw a random sample of 40 values from Pois(10) and compute x‾. (b) Repeat 10,000 times and draw a histogram of the means. (c) What mean and SD should the sampling distribution have? Check against your simulation.

    means <- replicate(10000, mean(rpois(40, 10))). A Poisson has variance equal to its mean.

    (c) Mean 10 and SE 10/40=0.5. The histogram is close to Gaussian, centred at 10 with SD close to 0.5.

  2. Fingerprint matching. An automated fingerprint system takes a Gaussian time with mean 12 seconds and SD 1.8 seconds to process one fingerprint. (a) What is the probability that one match takes less than 10 seconds? (b) For 16 independent matches, what is the probability that the average is less than 11.5 seconds? (c) That the total is less than 184 seconds? (d) Two of these answers are equal for a reason, and one is equal by coincidence. Which is which?

    For (b), SE=1.8/4=0.45. For (c), T has mean 192 and SD 1.8×4=7.2.

    (a) Pr(Z<−1.11)≈0.133. (b) Pr(Z<−0.5/0.45)=Pr(Z<−1.11)≈0.133. (c) 184/16=11.5, so it's the same event as (b): 0.133. (b) and (c) are equal for a reason (the same event). (a) equals them by coincidence: 2/1.8 and 0.5/0.45 happen to be the same number.

  3. Malware scans. Scan times are right-skewed with mean 90 seconds and SD 30 seconds. (a) Why can't we assume one scan time is approximately Gaussian? (b) For 50 scans, approximate the probability that the average exceeds 95 seconds. (c) That the total exceeds 4,800 seconds?

    SE=30/50≈4.24. For (c), convert the total to an average.

    (a) The population is skewed; the CLT is about averages, not single values. (b) z=5/4.24=1.18, so the probability is about 0.119. (c) 4800/50=96, so z=6/4.24=1.41 and the probability is about 0.079.

  4. Tutoring centre. The improvement in test scores for students who visit a tutoring centre is right-skewed, with mean 10 points and SD 10 points. This year 40 students visited. (a) What is the probability that their average improvement exceeds 12 points? (b) That their total improvement is less than 390 points?

    SE=10/40≈1.58; 390 total is an average of 9.75.

    (a) z=2/1.58=1.26, and Pr(Z>1.26)≈0.103. (b) z=−0.25/1.58=−0.16, and Pr(Z<−0.16)≈0.437.

  5. Chebyshev versus the CLT. A fair coin is tossed 100 times. (a) Use Chebyshev to bound the probability that the proportion of heads differs from 0.5 by 0.1 or more. (b) Use the Gaussian approximation. (c) Compute the exact probability with pbinom and compare.

    SE=0.05, so 0.1 is 2 SEs. "Differs by 0.1 or more" means 40 or fewer heads, or 60 or more.

    (a) ≤122=0.25. (b) Pr(|Z|≥2)≈0.046 (0.057 with the continuity correction, using 59.5). (c) 2 * (1 - pbinom(59, 100, 0.5)) gives 0.057. Chebyshev is safe but very loose; the CLT with continuity correction is excellent.

  6. Planning a poll. How many randomly chosen voters are needed for a 95% margin of error of ±2 percentage points, in the worst case p=0.5?

    Solve 1.960.25/n=0.02 for n.

    n=(1.96/0.02)2×0.25=2,401.

  7. The two hospitals. In "before you start", a large hospital has 45 births a day and a small one 15. Assuming each baby is a boy with probability 0.5, independently, find the probability that more than 60% of a day's babies are boys at each hospital.

    More than 60% of 15 is 10 or more; of 45, 28 or more. Use 1 - pbinom(9, 15, 0.5) and the like.

    Small: Pr(X≥10)≈0.151. Large: Pr(X≥28)≈0.068. Over a year, about 55 such days at the small hospital and 25 at the large one.

  8. Pricing insurance. An insurer sells 10,000 policies. Each has a 1% chance of a claim of $20,000 in a year, independently. (a) Find the mean and SD of the total claims. (b) What premium per policy makes the probability that total claims exceed total premiums only 1%?

    Total claims =20,000×N with N∼Bin(10000,0.01). Use the CLT with z=2.326 for the 99th percentile.

    (a) Mean $2,000,000; SD 20,00099≈$199,000 (since 10,000×0.01×0.99=99). (b) Premiums must cover 2,000,000+2.326×199,000≈$2,463,000, so about $246 per policy: $200 for the expected claim, plus $46 of "safety loading" to cover the randomness.

  9. Sum of 100 dice. One hundred fair dice are rolled. Approximate the probability that their total is 380 or more.

    One die has mean 3.5 and variance 3512 (Example 6.10). Use a continuity correction at 379.5.

    Mean 350, SD 100×35/12≈17.08; z=29.5/17.08=1.73, and Pr(Z≥1.73)≈0.042.

  10. The gambler's fallacy. A friend says: "This roulette wheel has come up black 8 times in a row, so by the law of averages red is more likely now." Explain, using the Law of Large Numbers correctly, why the friend is wrong, and what the law does guarantee.

    Distinguish the proportion of reds from the number of reds minus blacks.

    Spins are independent, so the next one is red with probability 1838 regardless of the past. The Law of Large Numbers guarantees only that the proportion of reds approaches 1838 over a long run, because future spins swamp the past, not because they compensate for it.

  11. Bootstrapping Michelson. Using the 20 measurements from the first of Michelson's experiments, m <- morley$Speed[morley$Expt == 1] (Section 2.3), (a) find the median, and (b) estimate its standard error with 2,000 bootstrap resamples. (c) Why might the bootstrap be less trustworthy here than for the 141 rivers?

    sd(replicate(2000, median(sample(m, replace = TRUE)))).

    (a) 940 (km/s above 299,000). (b) About 28 km/s (it varies a little with the seed). (c) With only 20 values, the sample is a rough picture of the population, and the resampled medians can take only a few distinct values.

  12. Biased but big. A website asks visitors to vote on whether the city should build a new stadium, and 40,000 people respond: 70% yes. (a) Compute the SE of the proportion, as if the respondents were a random sample. (b) Explain why the true margin of error is probably much larger than this suggests.

    0.7×0.3/40000. Who chooses to vote on a stadium poll?

    (a) About 0.0023 (0.23 points), giving a 95% margin of about ±0.45 points. (b) Respondents chose themselves. People who care strongly (often fans, or committed opponents) respond far more than others, so the sample is biased, and a larger self-selected sample is no better. The SE formula covers only random sampling error.

  13. Correlated measurements. n measurements each have variance σ2, and every pair has correlation ρ>0. (a) Show that Var[X‾]=σ2n(1+(n−1)ρ). (b) What happens as n→∞? (c) Interpret for repeated measurements on a slightly miscalibrated instrument.

    Use Var[∑Xi]=∑Var[Xi]+2∑i<jCov(Xi,Xj) (Section 6.6), with (n2) pairs.

    (a) Var[∑Xi]=nσ2+2(n2)ρσ2=nσ2(1+(n−1)ρ); divide by n2. (b) It tends to ρσ2, not 0: the Law of Large Numbers fails. (c) A shared calibration error affects every measurement the same way, and no amount of averaging removes it; you need an independent check, such as a second instrument or a reference sample.

  14. A quick CLT check. Simulate 10,000 means of samples of size n=30 from an exponential population with mean 5. What fraction of them lie within 1.96 SEs of 5? Repeat with n=5. Comment.

    m <- replicate(10000, mean(rexp(30, 0.2))), then mean(abs(m - 5) <= 1.96 * 5 / sqrt(30)).

    Both fractions come out close to 0.95 (about 0.953 for n=30 and 0.957 for n=5 in one run), but the misses are split unevenly. For n=5, about 4.3% of the means fall above the upper limit and almost none below the lower one; for n=30 it's about 3.5% and 1.2%. A Gaussian would split them 2.5% and 2.5%. The sampling distribution is still right-skewed, and skewness is the last feature to disappear. Check the tails separately with mean(m > 5 + 1.96 * 5 / sqrt(30)).