ECON 4370 / 6370 Computing for Economics

Lecture 10: Bootstrap and Resampling Methods

Zhan Gao

29 September 2026

One number is never enough

Estimate	Computed from	What it cannot tell us on its own
Average online spending of 10,000 households: $1,946 a year	one sample	how far the population average could be from it
Share of surveyed voters backing a candidate	one survey	whether another survey would agree
A regression coefficient	one dataset	whether the effect is distinguishable from zero

Every estimate comes from one sample. Today we ask how much it would move in another sample or to what extent we can quantify the uncertainty given the sample we have: first with theory (the law of large numbers and the central limit theorem), then by resampling the data we already have (the bootstrap).

Today’s route

1. Population and sample: two kinds of objects.
2. Law of large numbers and central limit theorem: what happens as n grows.
3. Standard errors and confidence intervals: quantifying uncertainty from one sample.
4. The bootstrap: resampling as an all-purpose substitute for formulas.
5. Online spending: the bootstrap on real data.

Other bootstrap intervals, a regression preview and two failure cases follow the main lesson.

Packages and data

Run the examples from the lectures/ directory.

library(tidyverse)
library(patchwork)   # side-by-side plots

Data	Unit	Used for
Simulated draws from an exponential distribution	fictitious online shoppers	sampling, LLN, CLT, bootstrap
data/web-browsers.csv	10,000 households, one year of online spending	the bootstrap on real data

Random numbers come from R’s pseudo-random generator, so set.seed() makes every chunk reproducible. The appendix also uses the boot package.

Population and sample

---

Two kinds of objects, and why the distinction matters

Population versus sample

• A population is the large pool of units we want to learn about: individuals, counties, firms, states.
  • Example: the proportion of US citizens who support a presidential candidate.
• Resource constraints mean the econometrician observes only a random sample of units and uses it to learn about the population.
  • Example: a survey of potential voters. Estimate the proportion above by the fraction of supporters in the survey.
• We model the population as a probability distribution and the sample as draws from that distribution.

Why bother with an abstract population?

The population lets us use the language of probability theory to

• compare procedures, for example LASSO against CART later in the course;
• assess uncertainty about estimates and predictions.

Sometimes the population is not a literal set of people. What is the population when we observe one time series of US GDP? Such conceptual questions matter, but we set them aside in this course.

Every object we meet in this course is either a population object or a sample object. It is critical to know which.

Sample objects estimate population objects

Sample object (computed from the data)	Population object (never observed)
average \bar{X}	expectation \mathbf{E}[X]
proportion	probability
histogram	density or mass function
model estimate \hat{\beta}	model parameter \beta

Sample objects are things we see directly in the data. Population objects are things we never see but want to learn about.

The rest of today asks how well a sample object approximates its population counterpart, and how to measure the gap.

A fictitious population: online spending

Suppose the amount an Amazon customer spends follows an exponential distribution with mean 1, in tens of dollars. dexp() is its density.

grid <- seq(0, 10, length.out = 1000)
dens <- data.frame(x = grid, y = dexp(grid))
p_dens <- ggplot(dens, aes(x, y)) +
  geom_line(colour = smu_red, linewidth = 1) +
  labs(x = "Spending (tens of dollars)",
       y = "Density")

Most customers spend little and a few spend a lot.

Areas under the density are probabilities

The shaded area is the population proportion of customers who spend at most $10, or equivalently the probability that a random customer spends at most $10.

pexp() is the cumulative distribution function: pexp(x) is the probability of spending at most x (times ten dollars).

pexp(1) - pexp(0)

## [1] 0.6321206

Is this a population or a sample quantity?

Answer: a population quantity. It is a feature of the distribution, computed without any data. Amazon might know the distribution from its transaction records, but the econometrician usually cannot.

Four functions for every distribution in R

Prefix	Meaning	Exponential	Normal
d	density (or mass) function	dexp()	dnorm()
p	cumulative probability \mathbf{P}(X \le x)	pexp()	pnorm()
q	quantile, the inverse of p	qexp()	qnorm()
r	random draws	rexp()	rnorm()

The same pattern covers unif, binom, pois, t, chisq and many more; see ?Distributions.

The d, p and q functions describe population objects. Only r produces a sample.

Drawing a random sample

set.seed(0)      # reproducible pseudo-random numbers
n <- 20
spend <- rexp(n)
round(spend, 2)

##  [1] 0.18 0.15 0.14 0.44 2.89 1.23 0.54 0.96 0.15 1.39 0.76 1.24 4.42 1.05 1.04 1.88 0.65 0.34 0.59
## [20] 2.36

This draws 20 customers from the fictitious population.

The econometrician sees only these 20 numbers, never the curve behind them.

rexp(n) uses rate 1 by default, so the mean is 1. Amazon might have a very good idea of the distribution from the transaction data it stores; we usually only get a random sample.

