6. Sampling distributions¶

Draw a random sample and compute its mean. Draw another, and the mean is different. A third, different again. If a statistic changes every time we sample, how can it tell us anything reliable? The answer is that the changes themselves follow a law. This chapter studies that law, the sampling distribution: its centre (the true value, for a well-chosen statistic), its spread (the standard error, \(\sigma/\sqrt{n}\) for a mean, which we prove) and its shape (nearly normal, by the central limit theorem, often for quite small samples). You will see it for proportions, watch it emerge in simulations, and end with the bootstrap, Bradley Efron's 1979 idea for finding a sampling distribution when no formula is available: resample your own data.

In this chapter

You will learn:

what a sampling distribution is, and why it is the bridge between one sample and the population
the mean and standard error of a sample mean, with proofs, including the finite population correction for sampling without replacement
how to use the central limit theorem in practice, when "\(n \ge 30\)" is and is not enough, and why skewness fades like \(1/\sqrt{n}\)
the sampling distribution of a sample proportion, the continuity correction, and the quick margin of error \(1/\sqrt{n}\)
how to see all of this by simulation in a few lines of Python
the bootstrap: the plug-in idea, its exact answer for the mean, how to use it for a median, and where it fails

You need: Chapter 3 (mean, variance, standard deviation) and Chapter 5 (simple random samples). From Probability: Chapter 6 (expectation, variance of sums), Chapter 10 (the normal law) and Chapter 12 (the central limit theorem). Gentler version: the Topics lesson on sampling distributions.

Levels: 6.1 is beginner; 6.2 is intermediate with an advanced part on finite populations; 6.3 and 6.4 are intermediate; 6.5 is beginner to intermediate; 6.6 is intermediate with advanced proofs.

6.1 Statistics vary from sample to sample¶

A small population we can see whole¶

Usually the population is huge and unknown. To understand what sampling does, start with one so small we can list every possible sample. A made-up population of four numbers:

\[ \{2,\ 4,\ 6,\ 8\}, \qquad \mu = 5, \qquad \sigma^2 = \frac{(2-5)^2 + (4-5)^2 + (6-5)^2 + (8-5)^2}{4} = 5 . \]

(For a whole population we divide by \(N\); the divisor \(n-1\) belongs to samples, see Chapter 7.)

Take a simple random sample of size 2, without replacement. There are \(\binom{4}{2} = 6\) possible samples, each with probability \(1/6\):

Sample	2, 4	2, 6	2, 8	4, 6	4, 8	6, 8
Mean \(\bar{x}\)	3	4	5	5	6	7

So \(\bar{x}\) is a random variable: before we sample, its value is uncertain, with the distribution

\(\bar{x}\)	3	4	5	6	7
probability	1/6	1/6	2/6	1/6	1/6

This table is the sampling distribution of \(\bar{x}\): the distribution of the values a statistic takes over all possible samples, each weighted by its probability. Two facts stand out.

Its centre is the truth: \(E[\bar{x}] = (3 + 4 + 5 + 5 + 6 + 7)/6 = 5 = \mu\).
It is less spread than the population: its variance is \(\frac{(3-5)^2 + (4-5)^2 + 0 + 0 + (6-5)^2 + (7-5)^2}{6} = \frac{10}{6} \approx 1.67\), against \(\sigma^2 = 5\).

The words¶

A parameter describes the population (\(\mu\), \(\sigma\), \(p\)). It is a fixed number.
A statistic is computed from the sample (\(\bar{x}\), \(s\), \(\hat{p}\), the median). It is a random variable, because the sample is random.
The sampling distribution of a statistic is its probability distribution over repeated random sampling.
The standard error (SE) of a statistic is the standard deviation of its sampling distribution. It answers: "how far from its centre is this statistic typically?"

Do not confuse three distributions that all appear in this chapter:

the population distribution: the values of all \(N\) units;
the data distribution: the \(n\) values in your one sample (what a histogram of your data shows);
the sampling distribution: the values a statistic would take over all possible samples. You never see it directly; you see one value from it.

Statistical inference, the rest of this book, is the art of reasoning from one draw of (3) back to (1).

6.2 The sampling distribution of the mean¶

The standard model¶

From now on the usual model is: the data \(X_1, \dots, X_n\) are independent and identically distributed (iid), each with mean \(\mu\) and variance \(\sigma^2\). This describes sampling with replacement from a finite population, or sampling from a population so large that removing a few units changes nothing, or repeated measurements of something with random error.

