The Bootstrap

Every confidence interval so far needed the sampling distribution of an estimator, and each was available only because a theorem supplied it: the central limit theorem for large samples, an exact normal result for normal data, the $t$ distribution for a studentized mean. Outside those cases the sampling distribution is unknown, and for many estimators it is intractable.

The bootstrap gets the sampling distribution by simulation. The obstacle is that simulating from the true distribution requires knowing it, which is the original problem. The bootstrap’s answer is to substitute the best available estimate of that distribution, which is the data.

The bootstrap principle. Approximate the sampling distribution of a sample statistic by resampling many data sets, with replacement, from an estimated distribution of the data-generating distribution, and computing the statistic on each resampled data set.

1. Why resampling with replacement

The data are a random sample from the population. The empirical distribution, which places mass $1/n$ on each observed value, is an estimate of the population distribution, and by the law of large numbers the empirical CDF converges to the true CDF at every point.

So drawing a sample from the empirical distribution approximates drawing a sample from the population. Drawing $n$ values from a distribution putting mass $1/n$ on each of $x_1, \dots, x_n$ is exactly sampling with replacement from the data.

With replacement is not a technical detail, it is the whole mechanism. Sampling $n$ values without replacement from $n$ observations returns the original data set every time, in a different order, and every bootstrap statistic would be identical. Replacement is what makes the resamples differ: from $\{1, 3, 4\}$ a bootstrap sample might be $\{3, 1, 1\}$.

No new information enters. A bootstrap sample contains no data point that was not already observed. The bootstrap does not improve the estimate; it estimates the variability of the estimate.

2. The empirical bootstrap

Also called the non-parametric bootstrap, because it assumes nothing about the shape of the population distribution.

Inputs. A data set of size $n$; the estimate of interest $\hat\mu$ computed from it; the number of bootstrap samples $B$.

Steps. For each $b = 1, 2, \dots, B$:

1. Draw a bootstrap sample of size $n$ with replacement from the data set.
2. Compute the bootstrap statistic $\hat\mu_b^*$ from that sample.
3. Compute the centred bootstrap statistic $(\hat\mu - \mu)_b^* = \hat\mu_b^* - \hat\mu$.

Output. $B$ copies of $\hat\mu_b^* - \hat\mu$, which approximate the sampling distribution of the difference between the estimator and the parameter.

mu11_hat <- mean(ward_11_score)
B <- 1000
n <- length(ward_11_score)
mu11_hat_dist <- numeric(B)
for (b in seq(B)) {
  boot <- sample(ward_11_score, n, replace = TRUE)
  mu11_hat_dist[b] <- mean(boot) - mu11_hat
}

Why the statistic is centred

Step 3 subtracts $\hat\mu$, and the reason is the point of the method.

What you want is the distribution of $\hat\mu - \mu$: how far the estimator falls from the truth. You cannot compute it, because $\mu$ is unknown.

The bootstrap world is a copy of the real one, one level down. In the real world the population has mean $\mu$ and the sample gives $\hat\mu$. In the bootstrap world the “population” is the observed data, whose mean is exactly $\hat\mu$, and each resample gives $\hat\mu_b^*$. So $\hat\mu$ plays the role of the unknown parameter, and

$\hat\mu_b^* - \hat\mu \quad\text{approximates}\quad \hat\mu - \mu$

The quantity you need has the unknown in it; the quantity you can compute replaces the unknown by something known and shifts everything down one level.

3. Worked example: estimating a tail probability

Ward 11 building evaluation scores, with $\hat\mu_{11}$ the sample mean. Estimate $P\left(|\bar{X}_{11} - \mu_{11}| > 0.5\right)$: the chance the estimate is off by more than half a point.

The distribution of $\bar{X}_{11} - \mu_{11}$ is unknown. Replace it with the bootstrap distribution of $\hat\mu_b^* - \hat\mu_{11}$ and compute the proportion exceeding 0.5 in absolute value:

mean(abs(mu11_hat_dist) > .5)
# [1] 0.731

So roughly 73% of the time the sample mean misses by more than 0.5. No distributional assumption was made, and no formula for the sampling distribution was needed.

4. The bootstrap confidence interval

The same idea supplies critical values for a confidence interval when the studentized mean’s distribution is unknown, which is the small-sample non-normal case.

The pivot is the studentized mean $\frac{\bar{X}_n - \theta}{S_n/\sqrt{n}}$. Its bootstrap analogue is computed from each resample:

$t_b^* = \frac{\bar{x}_n^* - \bar{x}_n}{s_n^*/\sqrt{n}}$

with the bootstrap sample’s own mean and standard deviation in the numerator and denominator, and the original $\bar{x}_n$ playing the role of the parameter.

Take empirical quantiles of the $B$ values of $t_b^*$ to get critical values $c_l^*$ and $c_u^*$, then invert:

$\left(\bar{x}_n - c_u^*\frac{s_n}{\sqrt{n}},\ \ \bar{x}_n - c_l^*\frac{s_n}{\sqrt{n}}\right)$