The sample distribution of spending

p_hist20 <- ggplot(data.frame(spend), aes(spend)) +
  geom_histogram(aes(y = after_stat(density)),
                 binwidth = 0.5, boundary = 0,
                 fill = smu_blue, colour = "white") +
  coord_cartesian(xlim = c(0, 10), ylim = c(0, 1)) +
  labs(x = "Spending (tens of dollars)",
       y = "Density")

mean(spend <= 1)

## [1] 0.55

mean(spend)

## [1] 1.119912

The bars between 0 and 1 give the sample fraction spending at most $10: 0.55 against the population value 0.632. The sample mean 1.12 sits closer to \mathbf{E}[X] = 1.

Overlay the population density

p_compare20 <- p_hist20 +
  geom_line(data = dens, aes(x, y),
            colour = smu_red, linewidth = 1)

The histogram is a sample object. The red density is the population object it estimates.

At n = 20 the two do not resemble each other much.

Let us repeat the experiment with larger samples.

Repeat with n = 20, 100, 500 and 2,500

set.seed(1)
sizes <- c(20, 100, 500, 2500)
samples <- lapply(sizes, rexp)            # a list of four samples
sample_plot <- function(s) {
  ggplot(data.frame(spend = s), aes(spend)) +
    geom_histogram(aes(y = after_stat(density)), binwidth = 0.5, boundary = 0,
                   fill = smu_blue, colour = "white") +
    geom_line(data = dens, aes(x, y), colour = smu_red, linewidth = 1) +
    coord_cartesian(xlim = c(0, 10), ylim = c(0, 1)) +
    labs(title = paste("n =", length(s)), x = "Spending (tens of dollars)", y = "Density")
}
t(sapply(samples, function(s) c(n = length(s), fraction = mean(s <= 1), mean = mean(s))))

##         n fraction      mean
## [1,]   20    0.550 1.0893254
## [2,]  100    0.600 0.9902038
## [3,]  500    0.624 0.9818994
## [4,] 2500    0.618 1.0400167

Larger samples look more like the population

The histogram tracks the density more closely as n grows, and the fraction and mean settle near 0.632 and 1: the law of large numbers at work.

Checkpoint: sample or population?

1. mean(spend <= 1) and pexp(1): which is a sample object and which a population object?
2. Amazon records every transaction for a year. Is the average of those transactions a sample or a population object?
3. With one million draws, would the histogram ever coincide with the density exactly?

Answers: (1) the first is a sample fraction, the second a population probability. (2) It depends on the question. If the population is this year’s actual customers, it is the population average. If this year’s transactions are one realisation of a process that could have turned out differently, it is a sample object and next year would give another draw. (3) No. A histogram of a finite sample is a sample object; each bar converges to the area under the density by the law of large numbers, but it never equals it exactly.

Law of large numbers and central limit theorem

---

What repeated sampling does to the sample mean

Notation

A random sample is n independent draws from the population distribution: \{X_i\}_{i=1}^n = \{X_1, \dots, X_n\}, \qquad \bar{X} = \frac{1}{n} \sum_{i=1}^n X_i .

X denotes an arbitrary draw. Its expectation \mathbf{E}[X] is the average of X in the population, and \text{Var}(X) = \mathbf{E}[(X - \mathbf{E}[X])^2] measures its spread. For the exponential distribution above, \mathbf{E}[X] = 1 and \text{Var}(X) = 1.

Two rules we will use:

Rule	Statement	Needs
Linearity	\mathbf{E}[aX + bY] = a\,\mathbf{E}[X] + b\,\mathbf{E}[Y]	nothing
Variance of a sum	\text{Var}(\sum_{i=1}^n X_i) = \sum_{i=1}^n \text{Var}(X_i)	independence

The law of large numbers

Law of large numbers (LLN). As n gets large, \bar{X} \approx \mathbf{E}[X].

More precisely: for any tolerance \varepsilon > 0, \mathbf{P}\left(\lvert \bar{X} - \mathbf{E}[X] \rvert > \varepsilon\right) \to 0 as n \to \infty.

In words, the average of a large random sample is close to the population expectation. This is what the simulation showed.

We can prove it in three steps: compute \mathbf{E}[\bar{X}], compute \text{Var}(\bar{X}), and combine them.

Step 1: the expectation of the sample mean

\begin{aligned} \mathbf{E}[\bar{X}] &= \mathbf{E}\left[\frac{1}{n} \sum_{i=1}^n X_i\right] = \frac{1}{n} \sum_{i=1}^n \mathbf{E}[X_i] \\ &= \frac{1}{n} \sum_{i=1}^n \mathbf{E}[X] = \mathbf{E}[X]. \end{aligned}

The second equality is linearity: the expectation of a sum is the sum of the expectations. The third uses that every X_i has the same distribution as X.

What does \mathbf{E}[\bar{X}] mean?

