Bootstrap lab

Three lines of R that replace pages of algebra

The whole bootstrap fits in one line of R:

replicate(B, statistic(sample(data, replace = TRUE)))

sample(data, replace = TRUE) draws a bootstrap sample of the same size as your data; replicate(B, ...) does it B times. Everything in this lab builds on that line.

Exercise 1: bootstrap the median

The 25 dice throws from the lesson are loaded as dice. Produce 2,000 bootstrap medians.

Inside median(), draw a resample of dice with replacement: sample(dice, replace = TRUE).

boot_medians <- replicate(2000, median(sample(dice, replace = TRUE)))
hist(boot_medians, breaks = seq(0.75, 6.25, by = 0.25), col = "#D6334A55",
     main = "2,000 bootstrap medians", xlab = "median")
boot_medians

Exercise 2: a percentile interval

The 95% percentile interval is simply the 2.5% and 97.5% quantiles of the bootstrap distribution.

c(0.025, 0.975)

ci_boot <- quantile(boot_means, probs = c(0.025, 0.975))
ci_t    <- t.test(dice)$conf.int
rbind(bootstrap = ci_boot, textbook_t = ci_t)

Exercise 3: bootstrap a regression slope

For regressions we resample whole rows (the “pairs bootstrap”), so each x stays with its own y. We use the heteroskedastic data from the regression lab, where the classical standard error was too small.

Same rule as always: TRUE.

boot_slopes <- replicate(2000, {
  rows <- sample(nrow(d), replace = TRUE)
  coef(lm(y ~ x, data = d[rows, ]))[["x"]]
})
c(bootstrap = sd(boot_slopes),
  classical = summary(lm(y ~ x, data = d))$coefficients["x", "Std. Error"])

Exercise 4: pick one bag

You must estimate the centre of a distribution and can use either the mean or the median. Which is better depends on the data. We’ll judge by mean squared error, \text{MSE} = \text{bias}^2 + \text{variance}, over many simulated samples where we know the true centre is 0.

The data here come from a t-distribution with 2 degrees of freedom: symmetric around 0, but with very heavy tails (occasional huge values). Compute each estimator’s MSE.

mean((medians - truth)^2)

truth <- 0
mse_mean   <- mean((means - truth)^2)
mse_median <- mean((medians - truth)^2)
c(mean = mse_mean, median = mse_median)

Now flip the situation. With light tails (close to normal data) the ranking reverses. Change df and rerun:

no estimator is best everywhere: it depends on the world you’re in

In real life you don’t know the truth, so you can’t run this simulation directly. But you can bootstrap from your data, treating your sample as the population, to compare the two estimators’ variability. That’s the “bootstrap in choosing estimators” idea from the slides.

Exercise 5: break the bootstrap

The lesson showed that the bootstrap fails for the sample maximum. Confirm the 63% figure: what share of bootstrap maxima are exactly equal to the original maximum?

The original estimate is max(x).

share_equal <- mean(boot_max == max(x))
c(share_equal = share_equal, theory = 1 - (1 - 1/200)^200)

And the fix, subsampling: draw m observations without replacement, with m much smaller than n.

Free play

Bootstrap anything you like. Some ideas: the 90th percentile of the incomes below; the correlation between two variables; the ratio of two means.

Back to the start: all lessons.