B <- 1000
alpha <- 0.1
t_boot <- numeric(B)
for (b in seq(B)) {
  boot <- sample(software, n, replace = TRUE)
  t_boot[b] <- (mean(boot) - x_bar) / (sd(boot) / sqrt(n))
}
cl <- quantile(t_boot, alpha / 2)
cu <- quantile(t_boot, 1 - alpha / 2)
x_bar - c(cu, cl) * s_n / sqrt(n)

Two features of that formula are easy to get wrong. Both endpoints subtract, and the upper critical value appears in the lower endpoint. Both follow from the inversion in Confidence Intervals: multiplying through by $-1$ reverses the inequalities and swaps the roles of the two critical values. For a symmetric pivot the swap is invisible because $c_l^* = -c_u^*$; the bootstrap distribution need not be symmetric, so here it matters.

Worked example: software reliability

Times between software failures, heavily right-skewed, sample mean 649.14. Two 90% intervals for the mean:

They differ, and the reason is the assumption each makes. The $t$ interval assumes the data are a random sample from a normal distribution, which is false here: the histogram is strongly skewed. It produces an interval symmetric about 649.14 by construction.

The bootstrap makes no such assumption. It resamples the actual data, so the skewness in the data is carried into the bootstrap distribution, and the resulting interval is not centred on the sample mean. The sample mean sits at 649.14 while the interval runs from 523 to 832, further above than below. That asymmetry is a feature of the data that the $t$ method cannot represent.

5. The parametric bootstrap

The empirical bootstrap resamples the data. The parametric bootstrap simulates from a fitted model instead.

Use it when you are willing to assume a distributional family but the sampling distribution of your statistic is still intractable.

Inputs. A data set of size $n$; estimates of the model parameters, for instance $\hat\mu$ and $\hat\sigma$; the number of bootstrap samples $B$.

Steps. For each $b = 1, \dots, B$:

1. Draw a sample of size $n$ independently from the estimated distribution $F_{\hat\theta}$, for instance $\mathcal{N}(\hat\mu, \hat\sigma^2)$.
2. Compute the bootstrap statistic $\hat\mu_b^*$.
3. Compute the centred statistic $\hat\mu_b^* - \hat\mu$.

mu11_hat <- mean(ward_11_score)
sd_mu11_hat <- sd(ward_11_score)
B <- 1000
n <- length(ward_11_score)
mu11_hat_dist_norm <- numeric(B)
for (b in seq(B)) {
  boot <- rnorm(n, mean = mu11_hat, sd = sd_mu11_hat)
  mu11_hat_dist_norm[b] <- mean(boot) - mu11_hat
}

The difference in one line: the parametric bootstrap can generate values that were never observed, and the empirical bootstrap cannot. rnorm produces new numbers from a fitted curve; sample only reshuffles what is already there.

That is the trade. The parametric version smooths over gaps in a small sample and can extrapolate into the tails, and it is wrong if the assumed family is wrong. The empirical version assumes nothing and cannot produce any value outside the observed range, which makes it poor for tail quantities in small samples.

For the exponential, the fitted rate comes from the maximum likelihood estimator $\hat\lambda = 1/\bar{x}_n$, and the resampling step is rexp(n, rate = 1/mu_hat).

6. What the bootstrap does not do

It does not improve the estimate. Averaging bootstrap statistics does not give a better estimate than the original sample statistic.

The simulation makes the point. Draw a population of ten million from a geometric distribution, take a sample of 100, then bootstrap:

mu               # true mean       98.92987
mu_hat           # sample mean    112.18
mean(boot_dist)  # bootstrap mean 112.8448

The sample mean missed the truth by 13. The bootstrap mean is 112.84, which reproduces the sample’s error rather than correcting it, and could not do otherwise: the bootstrap resamples the sample, so it inherits everything the sample got wrong. As $B$ grows the bootstrap mean converges to $\hat\mu$, not to $\mu$.

Use the estimate from the whole data set. The bootstrap distribution tells you how variable that estimate is, and nothing about which direction it is wrong in.

$B$ and $n$ do different jobs. $B$ controls the simulation error and is free, so make it large. $n$ controls the statistical error and is fixed by the data. A bootstrap with $B = 10^6$ on 15 observations still only knows what those 15 observations know.

Vocabulary to deploy

• Bootstrap principle: approximate a sampling distribution by resampling from an estimate of the data-generating distribution.
• Empirical or non-parametric bootstrap: resample the data with replacement.
• Parametric bootstrap: simulate from a fitted parametric model.
• Bootstrap sample, bootstrap statistic $\hat\mu_b^*$, centred bootstrap statistic $\hat\mu_b^* - \hat\mu$.
• Bootstrap critical values $c_l^*$, $c_u^*$, as empirical quantiles of the bootstrap studentized statistic.
• $B$, the number of bootstrap samples, controlling simulation error rather than statistical error.