• Draw a random sample of size n once: one value of \bar{X}.
• Repeat B times: B separate draws of \bar{X}.
• Average over many, many repetitions: the result is \mathbf{E}[\bar{X}].

Step 2: the variance of the sample mean

A single \bar{X} need not be close to \mathbf{E}[X] when n is small; it fluctuates around it. How much is measured by its variance:

\begin{aligned} \text{Var}(\bar{X}) &= \text{Var}\left(\frac{1}{n} \sum_{i=1}^n X_i\right) = \frac{1}{n^2}\, \text{Var}\left(\sum_{i=1}^n X_i\right) \\ &= \frac{1}{n^2} \sum_{i=1}^n \text{Var}(X_i) = \frac{1}{n}\, \text{Var}(X). \end{aligned}

The second equality pulls the constant 1/n out as its square. The third uses independence: the variance of a sum of independent draws is the sum of the variances.

As n grows, \text{Var}(\bar{X}) shrinks towards zero. What does this say about the value of \bar{X}?

Step 3: combine the two

Chebyshev’s inequality bounds the probability of a large deviation by the variance:

\mathbf{P}\left(\lvert \bar{X} - \mathbf{E}[X] \rvert > \varepsilon\right) \le \frac{\text{Var}(\bar{X})}{\varepsilon^2} = \frac{\text{Var}(X)}{n\, \varepsilon^2} \longrightarrow 0 .

To recap: the typical value of \bar{X} is \mathbf{E}[X], because \mathbf{E}[\bar{X}] = \mathbf{E}[X], and its spread shrinks with n, because \text{Var}(\bar{X}) = n^{-1}\,\text{Var}(X). Together they give \bar{X} \approx \mathbf{E}[X].

Why we care. To learn the expectation of a population distribution, such as the average amount a customer spends, a large random sample does very well. At n = 2{,}500 the standard deviation of \bar{X} is 1/\sqrt{2500} = 0.02.

Chebyshev: for any random variable Z with finite variance, P(|Z - E[Z]| > eps) <= Var(Z)/eps^2. The LLN holds under the weaker condition that E|X| is finite, but the variance route is the elementary proof.

The sample mean has its own distribution

By now we understand the population distribution of X: monitor a random Amazon customer who just filled a basket and record the total. Do this for the universe of customers and the histogram of totals is the distribution of X.

Like X, \bar{X} is a random variable, so it has a population distribution too. How do we think of it?

• Draw a random sample of size n once: one \bar{X}. This is a single draw from the distribution of \bar{X}.
• Repeat B times: B draws from that distribution.
• Pick B very large and plot the histogram of all the draws: this is the population distribution of \bar{X}, also called its sampling distribution.

The central limit theorem

Central limit theorem (CLT). When n is large, the distribution of \bar{X} is close to a normal distribution with the expectation and variance we just derived: \bar{X} \approx \mathcal{N}\!\left(\mathbf{E}[X],\; n^{-1}\,\text{Var}(X)\right).

Formally, \sqrt{n}\,(\bar{X} - \mathbf{E}[X]) converges in distribution to \mathcal{N}(0, \text{Var}(X)), whatever the shape of the distribution of X, provided \text{Var}(X) is finite.

The exponential distribution is far from normal, so this is a strong claim. Let us check it by simulation.

A Monte Carlo experiment

Follow the recipe: draw a sample of size n and compute \bar{X}; do this B times; compare the histogram of \bar{X}_1, \dots, \bar{X}_B with the normal density \mathcal{N}(1, 1/n).

set.seed(2)
n <- 100       # sample size
B <- 1000      # number of repetitions
means <- rep(0, B)               # placeholder to fill
for (b in 1:B) {
  s <- rexp(n)
  means[b] <- mean(s)
}
c(mean = mean(means), variance = var(means), `1/n` = 1 / n)

##        mean    variance         1/n 
## 0.998033537 0.009626014 0.010000000

The average of the sample means is close to \mathbf{E}[X] = 1 and their variance is close to \text{Var}(X)/n = 0.01.

Plot the sample means against the normal density

mean_plot <- function(m, n) {
  grid <- seq(0.6, 1.4, length.out = 400)
  normal <- data.frame(x = grid, y = dnorm(grid, mean = 1, sd = 1 / sqrt(n)))
  ggplot(data.frame(m), aes(m)) +
    geom_histogram(aes(y = after_stat(density)), bins = 30,
                   fill = smu_blue, colour = "white") +
    geom_line(data = normal, aes(x, y), colour = smu_red, linewidth = 1) +
    coord_cartesian(xlim = c(0.6, 1.4)) +
    labs(title = paste0("n = ", n, ", B = ", length(m)),
         x = "Sample mean", y = "Density")
}

The red curve is the CLT’s prediction \mathcal{N}(1, 1/n), not a curve fitted to the histogram.