Theorem (mean and standard error of \(\bar{X}\)). If \(X_1, \dots, X_n\) are iid with mean \(\mu\) and variance \(\sigma^2\), then the sample mean \(\bar{X} = \frac{1}{n}\sum_{i=1}^n X_i\) satisfies

\[ E[\bar{X}] = \mu, \qquad \operatorname{Var}(\bar{X}) = \frac{\sigma^2}{n}, \qquad \operatorname{SE}(\bar{X}) = \frac{\sigma}{\sqrt{n}} . \]

Proof. Expectation is linear (Probability, Chapter 6), so \(E[\bar{X}] = \frac{1}{n}\sum E[X_i] = \frac{1}{n} \cdot n\mu = \mu\).

For the variance, two rules: \(\operatorname{Var}(cY) = c^2\operatorname{Var}(Y)\), and the variance of a sum of independent variables is the sum of their variances (the covariances are zero). So

\[ \operatorname{Var}(\bar{X}) = \frac{1}{n^2}\operatorname{Var}\Big(\sum_{i=1}^n X_i\Big) = \frac{1}{n^2}\sum_{i=1}^n \operatorname{Var}(X_i) = \frac{1}{n^2} \cdot n\sigma^2 = \frac{\sigma^2}{n}. \]

Take the square root for the standard error. \(\blacksquare\)

Three lessons are packed into \(\sigma/\sqrt{n}\).

Averages are more stable than single values. That is why we average.
The square root law. To halve the standard error you need four times the data; to cut it by ten, a hundred times. Precision is expensive.
The population size does not appear. A random sample of 1,000 is about as precise for a country of 40 million as for a city of 400,000. What matters is \(n\), not \(n/N\) (as long as \(n\) is a small fraction of \(N\)).

Example 6.1 (a die). One roll of a fair die has \(\mu = 3.5\) and \(\sigma^2 = 35/12\), so \(\sigma \approx 1.708\). The mean of \(n\) rolls has standard error \(1.708/\sqrt{n}\):

\(n\)	1	4	25	100
SE of \(\bar{X}\)	1.708	0.854	0.342	0.171

For \(n = 2\) the sampling distribution is easy to list: the mean of two dice takes the values \(1, 1.5, \dots, 6\) with counts \(1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1\) out of 36, a triangle. Already one roll's flat distribution has become peaked in the middle.

The estimated standard error¶