The histogram is bell-shaped at n = 100

A thousand sample means, each from 100 draws of a skewed population, follow the normal density closely.

At n = 500 the distribution concentrates

set.seed(3)
means500 <- replicate(B, mean(rexp(500)))    # the same loop in one line
c(mean = mean(means500), variance = var(means500), `1/n` = 1 / 500)

##        mean    variance         1/n 
## 0.999790718 0.001994272 0.002000000

Same axes: the spread shrinks by \sqrt{5}, from a standard deviation of 0.10 to 0.045.

Checkpoint: LLN versus CLT

1. Which theorem says the histogram of \bar{X} narrows as n grows, and which says it becomes bell-shaped?
2. Step 2 used independence. Where would the argument break for a time series of quarterly GDP?
3. The variance of the sample means was about 0.0096 at n = 100. What would you expect at n = 400?

Answers: (1) the LLN gives the narrowing, because the variance shrinks like 1/n; the CLT gives the shape. (2) Consecutive quarters are correlated, so the variance of the sum picks up covariance terms; the sample mean still converges under weak dependence, but its variance is no longer Var(X)/n. (3) About a quarter of it, 0.0025, since Var(X)/n = 1/400.

Standard errors and confidence intervals

---

Turning the CLT into a statement about one sample

From the CLT to a standard error

In practice we see one sample and want to know how far \bar{X} could be from \mathbf{E}[X].

The CLT says the spread of \bar{X} is \text{SD}(\bar{X}) = \sqrt{\text{Var}(X)/n}. We do not know \text{Var}(X), so we replace it by its sample analogue:

\text{SE}(\bar{X}) = \frac{\widehat{\text{SD}}(X)}{\sqrt{n}}, \qquad \widehat{\text{SD}}(X) = \sqrt{\frac{1}{n-1} \sum_{i=1}^n (X_i - \bar{X})^2},

where \widehat{\text{SD}}(X) is the sample standard deviation, sd() in R.

A standard error is the estimated standard deviation of an estimator across repeated samples.

A 95% confidence interval

\bar{X} \pm 2 \cdot \text{SE}(\bar{X})

Why …	Because …
… normal?	\bar{X} is approximately normal by the CLT.
… 2?	The middle 95% of a normal distribution lies within 1.96 standard deviations of its mean, and 1.96 \approx 2.
… this SE?	\text{SE}(\bar{X}) approximates \text{SD}(\bar{X}) = \text{Var}(\bar{X})^{1/2}.

qnorm(0.975)     # the 97.5% quantile of the standard normal

## [1] 1.959964

Two layers of approximation

1. The normal distribution stands in for the exact distribution of \bar{X} (the CLT).
2. The estimated \text{SE}(\bar{X}) stands in for the true \text{SD}(\bar{X}).

If we knew the distribution of \bar{X}, neither would be needed. Compute the interval for one sample of 500 customers:

set.seed(5)
n <- 500
spend <- rexp(n)
se <- sd(spend) / sqrt(n)
c(estimate = mean(spend), se = se)

##   estimate         se 
## 1.05845529 0.04701586

mean(spend) + c(-2, 2) * se

## [1] 0.9644236 1.1524870

The interval covers the population mean of 1, which in a real application we would not know.

What does “95%” promise?

The interval is random and \mathbf{E}[X] = 1 is fixed. Over repeated samples, 95% of the intervals should cover 1. In the R world we can check:

covers <- function(n, R = 2000) {
  set.seed(4)
  mean(replicate(R, {
    s <- rexp(n)
    abs(mean(s) - 1) <= 2 * sd(s) / sqrt(n)
  }))
}
ns <- c(10, 20, 100, 500)
sapply(setNames(ns, ns), covers)

##     10     20    100    500 
## 0.8715 0.9185 0.9555 0.9625

With a skewed population and n = 10, the normal approximation is poor and the interval covers too rarely. By n = 100 it keeps its promise.

Using 2 instead of 1.96 makes the nominal level 95.4%, so coverage slightly above 0.95 at n = 500 is expected. Each coverage rate is itself a sample proportion from 2,000 repetitions, with a standard error of about 0.005.

One sample, one sample mean

In the fictitious world of our R code we can draw \bar{X} B times for a very large B. The histogram of those draws is the distribution of \bar{X}, and we can read standard errors and quantiles straight off it.

In reality we never get to do this. We draw one sample \{X_i\}_{i=1}^n and compute one \bar{X}.

We cannot draw more samples from the population. Can we draw from something that resembles it?

The bootstrap

---

Resample the sample to mimic resampling the population

The bootstrap idea

• The LLN says that when n is large, the sample distribution of X (the histogram) is close to the population distribution (the density).
• If we want to draw from the population, we can do almost as well by drawing from the sample instead.
• The bootstrap recomputes the estimator on many such draws and uses the spread of the results as the sampling distribution of the estimator.

The bootstrap approximates the population distribution of an estimator (here \bar{X}) by the distribution of the same estimator across resamples of the data.

It is all-purpose: the same recipe works for means, medians, ratios, regression coefficients and much more.

Why yet another way of getting CIs?

• For simple estimators like the mean, the standard-error formula is easy. For complicated estimators, the bootstrap is often the easiest way to get a standard error.
• Bootstrap CIs preview related methods later in the course: random forests use bagging, short for bootstrap aggregation.
• Conceptually, the bootstrap sharpens the distinction between sample and population distributions.

The ideal: repeated draws from the population

• Top: the population density.
• Middle: five samples drawn from it, each producing one estimate \hat{\beta}_b. Here \hat{\beta} plays the role of \bar{X}.
• Bottom: the histogram of the estimates is the sampling distribution.

This is exactly what the Monte Carlo loop did. It requires the population, which we never have.

[Sources] - Taddy (2019), Business Data Science, chapter 1, figures on the bootstrap.

The bootstrap: repeated draws from the sample

• Top: the population density is replaced by the histogram of the one sample we have.
• Middle: five bootstrap samples drawn with replacement from it, each producing \hat{\beta}_b.
• Bottom: the histogram of the bootstrap estimates approximates the sampling distribution.

Everything below the top panel is the same recipe as before. Only the source of the draws has changed.

The bootstrap algorithm

1. Start with a random sample \{X_i\}_{i=1}^n.
2. For each b = 1, \dots, B:
   1. Resample with replacement n observations from \{X_i\}_{i=1}^n. Call this bootstrap sample \{X_i^b\}_{i=1}^n.
   2. Compute the estimator on the bootstrap sample. For the mean, call it \bar{X}_b.
3. \{\bar{X}_b\}_{b=1}^B is the bootstrap approximation to the population distribution of \bar{X}.

“With replacement” matters: a bootstrap sample repeats some observations and omits others. Without replacement, every resample would be the original sample in a different order.

Bootstrap standard errors and confidence intervals

The bootstrap standard error \text{SE}_{\text{boot}} is the standard deviation of \{\bar{X}_b\}_{b=1}^B, an alternative to \widehat{\text{SD}}(X)/\sqrt{n}.

Two 95% confidence intervals:

Interval	Formula	Uses
Normal	\bar{X} \pm 2 \cdot \text{SE}_{\text{boot}}	the bootstrap spread and the normal shape
Quantile	\bar{X} \pm q_{0.95}, where q_{0.95} is the 0.95 quantile of \lvert \bar{X}_b - \bar{X} \rvert over b	the bootstrap distribution directly

The quantile interval is preferred, because a few extreme bootstrap draws can inflate \text{SE}_{\text{boot}}. Standard errors are nonetheless reported more often because they are simpler.

One bootstrap sample in R

We have the sample of n = 500 customers in spend. Put X_1, \dots, X_n in a bag and draw from it n times with replacement:

set.seed(6)
idx <- sample.int(n, replace = TRUE)   # n indices drawn from 1:n with replacement
sort(idx)[1:10]                        # some indices repeat, others never appear

##  [1]  3  4  5  6  6  6  7  8 10 11

boot_sample <- spend[idx]
mean(boot_sample)                      # one bootstrap draw of the mean

## [1] 1.037495

length(unique(idx)) / n

## [1] 0.62

Each observation is left out of a bootstrap sample with probability (1 - 1/n)^n \approx e^{-1} = 0.37, so about 63% of the original observations appear at least once.

Repeat B = 1,000 times

set.seed(7)
boot_means <- rep(0, B)
for (b in 1:B) {
  idx <- sample.int(n, replace = TRUE)
  boot_means[b] <- mean(spend[idx])
}
c(sample_mean = mean(spend), boot_se = sd(boot_means))

## sample_mean     boot_se 
##   1.0584553   0.0468566

This is the CLT loop with one change: spend[idx] replaces rexp(n). We draw from the sample instead of the population.

In the R world we can also draw the truth for comparison:

set.seed(8)
true_means <- replicate(B, mean(rexp(n)))

Bootstrap versus truth: build the picture

grid <- seq(0.85, 1.25, length.out = 400)
draws <- bind_rows(
  data.frame(source = "Bootstrap: resample the sample", m = boot_means),
  data.frame(source = "Truth: redraw from the population", m = true_means))
normals <- bind_rows(
  data.frame(source = "Bootstrap: resample the sample", x = grid,
             y = dnorm(grid, mean(spend), sd(spend) / sqrt(n))),
  data.frame(source = "Truth: redraw from the population", x = grid,
             y = dnorm(grid, 1, 1 / sqrt(n))))
p_compare <- ggplot(draws, aes(m)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = smu_blue, colour = "white") +
  geom_line(data = normals, aes(x, y), colour = smu_red, linewidth = 1) +
  geom_vline(xintercept = 1, linetype = "dashed") +
  facet_wrap(~ source) +
  labs(x = "Sample mean", y = "Density")