The formula \(\sigma/\sqrt{n}\) needs \(\sigma\), which we usually do not know. We replace it with the sample standard deviation \(s\) and report the estimated standard error \(s/\sqrt{n}\). It is the number printed next to a mean in almost every scientific table. Replacing \(\sigma\) by \(s\) costs something when \(n\) is small; Chapter 8 shows exactly what (Student's \(t\)).

Example 6.2. A made-up sample of 25 measurements has \(\bar{x} = 41.2\) and \(s = 6.0\). The estimated standard error is \(6.0/\sqrt{25} = 1.2\). The data vary by about 6 around their mean; the mean would vary by only about 1.2 from one sample of 25 to another.

Sampling without replacement: the finite population correction¶

Level: advanced.

In the four-number example, the variance of \(\bar{x}\) was \(10/6\), not \(\sigma^2/n = 5/2\). The difference is that we sampled without replacement, so the draws were not independent: once 8 is drawn, it cannot be drawn again, and the next draw tends to be smaller.

Theorem (finite population correction). For a simple random sample of size \(n\) from a population of \(N\) values with mean \(\mu\) and variance \(\sigma^2\) (divisor \(N\)),

\[ E[\bar{X}] = \mu, \qquad \operatorname{Var}(\bar{X}) = \frac{\sigma^2}{n}\cdot\frac{N-n}{N-1}. \]

Proof. Let \(X_i\) be the value of the \(i\)-th unit drawn. By symmetry each \(X_i\) is a uniformly random member of the population, so \(E[X_i] = \mu\) and \(\operatorname{Var}(X_i) = \sigma^2\), which gives \(E[\bar{X}] = \mu\). Also by symmetry every pair \(X_i, X_j\) with \(i \ne j\) has the same covariance; call it \(c\).

To find \(c\), imagine drawing all \(N\) units (\(n = N\)). Then \(X_1 + \dots + X_N\) is the population total, a fixed number, so its variance is zero:

\[ 0 = \operatorname{Var}\Big(\sum_{i=1}^N X_i\Big) = N\sigma^2 + N(N-1)c, \qquad\text{so}\qquad c = -\frac{\sigma^2}{N-1}. \]

The covariance of two draws depends only on the population, not on how many draws come after, so the same \(c\) holds for any \(n\). Then

\[ \operatorname{Var}\Big(\sum_{i=1}^n X_i\Big) = n\sigma^2 + n(n-1)\Big(-\frac{\sigma^2}{N-1}\Big) = n\sigma^2\Big(1 - \frac{n-1}{N-1}\Big) = n\sigma^2\,\frac{N-n}{N-1}, \]

and dividing by \(n^2\) gives the result. \(\blacksquare\)

Check: for \(\{2,4,6,8\}\) with \(n = 2\): \(\frac{5}{2}\cdot\frac{4-2}{4-1} = \frac{5}{2}\cdot\frac{2}{3} = \frac{5}{3} \approx 1.67\), as we found by listing.

The factor \(\frac{N-n}{N-1}\), roughly \(1 - n/N\), is the finite population correction (fpc). It matters only when the sample is a sizeable fraction of the population. For a sample of 1,000 from 400,000 it is \(0.9975\); we ignore it. When \(n = N\) it is zero: a census has no sampling error.

6.3 The central limit theorem in practice¶

The statement we use¶

The theorem above gives the centre and spread of \(\bar{X}\). The central limit theorem (CLT) gives the shape. It is proved in Probability, Chapter 12; here is the version statisticians use.

Central limit theorem. If \(X_1, X_2, \dots\) are iid with mean \(\mu\) and finite variance \(\sigma^2 > 0\), then for every number \(z\),

\[ P\Big(\frac{\bar{X} - \mu}{\sigma/\sqrt{n}} \le z\Big) \to \Phi(z) \quad \text{as } n \to \infty, \]

where \(\Phi\) is the standard normal cumulative distribution function.

In words: for large \(n\), \(\bar{X}\) is approximately normal with mean \(\mu\) and standard deviation \(\sigma/\sqrt{n}\), whatever the shape of the population, as long as its variance is finite. We write \(\bar{X} \approx N(\mu, \sigma^2/n)\).

This is why the normal curve is everywhere in statistics. Individual measurements are often far from normal; averages and totals of many of them are close to normal.

Using it¶

Example 6.3 (made-up waiting times). Waiting times at a counter are exponential with mean 2 minutes (and so standard deviation 2). What is the chance that the average of 40 waits exceeds 2.5 minutes?

The SE is \(2/\sqrt{40} \approx 0.316\), so \(z = (2.5 - 2)/0.316 \approx 1.58\), and the CLT gives \(P(\bar{X} > 2.5) \approx 1 - \Phi(1.58) \approx 0.057\).

Here we can check exactly. A sum of \(n\) independent exponentials with mean 2 has a gamma distribution (Probability, Chapter 11), and the exact answer, computed with that gamma law, is about \(0.065\). The approximation is decent but not perfect, because the exponential is strongly skewed.

Example 6.4 (dice, with a continuity correction). Roll 25 dice. What is the chance their average is at least 4, that is, their sum is at least 100? The sum has mean \(87.5\) and standard deviation \(5 \times 1.708 = 8.54\).

Without correction: \(P(S \ge 100) \approx 1 - \Phi\big((100 - 87.5)/8.54\big) = 1 - \Phi(1.46) \approx 0.072\).
The sum takes only whole values, so "at least 100" is better matched by the normal area above 99.5 (the continuity correction): \(1 - \Phi\big((99.5 - 87.5)/8.54\big) = 1 - \Phi(1.41) \approx 0.080\).
Exact (adding up all \(6^{25}\) outcomes by computer convolution): \(0.0803\).

The corrected approximation is excellent. The die is symmetric, so the CLT works fast.

How large must \(n\) be?¶

Textbooks often say "\(n \ge 30\)". That is a rough rule, not a theorem, and it can be badly wrong in both directions.

For a symmetric population with light tails (a die, a uniform, a normal), \(n = 5\) or \(10\) is already fine.
For a skewed population, convergence is slower. The table compares the exact probability that the mean of \(n\) exponential waits lands above \(\mu + 1.645\,\sigma/\sqrt{n}\) (which the normal approximation says is 0.05) and below \(\mu - 1.645\,\sigma/\sqrt{n}\) (also 0.05 by the normal):

\(n\)	upper tail, exact	lower tail, exact
5	0.067	0.011
10	0.064	0.025
30	0.059	0.037
100	0.055	0.044

Even at \(n = 100\) the two tails are unbalanced. For populations with heavier skew (incomes, insurance claims, file sizes) you may need thousands.

For populations with infinite variance (the Cauchy law, or very heavy-tailed data), the CLT does not apply at all: see the failures in Probability, Chapter 12.

Why skewness fades like \(1/\sqrt{n}\)¶

Level: intermediate.

The skewness of a random variable \(Y\) with mean \(m\) and standard deviation \(\tau\) is \(\gamma(Y) = E[(Y - m)^3]/\tau^3\): zero for a symmetric law, positive when the right tail is longer. For the exponential law it is 2.

Proposition. If \(X_1, \dots, X_n\) are iid with skewness \(\gamma\), then \(\bar{X}\) has skewness \(\gamma/\sqrt{n}\).

Proof. Skewness does not change when we shift or rescale by a positive number, so we may study the sum \(S = \sum Y_i\) of the centred variables \(Y_i = X_i - \mu\). Expand \(E[S^3] = \sum_{i,j,k} E[Y_iY_jY_k]\). If one index differs from the other two, independence lets that factor split off with expectation \(E[Y_i] = 0\). So only the terms \(i = j = k\) survive, and \(E[S^3] = n\,E[Y_1^3] = n\gamma\sigma^3\). The variance of \(S\) is \(n\sigma^2\), so

\[ \gamma(S) = \frac{n\gamma\sigma^3}{(n\sigma^2)^{3/2}} = \frac{\gamma}{\sqrt{n}}. \qquad \blacksquare \]

For exponential waits, the mean of 40 still has skewness \(2/\sqrt{40} \approx 0.32\), which is why Example 6.3 was off. A useful habit: the more skewed your data look, the more you should distrust the normal approximation for small \(n\), and the more attractive the bootstrap of 6.6 becomes.

Other statistics¶

The CLT is about means, but many statistics behave like means for large \(n\) and are approximately normal too. One example, stated without proof: for a continuous population with density \(f\) and median \(m\) (with \(f(m) > 0\)), the sample median is approximately normal with mean \(m\) and standard error \(\frac{1}{2f(m)\sqrt{n}}\). For a normal population this is \(\sqrt{\pi/2}\,\sigma/\sqrt{n} \approx 1.25\,\sigma/\sqrt{n}\): the median is about 25 percent noisier than the mean for normal data, the price of its robustness (Chapter 3). Formulas like this need \(f(m)\), which is hard to estimate; the bootstrap avoids it.

6.4 Proportions¶

A proportion is a mean¶

Ask \(n\) randomly chosen people a yes/no question. Code each answer \(X_i = 1\) for yes and \(0\) for no. If a fraction \(p\) of the population would say yes, each \(X_i\) is a Bernoulli variable with \(E[X_i] = p\) and \(\operatorname{Var}(X_i) = p(1-p)\). The sample proportion is

\[ \hat{p} = \frac{\text{number of yes}}{n} = \bar{X}. \]

So everything in 6.2 and 6.3 applies at once. The count \(n\hat{p}\) has a binomial distribution (Probability, Chapter 5), and

\[ E[\hat{p}] = p, \qquad \operatorname{SE}(\hat{p}) = \sqrt{\frac{p(1-p)}{n}} . \]

For large \(n\), \(\hat{p} \approx N\big(p,\ p(1-p)/n\big)\). A common rule says the approximation is reasonable when both \(np\) and \(n(1-p)\) are at least 10 (some books say 5); like "\(n \ge 30\)", it is a guide, not a guarantee.

Example 6.5. In a population 30 percent own a bicycle (made up). In a random sample of 200, what is the chance that 25 percent or fewer own one?

SE: \(\sqrt{0.3 \times 0.7/200} \approx 0.0324\).
Normal approximation: \(\Phi\big((0.25 - 0.30)/0.0324\big) = \Phi(-1.54) \approx 0.061\).
With continuity correction (\(\hat{p} \le 0.25\) means count \(\le 50\), so use 50.5): \(\Phi\big((50.5/200 - 0.3)/0.0324\big) = \Phi(-1.47) \approx 0.071\).
Exact binomial sum: \(P(\text{count} \le 50) \approx 0.0695\).

The largest standard error, and the quick margin of error¶

The function \(p(1-p)\) is largest at \(p = \tfrac12\), where it equals \(\tfrac14\) (since \(p(1-p) = \tfrac14 - (p - \tfrac12)^2\)). So for every \(p\),

\[ \operatorname{SE}(\hat{p}) \le \frac{1}{2\sqrt{n}} . \]

About 95 percent of a normal distribution lies within 2 standard deviations of its centre, so \(\hat{p}\) is within \(2 \cdot \frac{1}{2\sqrt{n}} = \frac{1}{\sqrt{n}}\) of \(p\) about 95 percent of the time, whatever \(p\) is. This is the pollster's quick margin of error: for \(n = 1000\), \(1/\sqrt{1000} \approx 0.032\), the familiar "plus or minus 3 points". It covers only random sampling error, not any of the biases of Chapter 5. Chapter 8 makes it precise.

6.5 Simulation: seeing it happen¶

A sampling distribution is a statement about all the samples you could have drawn. A computer can draw thousands of them in a second, so we can look at the distribution directly instead of trusting a formula.

The first person to do this was William Sealy Gosset ("Student"), in his 1908 paper on the \(t\) distribution (Chapter 8). With no computer, he wrote the heights and finger measurements of 3,000 criminals (from a published anthropometric data set) on pieces of card, shuffled them, and drew 750 samples of 4, to check his formula for the distribution of a standardised mean. Today it takes three lines.

Means of exponential waits¶

import random, statistics as st
random.seed(1)
for n in (1, 5, 40):
    means = [st.fmean(random.expovariate(0.5) for _ in range(n))   # mean-2 exponential waits
             for _ in range(10000)]                                 # 10,000 samples of size n
    print(n, round(st.fmean(means), 3), round(st.pstdev(means), 3))

The output, with the theoretical standard error \(2/\sqrt{n}\) added:

\(n\)	average of the 10,000 means	SD of the 10,000 means	theory \(2/\sqrt{n}\)
1	1.999	2.009	2.000
5	2.018	0.900	0.894
40	1.998	0.315	0.316

The simulated centre matches \(\mu = 2\) and the simulated spread matches \(\sigma/\sqrt{n}\), up to the small error that comes from using 10,000 samples rather than infinitely many. A histogram of the means list (try statistics.quantiles or a plotting library) shows the shapes in the figure of 6.3: skewed for \(n = 5\), nearly normal for \(n = 40\).

Proportions¶

random.seed(7)
props = [sum(random.random() < 0.3 for _ in range(200)) / 200 for _ in range(10000)]
print(st.fmean(props), st.pstdev(props), sum(p <= 0.25 for p in props) / 10000)

This prints about \(0.300\), \(0.032\) and \(0.066\): the centre \(p = 0.3\), the SE \(0.0324\), and a tail probability close to the exact \(0.0695\) of Example 6.5.

How accurate is a simulation?¶

A simulation is itself a sample, of size \(B\) (here \(B = 10{,}000\)), from the sampling distribution. A probability estimated as a fraction of \(B\) simulated samples has standard error \(\sqrt{q(1-q)/B}\), where \(q\) is the true probability. For \(q \approx 0.07\) and \(B = 10{,}000\) that is about \(0.0026\), so \(0.066\) against \(0.0695\) is a perfectly ordinary difference. To get one more decimal place, multiply \(B\) by 100. This Monte Carlo error is separate from, and usually much smaller than, the sampling error we are studying.

What simulation is good for¶

Simulation answers "what if" questions that formulas cannot: the sampling distribution of a median, of a trimmed mean, of a ratio, of a correlation; the effect of skewness, outliers or a small \(n\); whether a method keeps its promised error rate. Statisticians use it constantly to check methods before trusting them. Its one limitation is that it needs the population to be known. With real data we do not know the population, which leads to the bootstrap.

6.6 The bootstrap idea¶

Level: intermediate, with advanced proofs.

The problem¶

You have one sample, \(x_1, \dots, x_n\). You compute a statistic, perhaps the median, and you want its standard error. For the mean there is a formula, \(s/\sqrt{n}\). For the median the formula in 6.3 needs the population density at the median, which you do not know. For a trimmed mean, a ratio, or the largest eigenvalue of a correlation matrix there may be no usable formula at all.

If you knew the population, you could simulate: draw many samples of size \(n\), compute the statistic for each, and look at the spread. You do not know the population. But you have the next best thing: the sample itself.

The plug-in principle¶

The empirical distribution \(\hat{F}\) of the data is the distribution that puts probability \(1/n\) on each observed value \(x_1, \dots, x_n\). It is our best guess of the population. The bootstrap says: treat \(\hat{F}\) as if it were the population, and simulate from it.

Drawing a sample of size \(n\) from \(\hat{F}\) means drawing \(n\) times with replacement from the data. Such a sample is a bootstrap resample, written \(x_1^*, \dots, x_n^*\). Some data values appear in it more than once, others not at all.

The bootstrap standard error.

Draw \(B\) resamples (typically \(B\) = 1,000 to 10,000), each of size \(n\), with replacement from the data.
For each, compute the statistic: \(\hat{\theta}^*_1, \dots, \hat{\theta}^*_B\).
The standard deviation of these \(B\) values is the bootstrap estimate of the standard error.

The name comes from the phrase "to pull oneself up by one's bootstraps": the data are used to assess themselves. Bradley Efron introduced the method in "Bootstrap methods: another look at the jackknife" in the Annals of Statistics in 1979. It was one of the first statistical methods designed around cheap computing, and it changed everyday practice.

Example¶

Example 6.6 (made-up data). Ten walking times, in minutes:

\[ 12,\ 15,\ 9,\ 22,\ 14,\ 31,\ 11,\ 17,\ 13,\ 19 . \]

The mean is \(16.3\), the median is \(14.5\), \(s \approx 6.447\), and the usual standard error of the mean is \(s/\sqrt{10} \approx 2.04\).

import random, statistics as st
x = [12, 15, 9, 22, 14, 31, 11, 17, 13, 19]
random.seed(2026)
boot_means, boot_medians = [], []
for _ in range(2000):
    r = random.choices(x, k=len(x))          # one resample, with replacement
    boot_means.append(st.fmean(r))
    boot_medians.append(st.median(r))
print(st.stdev(boot_means), st.stdev(boot_medians))

This prints about \(1.98\) for the mean and about \(2.07\) for the median. The bootstrap SE of the mean agrees well with the formula. For the median we had no formula, and now we have an answer.

What the bootstrap gives for the mean, exactly¶

Proposition. With infinitely many resamples (\(B \to \infty\)), the bootstrap variance of the sample mean equals \(\hat{\sigma}^2/n\), where \(\hat{\sigma}^2 = \frac{1}{n}\sum_{i=1}^n (x_i - \bar{x})^2\) is the variance of the data computed with divisor \(n\).

Proof. Under \(\hat{F}\), each resampled value \(x_j^*\) is one of \(x_1, \dots, x_n\) with probability \(1/n\) each, independently. So \(E^*[x_j^*] = \bar{x}\) and \(\operatorname{Var}^*(x_j^*) = \frac{1}{n}\sum_i (x_i - \bar{x})^2 = \hat{\sigma}^2\) (the star means "over resampling, with the data held fixed"). The resample mean \(\bar{x}^*\) is a mean of \(n\) iid draws, so by the theorem of 6.2 its variance is \(\hat{\sigma}^2/n\). \(\blacksquare\)

For Example 6.6, \(\hat{\sigma}/\sqrt{n} = \sqrt{9/10}\cdot s/\sqrt{10} \approx 1.934\). The simulated \(1.98\) differs from it only by Monte Carlo error. The ideal bootstrap is slightly smaller than \(s/\sqrt{n}\) because it uses divisor \(n\); for \(n\) of 20 or more the difference hardly matters.

Some facts about resamples¶

How many different resamples are there? A resample is determined by how many times each of the \(n\) data points is used: \(n\) non-negative counts adding to \(n\). By "stars and bars" (Probability, Chapter 2) there are \(\binom{2n-1}{n}\) such count patterns, which is \(92{,}378\) for \(n = 10\). They are not equally likely: the pattern using each point once has probability \(10!/10^{10} \approx 0.00036\), while a pattern using one point ten times has probability \(10^{-10}\). For \(n = 20\) there are about \(6.9 \times 10^{10}\) patterns, which is why we sample \(B\) of them at random rather than listing them all.

How much of the data is in a resample? A given data point is missed by one draw with probability \(1 - 1/n\), and by all \(n\) independent draws with probability \((1 - 1/n)^n\). This tends to \(e^{-1} \approx 0.368\), a standard limit from calculus. So a resample contains, on average, about 63.2 percent of the distinct data points (\(65.1\) percent for \(n = 10\)).

When the bootstrap fails¶

The bootstrap works when the statistic depends smoothly on the whole distribution, as means, medians (with care), correlations and regression coefficients do. It can fail badly when the statistic depends on the extreme of the data.

Example: the maximum. Suppose the data are uniform on \([0, \theta]\) and we estimate \(\theta\) by the sample maximum \(M\) (we meet this estimator in the German tank problem of Chapter 7). The true sampling distribution of \(M\) is continuous: \(P(M = \theta) = 0\). But a resample's maximum equals the data maximum whenever the largest data point is drawn at least once, which happens with probability \(1 - (1 - 1/n)^n \to 1 - e^{-1} \approx 0.632\). So the bootstrap distribution puts about 63 percent of its weight on a single value, and looks nothing like the truth.

Other failures: very small samples (the data are a poor picture of the population), heavy-tailed populations with infinite variance, and data that are not independent (resampling single time points destroys the dependence; block bootstraps exist for time series). Chapter 15 studies the bootstrap in depth, and Chapter 8 uses it to build confidence intervals.

Two sources of error, kept apart¶

A bootstrap answer has two errors.

Bootstrap (statistical) error: \(\hat{F}\) is not the population. This depends on \(n\), and no amount of computing removes it.
Monte Carlo error: we used \(B\) resamples, not all of them. This shrinks like \(1/\sqrt{B}\) and is cheap to reduce.

For a standard error, \(B = 1{,}000\) to \(2{,}000\) is usually plenty; for the tails of a distribution (confidence intervals), more is better.

Check yourself¶

1. A population has mean 50 and standard deviation 12. For random samples of size 36 (from a large population), give the mean and standard error of the sample mean.

\(E[\bar{X}] = 50\) and \(\operatorname{SE} = 12/\sqrt{36} = 2\).

2. How large a sample is needed to make the standard error of a mean one quarter of what it is with \(n = 25\)?

The standard error is proportional to \(1/\sqrt{n}\). To divide it by 4 we must multiply \(n\) by \(16\): \(n = 400\).

3. From the population \(\{1, 2, 3, 4, 5\}\), list the sampling distribution of the mean of a simple random sample of size 2, and check the finite population correction formula.

The 10 samples have means \(1.5, 2, 2.5, 3, 2.5, 3, 3.5, 3.5, 4, 4.5\) (for \(\{1,2\}, \{1,3\}, \{1,4\}, \{1,5\}, \{2,3\}, \{2,4\}, \{2,5\}, \{3,4\}, \{3,5\}, \{4,5\}\)). Their mean is 3, and their variance is \((2.25 + 1 + 0.25 + 0 + 0.25 + 0 + 0.25 + 0.25 + 1 + 2.25)/10 = 7.5/10 = 0.75\). The population has \(\mu = 3\) and \(\sigma^2 = (4+1+0+1+4)/5 = 2\). The formula gives \(\frac{2}{2}\cdot\frac{5-2}{5-1} = 0.75\). They agree.

4. A population has standard deviation 15. Using the CLT, find the approximate probability that the mean of a random sample of 36 is more than 5 away from the population mean.

\(\operatorname{SE} = 15/6 = 2.5\), so 5 is 2 standard errors. \(P(|Z| > 2) = 2(1 - \Phi(2)) \approx 0.0455\).

5. Why can't the CLT be used to say that the mean of 10 incomes from a very skewed income distribution is approximately normal? What would you do instead?

The CLT is a limit statement; how fast it takes hold depends on the shape. A strongly skewed population gives sample means with skewness \(\gamma/\sqrt{n}\), still large for \(n = 10\) if \(\gamma\) is large. Use a larger sample, a transformation (such as the logarithm), or the bootstrap, and check by simulation.

6. A poll of 400 people finds 52 percent in favour. Give the estimated standard error of the proportion and the quick margin of error.

\(\sqrt{0.52 \times 0.48/400} \approx 0.0250\). The quick margin is \(1/\sqrt{400} = 0.05\), that is, plus or minus 5 points. (Two estimated standard errors, \(0.0500\), agree.)

7. Prove that \(p(1-p) \le 1/4\) for every \(p\) between 0 and 1, and say when equality holds.

\(\frac14 - p(1-p) = \frac14 - p + p^2 = (p - \frac12)^2 \ge 0\), with equality exactly when \(p = \frac12\).

8. A simulation estimates a probability as 0.20 from \(B = 2{,}500\) simulated samples. What is the Monte Carlo standard error? How large must \(B\) be to make it about 0.002?

\(\sqrt{0.2 \times 0.8/2500} = 0.4/50 = 0.008\). To divide it by 4 we need \(16\) times as many: \(B = 40{,}000\).

9. The data are 3, 7, 8. List all the distinct bootstrap resamples (as count patterns), give the number, and compute the ideal bootstrap variance of the mean.

Count patterns \((a, b, c)\) with \(a + b + c = 3\): there are \(\binom{5}{3} = 10\): \((3,0,0), (0,3,0), (0,0,3), (2,1,0), (2,0,1), (1,2,0), (0,2,1), (1,0,2), (0,1,2), (1,1,1)\). The data have mean 6 and divisor-\(n\) variance \(\hat{\sigma}^2 = (9 + 1 + 4)/3 = 14/3\). The ideal bootstrap variance of the mean is \(\hat{\sigma}^2/n = 14/9 \approx 1.556\).

10. For \(n = 5\), what is the probability that a particular data point is missing from a bootstrap resample? What is the expected number of distinct points in a resample?

\((1 - 1/5)^5 = 0.8^5 = 0.32768\). Each point is present with probability \(0.67232\), so by linearity the expected number of distinct points is \(5 \times 0.67232 \approx 3.36\).

11. Show that the sum of the third central moments argument in 6.3 fails without independence: give two identically distributed variables with \(E[(Y_1 + Y_2)^3] \ne E[Y_1^3] + E[Y_2^3]\).

Take \(Y_1 = Y_2 = Y\), any centred variable with \(E[Y^3] \ne 0\) (for example \(Y = X - 1\) with \(X\) exponential with mean 1, so \(E[Y^3] = 2\)). Then \(E[(2Y)^3] = 8E[Y^3] = 16\), while \(E[Y_1^3] + E[Y_2^3] = 4\). The cross terms \(3E[Y_1^2Y_2] + 3E[Y_1Y_2^2]\) do not vanish.

12. Data are uniform on \([0, \theta]\), \(n = 20\). What is the probability that a bootstrap resample has the same maximum as the data? Why does this show the bootstrap fails here?

\(1 - (1 - 1/20)^{20} = 1 - 0.95^{20} \approx 1 - 0.358 = 0.642\). So about 64 percent of the bootstrap distribution of the maximum sits on one value, while the true sampling distribution of the maximum is continuous with no atoms. The bootstrap picture of the variability is badly wrong.

What you learned¶

A statistic is a random variable; its sampling distribution describes its values over all possible samples, and its standard deviation is the standard error.
For iid data, \(E[\bar{X}] = \mu\) and \(\operatorname{SE}(\bar{X}) = \sigma/\sqrt{n}\); without replacement multiply the variance by \((N-n)/(N-1)\). We estimate the SE by \(s/\sqrt{n}\).
The CLT makes \(\bar{X}\) approximately normal for large \(n\); how large depends on skewness, which fades like \(1/\sqrt{n}\). "\(n \ge 30\)" is a rough guide only.
A proportion is a mean of zeros and ones, with SE \(\sqrt{p(1-p)/n} \le 1/(2\sqrt{n})\); the quick 95 percent margin of error is \(1/\sqrt{n}\).
Simulation lets you see sampling distributions directly, with a Monte Carlo error of order \(1/\sqrt{B}\).
The bootstrap resamples the data with replacement to estimate a sampling distribution; for the mean it reproduces \(\hat{\sigma}/\sqrt{n}\), and it fails for statistics like the maximum.

Next: Chapter 7, what makes an estimator good, and two general ways to build one: the method of moments and maximum likelihood.