The red curves are the normal approximations: \mathcal{N}(\bar{X}, \text{SE}^2) for the bootstrap and the CLT’s \mathcal{N}(1, 1/n) for the truth. The dashed line marks \mathbf{E}[X] = 1.

Bootstrap versus truth: the picture

Nearly the same shape and spread. The bootstrap distribution is centred at the sample mean \bar{X} = 1.06 rather than at 1, because it is built from the sample.

Read off the standard errors and intervals

c(bootstrap = sd(boot_means), formula = sd(spend) / sqrt(n),
  monte_carlo = sd(true_means), exact = 1 / sqrt(n))

##   bootstrap     formula monte_carlo       exact 
##  0.04685660  0.04701586  0.04604075  0.04472136

q95 <- unname(quantile(abs(boot_means - mean(spend)), 0.95))
ci <- rbind(boot_normal = mean(spend) + c(-2, 2) * sd(boot_means),
            boot_quantile = mean(spend) + c(-1, 1) * q95,
            formula = mean(spend) + c(-2, 2) * sd(spend) / sqrt(n))
colnames(ci) <- c("lower", "upper")
ci

##                   lower    upper
## boot_normal   0.9647421 1.152168
## boot_quantile 0.9689625 1.147948
## formula       0.9644236 1.152487

Four routes to the same spread, and three intervals that agree to within 0.005. The truth is close to normal by the CLT, so the bootstrap distribution is too.

Checkpoint: bootstrap mechanics

1. What would happen if step 2a resampled without replacement?
2. Raising B from 1,000 to 100,000: does \text{SE}_{\text{boot}} get closer to the true \text{SD}(\bar{X})?
3. The bootstrap intervals are centred at \bar{X} = 1.06, not at 1. Is that a problem?

Answers: (1) Every bootstrap sample would be a permutation of the original, every bootstrap mean would equal the sample mean, and SE_boot would be zero. (2) A larger B reduces the Monte Carlo error in computing the bootstrap distribution, but the gap between the sample and the population is fixed by n; only more data shrinks it. (3) No. In practice we never know E[X]; the interval is designed to cover it in 95% of samples, and this one does.

Online spending

---

The bootstrap on 10,000 households

The data

browser <- read.csv("./data/web-browsers.csv")
dim(browser)

## [1] 10000     7

head(browser, 3)

##   id anychildren broadband hispanic  race region spend
## 1  1           0         1        0 white     MW   424
## 2  2           1         1        0 white     MW  2335
## 3  3           1         1        0 white     MW   279

summary(browser$spend)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       1     162     510    1946    1523  401338

One row per household: annual online spending in dollars plus a few demographics. The mean is almost four times the median, the sign of a heavy right tail.

[Sources] - Taddy (2019), Business Data Science, chapter 1. The same file returns in lecture 10.

Spending is heavily skewed

p_raw <- ggplot(browser, aes(spend)) +
  geom_histogram(bins = 30, fill = smu_blue,
                 colour = "white") +
  labs(x = "Spending (dollars)", y = "Households")
p_log <- ggplot(browser, aes(log(spend))) +
  geom_histogram(bins = 30, fill = smu_gold,
                 colour = "white") +
  labs(x = "log(spending)", y = "Households")

In dollars, a handful of households with enormous spending stretch the axis. In logs the distribution is roughly bell-shaped.

The CLT still applies to the mean in dollars, but the heavy tail is why a large n helps.

Bootstrap the mean of spending

set.seed(9)
n <- nrow(browser)
boot_spend <- rep(0, B)
for (b in 1:B) {
  idx <- sample.int(n, replace = TRUE)
  boot_spend[b] <- mean(browser$spend[idx])
}
se_formula <- sd(browser$spend) / sqrt(n)
c(mean = mean(browser$spend), boot_se = sd(boot_spend), formula_se = se_formula)

##       mean    boot_se formula_se 
## 1946.43930   79.47364   80.38610

The same loop as before, with browser$spend in place of the simulated sample. The two standard errors differ by about a dollar.

The confidence intervals

xbar <- mean(browser$spend)
dev <- abs(boot_spend - xbar)
q95 <- unname(quantile(dev, 0.95))
ci <- rbind(
  boot_normal = xbar + c(-2, 2) * sd(boot_spend),
  boot_quantile = xbar + c(-1, 1) * q95,
  formula = xbar + c(-2, 2) * se_formula)
colnames(ci) <- c("lower", "upper")
round(ci)

##               lower upper
## boot_normal    1787  2105
## boot_quantile  1792  2101
## formula        1786  2107

The average household spends about $1,950 a year online, and the population average is plausibly between $1,790 and $2,105.

A general bootstrap function

The loop changes only in the statistic it computes, so wrap it:

boot_stat <- function(x, stat, B = 1000) {
  replicate(B, stat(x[sample.int(length(x), replace = TRUE)]))
}

The median of spending has no simple standard-error formula: it involves the density of spending at the median, which we do not know. The bootstrap needs one line:

set.seed(10)
boot_med <- boot_stat(browser$spend, median)
med <- median(browser$spend)
c(median = med, boot_se = sd(boot_med))

##     median    boot_se 
## 510.000000   9.893141

med + c(-1, 1) * unname(quantile(abs(boot_med - med), 0.95))

## [1] 491 529

Beyond one variable: a ratio of two averages

Households with broadband spend on average about twice as much as those without. To bootstrap the ratio, resample rows, so that each household keeps its own spend and broadband:

ratio <- function(d) mean(d$spend[d$broadband == 1]) / mean(d$spend[d$broadband == 0])
set.seed(11)
boot_ratio <- replicate(B, ratio(browser[sample.int(n, replace = TRUE), ]))
c(ratio = ratio(browser), boot_se = sd(boot_ratio))

##     ratio   boot_se 
## 2.0399882 0.1884974

ratio(browser) + c(-1, 1) * unname(quantile(abs(boot_ratio - ratio(browser)), 0.95))

## [1] 1.663042 2.416935

A delta-method formula for the standard error of a ratio exists, but it is tedious to derive. The bootstrap is three lines. Resampling rows is also how the bootstrap handles regression coefficients (appendix).

When can we use the bootstrap?

• The bootstrap typically works well where the usual interval \bar{X} \pm 2 \cdot \text{SE} also works well: smooth estimators and a large sample. Otherwise it typically fails too.
• So why use it? Deriving standard-error formulas is often harder than coding up the bootstrap.
• Later in the course, LASSO and other high-dimensional methods involve penalties for model complexity and data-driven model selection. Standard CIs fail there, and so does the plain bootstrap unless it is adjusted.
• The bootstrap remains useful in many settings, and it is a gateway to other simulation-based tools: cross-validation (lecture 12) and bagging (lecture 14).

The appendix shows a textbook failure: the sample maximum, where the bootstrap distribution puts about two thirds of its mass on a single value.

Checkpoint: design a bootstrap

You want a 95% interval for (a) the share of households spending more than $5,000 a year and (b) the difference in median spending between households with and without children.

1. What is the statistic function in each case?
2. Which of the two has an easy formula for its standard error?
3. How would you check that B = 1{,}000 replications are enough?

Answers: (1) mean(spend > 5000) applied to a resampled vector; for (b), a function of a resampled data frame that returns median(spend[anychildren == 1]) - median(spend[anychildren == 0]). (2) The share: a proportion p has standard error sqrt(p(1 - p)/n). Medians have no simple formula. (3) Re-run with another seed or with B = 5000 and compare SE_boot; if the answer moves at the precision you report, increase B. The appendix shows that the Monte Carlo error is small at B = 1000.

Sampling variation is the price of using a sample

Two kinds of objects. Sample objects (averages, fractions, histograms, estimates) approximate population objects (expectations, probabilities, densities, parameters). Know which is which.

Two theorems. The law of large numbers says \bar{X} \approx \mathbf{E}[X] because \text{Var}(\bar{X}) = \text{Var}(X)/n shrinks. The central limit theorem adds the shape: approximately normal.

One sample, one interval. \text{SE}(\bar{X}) = \widehat{\text{SD}}(X)/\sqrt{n}, and \bar{X} \pm 2 \cdot \text{SE} covers the truth in about 95% of samples once n is moderate.

Resample when there is no formula. The bootstrap draws with replacement from the sample, recomputes the estimator, and reads the standard error and quantiles off the result. It works where the usual interval works.

Reading and data sources

• Taddy (2019), Business Data Science, McGraw-Hill, chapter 1: sampling, the bootstrap, and the web-browser data.
• Efron and Tibshirani (1993), An Introduction to the Bootstrap, Chapman and Hall: the standard reference.
• R documentation: Distributions, sample, replicate, set.seed.
• The boot package (Canty and Ripley) implements the bootstrap and several confidence-interval methods.

Additional details

---

More bootstrap intervals, the boot package, a regression preview, and two failure cases

Other bootstrap confidence intervals

The bootstrap distribution supports several intervals. For the simulated sample of 500 customers:

q <- quantile(boot_means, c(0.025, 0.975))
q95_sim <- unname(quantile(abs(boot_means - mean(spend)), 0.95))
ci <- rbind(boot_normal = mean(spend) + c(-2, 2) * sd(boot_means),
            symmetric_quantile = mean(spend) + c(-1, 1) * q95_sim,
            percentile = q,
            basic = 2 * mean(spend) - rev(q))
colnames(ci) <- c("lower", "upper")
ci

##                        lower    upper
## boot_normal        0.9647421 1.152168
## symmetric_quantile 0.9689625 1.147948
## percentile         0.9723376 1.152070
## basic              0.9648407 1.144573

• Percentile: the 2.5% and 97.5% quantiles of \bar{X}_b. Simple, but it does not correct for bias.
• Basic: the percentile interval reflected around \bar{X}, from 2\bar{X} - q_{0.975} to 2\bar{X} - q_{0.025}.
• With a symmetric, nearly unbiased bootstrap distribution the four intervals almost coincide, as here.

The boot package

boot() takes the data and a function of the data and an index vector; boot.ci() computes the intervals.

library(boot)
set.seed(9)
b <- boot(browser$spend, statistic = function(x, i) mean(x[i]), R = 1000)
c(original = b$t0, boot_se = sd(b$t))

##   original    boot_se 
## 1946.43930   81.84329

boot.ci(b, type = c("norm", "perc", "basic"))

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = b, type = c("norm", "perc", "basic"))
## 
## Intervals : 
## Level      Normal              Basic              Percentile     
## 95%   (1786, 2107 )   (1770, 2091 )   (1802, 2123 )  
## Calculations and Intervals on Original Scale

Preview: bootstrapping a regression coefficient

Lecture 10 regresses log spending on broadband access. The pairs bootstrap resamples rows and refits the model:

fit <- lm(log(spend) ~ broadband, data = browser)
coef(summary(fit))[, c("Estimate", "Std. Error")]

##             Estimate Std. Error
## (Intercept) 5.732062 0.03957177
## broadband   0.555759 0.04356445

set.seed(12)
boot_beta <- replicate(B, {
  d <- browser[sample.int(nrow(browser), replace = TRUE), ]
  coef(lm(log(spend) ~ broadband, data = d))["broadband"]
})
c(formula_se = coef(summary(fit))["broadband", "Std. Error"], boot_se = sd(boot_beta))

## formula_se    boot_se 
## 0.04356445 0.04427595

The regression standard error assumes errors with constant variance; the pairs bootstrap does not. Here the two agree.

How many bootstrap replications?

B controls only the Monte Carlo error of the bootstrap approximation. Re-run the mean of spending with different B:

set.seed(14)
sapply(c(B100 = 100, B1000 = 1000, B10000 = 10000),
       function(B) sd(boot_stat(browser$spend, mean, B)))

##     B100    B1000   B10000 
## 84.27191 81.58051 80.48366

A rule of thumb: the Monte Carlo standard error of \text{SE}_{\text{boot}} is roughly \text{SE}_{\text{boot}}/\sqrt{2B}, about 2% of \text{SE}_{\text{boot}} at B = 1{,}000.

Quantile intervals need more replications than standard errors, because they depend on the tails of the bootstrap distribution. B = 1{,}000 or more is common for a 95% interval.

Where the bootstrap fails: the sample maximum

The maximum of m draws from \text{Uniform}(0, 1) estimates the upper end point 1. It is not a smooth function of the data, and the bootstrap breaks:

set.seed(13)
m <- 100
u <- runif(m)
boot_max <- replicate(B, max(u[sample.int(m, replace = TRUE)]))
true_max <- replicate(B, max(runif(m)))
c(sample_max = max(u), share_at_max = mean(boot_max == max(u)),
  boot_sd = sd(boot_max), true_sd = sd(true_max))

##   sample_max share_at_max      boot_sd      true_sd 
##  0.999366581  0.647000000  0.012257865  0.009201846

About 65% of the bootstrap maxima equal the sample maximum itself, close to the 63% chance that the largest observation is drawn at least once. The bootstrap spread is a third too large, and its shape is wrong.

The bootstrap distribution of the maximum is lumpy

The true sampling distribution is smooth and skewed towards 1; the bootstrap version is a spike. Estimators that select among models, such as LASSO, are non-smooth in a related way.

Where the law of large numbers fails: infinite variance

Both theorems need moments. The Cauchy distribution has no expectation, and averages of even 100,000 draws do not settle down:

set.seed(15)
replicate(5, mean(rcauchy(1e5)))

## [1]  0.4744233  5.5317096 -2.2550010  6.8630508 -0.5150997

Heavy tails in economic data, such as firm sizes, city populations or wealth, are usually less extreme than this, but they slow the CLT down: compare the coverage at n = 10 earlier. Log transformations, as with log(spend), tame the tails.

Install the required packages

Run once in your R environment if needed:

install.packages(c("tidyverse", "patchwork", "boot"))

boot ships with R as a recommended package, so it is usually installed already. The web-browser data are in lectures/data/ and the schematic figures in lectures/images/, so the deck renders without network access.

Return to the question

Where did the estimate come from, and how much would it move in another sample?

The standard error answers with a formula when one exists. The bootstrap answers by resampling when one does not.

Return to the main takeaway