5  Sampling & Resampling

---

We rarely observe an entire population; instead we work with a sample and ask what its statistics tell us about the population that produced it.

This chapter asks three questions.

• How much does a statistic change from one sample to the next?
  • Simple random sample: this is the baseline sampling scheme, where every subset of size \(n\) has the same probability of being chosen.
  • Sampling distribution: this shows how a statistic, such as the mean, varies across samples from one population.
• What happens to the sample mean and other statistics as samples grow?
  • Central Limit Theorem: this says the sampling distribution of \(\hat{m}\) becomes approximately Normal.
  • Law of Large Numbers: this says the sampling distribution of \(\hat{m}\) tightens around the population mean.
  • Fundamental Theorem of Statistics: this says the ECDF converges to the true CDF, so statistics built from it, such as quantiles, converge too.
• How sure are we?
  • Bootstrap and jackknife: with only one sample, these resample the data to approximate the sampling distribution.
  • Standard error: this is the standard deviation of a sampling distribution, which reports precision as one number.

5.1 Sampling

In practice we almost never observe the entire group we care about; we observe a slice and have to make judgments from that slice.

ImportantKey Definition

A population is the entire collection of units of interest.

A sample is a subset of the population that we actually observe.

A simple random sample is a sample where every subset of size \(n\) has the same probability of being chosen.

The simple random sample is useful as the default baseline. Every other sampling scheme is described by how it departs from it. The mathematical results in the rest of the chapter also assume a simple random sample.

NoteMust Know

List the samples of size two from the small population \(\{1, 2, 3, 4\}\).

#All possible samples of two from a bag of numbers {1, 2, 3, 4} without replacement
# {1, 2} {1, 3} {1, 4}
#        {2, 3} {2, 4}
#               {3, 4}

#combn() lists every possible sample, one sample per column
population_values <- c(1, 2, 3, 4) # population
combn(population_values, 2)   
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    1    1    2    2    3
## [2,]    2    3    4    3    4    4

# Simple random sample (no duplicates, equal probability)
sample(population_values, 2, replace=FALSE) #sample
## [1] 1 3

Often, we think of the population as being infinitely large. This is an approximation that makes mathematical and computational work much simpler.

Intuition for infinite populations: imagine drawing names from a giant urn. If the urn has only \(10\) names, then removing one name slightly changes the composition of the urn, and the probabilities shift for the next name you draw. Now imagine the urn has \(100\) billion names, so that removing one makes no noticeable difference. We can pretend the composition never changes: each draw is essentially identical and independent (iid). We can actually guarantee the names are iid by putting any names drawn back into the urn (sampling with replacement).

NoteMust Know

Sampling with replacement from the small population \(\{1, 2, 3, 4\}\) allows duplicates, so every draw has the same four possible values.

#All possible samples of two from a bag of numbers {1, 2, 3, 4} with replacement
#{1, 1} {1, 2} {1, 3} {1, 4}
#       {2, 2} {2, 3} {2, 4}
#              {3, 3} {3, 4}
#                     {4, 4}

# Simple random sample (duplicates, equal probability)
population_values <- c(1, 2, 3, 4) # population
sample(population_values, 2, replace=TRUE)
## [1] 1 2

Population Theory

In practice, we never see the whole population, so we cannot check how close a sample statistic is to the truth. To make such checks possible, this chapter pretends that one dataset is the entire population. It is really a sample itself, but pretending lets us compute the true population mean and compare it with what one sample tells us.

## MAIN EXAMPLE
# The data come from the Ecdat package, which you need to install once.
# install.packages('Ecdat') # only once
data(Wages1, package='Ecdat')

# Pretend these 3294 hourly wages are the whole population
population_wages <- Wages1[, 'wage'] # hourly wage in 1980 dollars

# One simple random sample of 50 workers
# Each time you run, you get a different sample of workers.
x_sim <- sample(population_wages, 50, replace=FALSE)

Permutations and Combinations

For small populations, it can be useful to count how many samples there are. Factorials are used for counting permutations: different ways of rearranging \(N\) distinct objects into a sequence. There are \(N! = N \times (N-1) \times \cdots \times 2 \times 1\) orderings, read “\(N\) factorial”. When we select \(n\) items from \(N\) and the order matters, there are \(N \times (N-1) \times \cdots \times (N-n+1) = \frac{N!}{(N-n)!}\) permutations.

When the order does not matter, we divide out the \(n!\) rearrangements and count combinations. The binomial coefficient counts these subsets of \(n\) elements from a set with \(N\) elements, and is mathematically defined as \(\tbinom{N}{n}=\frac{N!}{n!(N-n)!}\).

NoteMust Know

How many ways are there to order the numbers \(\{1, 2, 3\}\)?

# All Permutations
#(1, 2, 3) (1, 3, 2)
#(2, 1, 3) (2, 3, 1)
#(3, 2, 1) (3, 1, 2)
factorial(3)
## [1] 6

How many possible samples of two are there from a population with this data: \(\{7, 10, 122, 55\}\)?

#Ways to draw n=2 from a set with N=4
#{7, 10}  {7, 122}   {7, 55}
#        {10, 122}  {10, 55}
#                  {122, 55}

#combn() lists every possible sample, one sample per column
combn(c(7, 10, 122, 55), 2)   
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    7    7    7   10   10  122
## [2,]   10  122   55  122   55   55

# choose() counts them
choose(4, 2)
## [1] 6

The answer depends only on \(N=4\) and \(n=2\), not on the values themselves.

TipTest Yourself

A small firm has four workers, labeled \(\{1, 2, 3, 4\}\), and must name one manager and one assistant manager. We write an ordered pair in parentheses, \((\text{manager}, \text{assistant})\), and an unordered set in braces. Order matters here: \((1, 2)\) makes worker \(1\) the manager, while \((2, 1)\) makes worker \(2\) the manager. There are \(4\) choices for manager and then \(3\) for assistant, so there are \(4 \times 3 = 12\) permutations.

#All (manager, assistant) pairs from {1, 2, 3, 4}
#(1, 2) (1, 3) (1, 4)
#(2, 1) (2, 3) (2, 4)
#(3, 1) (3, 2) (3, 4)
#(4, 1) (4, 2) (4, 3)

# N x (N-1) x ... x (N-n+1) with N=4 and n=2
4 * 3
## [1] 12

# Same count written with factorials: N!/(N-n)!
factorial(4) / factorial(4 - 2)
## [1] 12

Now suppose the firm instead picks two workers for a committee with no roles, so order does not matter. Each committee appears twice in the list above: \((1, 2)\) and \((2, 1)\) are both the committee \(\{1, 2\}\). Dividing the \(12\) permutations by the \(2! = 2\) orderings of each pair leaves \(\binom{4}{2} = 6\) combinations.

# Permutations divided by the orderings within each pair
(4 * 3) / factorial(2)
## [1] 6
choose(4, 2)
## [1] 6

TipTest Yourself

License plates. A plate has the format AB1 23C: three letter slots and three number slots, with no letter and no number repeated. The letters can be filled \(26 \times 25 \times 24\) ways and the numbers \(10 \times 9 \times 8\) ways, so by the multiplication rule there are \((26 \times 25 \times 24)(10 \times 9 \times 8)\) possible plates.

Bank-account sample. From \(20\) accounts we draw a sample of \(3\), and the order of selection does not matter. The number of possible samples is \(\binom{20}{3}\).

HR committee. In the morning, \(2\) of \(25\) employees are chosen for distinct roles (manager, then assistant manager), so order matters: \(25 \times 24\) ways. In the afternoon, \(5\) of \(23\) employees are chosen for termination, where order does not matter: \(\binom{23}{5}\) ways. The two decisions are independent, so the multiplication rule gives \((25 \times 24)\binom{23}{5}\) outcomes for the day.

# License plates: 3 distinct letters, 3 distinct numbers
(26 * 25 * 24) * (10 * 9 * 8)
## [1] 11232000

# Bank-account sample of 3 from 20 (order irrelevant)
choose(20, 3)
## [1] 1140

# HR committee: morning ordered, afternoon unordered
(25 * 24) * choose(23, 5)
## [1] 20189400

5.2 Sampling Distributions

A statistic computed on one sample is just one number; we want to know how much that number would change if we drew a different sample.

ImportantKey Definition

The sampling distribution of a statistic is the probability distribution of the statistic’s value across all possible samples of size \(n\) from the same population. It tells us how much a statistic varies from one sample to the next.

The sampling distribution is useful for quantifying uncertainty: its spread shows how much the statistic typically moves from sample to sample, and its center shows where it is targeted on average. For example, the sampling distribution of the mean shows how \(\hat{m}\) varies from sample to sample; it can also be called the probability distribution of the sample mean.

NoteMust Know

Given ages for population of \(4\) students, compute the sampling distribution for the mean with samples of \(n=2\).

population_values <- c(18, 20, 22, 24) # Ages for student population
# six possible samples
m1_hat <- mean( population_values[c(1, 2)] ) #{1, 2}
m2_hat <- mean( population_values[c(1, 3)] ) #{1, 3}
m3_hat <- mean( population_values[c(1, 4)] ) #{1, 4}
m4_hat <- mean( population_values[c(2, 3)] ) #{2, 3}
m5_hat <- mean( population_values[c(2, 4)] ) #{2, 4}
m6_hat <- mean( population_values[c(3, 4)] ) #{3, 4}
# sampling distribution
sample_means <- c(m1_hat, m2_hat, m3_hat, m4_hat, m5_hat, m6_hat)
hist(sample_means,
    freq=FALSE, breaks=100,
    main=NA, border=NA)

TipTest Yourself

Now compute the sampling distribution for the median with samples of \(n=3\). There are four possible samples, \(\{18,20,22\}\), \(\{18,20,24\}\), \(\{18,22,24\}\), and \(\{20,22,24\}\), with medians \(20\), \(20\), \(22\), and \(22\).

population_values <- c(18, 20, 22, 24) # Ages for student population
# four possible samples
med1_hat <- median( population_values[c(1, 2, 3)] )
med2_hat <- median( population_values[c(1, 2, 4)] )
med3_hat <- median( population_values[c(1, 3, 4)] )
med4_hat <- median( population_values[c(2, 3, 4)] )
sample_medians <- c(med1_hat, med2_hat, med3_hat, med4_hat)
sample_medians
## [1] 20 20 22 22

NoteMust Know

To consider an infinite population, draw from a uniform distribution instead of a list of values.

# Three Sample Example from infinite population
x1_sim <- runif(100)
m1_hat <- mean(x1_sim)
x2_sim <- runif(100)
m2_hat <- mean(x2_sim)
x3_sim <- runif(100)
m3_hat <- mean(x3_sim)
sample_means <- c(m1_hat, m2_hat, m3_hat)
sample_means
## [1] 0.5058161 0.5143354 0.4557783

# An Equivalent Approach: fill vector in a loop
sample_means <- rep(NA, 3)
for(i in seq(sample_means)){
    x_sim <- runif(100)
    m_hat <- mean(x_sim)
    sample_means[i] <- m_hat
}
sample_means
## [1] 0.5004781 0.5209211 0.5254840

For more on loops, see Loops.

Each sample has its own histogram and its own mean.

# Three Sample Example w/ Visual
par(mfrow=c(1, 3))
for(b in 1:3){
    x_sim <- runif(100)
    m_hat <-  mean(x_sim)
    hist(x_sim,
        breaks=seq(0, 1, by=.1), #for comparability
        freq=FALSE, main=NA, border=NA)
    abline(v=m_hat, col=rgb(1, 0, 0, .8), lwd=2)
    title(paste0('mean= ', round(m_hat, 2)),  font.main=1)
}

Expanding the loop to many samples gives the sampling distribution of the mean.

# Many sample example
sample_means <- rep(NA, 500)
for(i in seq_along(sample_means)){
    x_sim <- runif(1000)
    m_hat <- mean(x_sim)
    sample_means[i] <- m_hat
}
hist(sample_means,
    breaks=seq(0.45, 0.55, by=.001),
    border=NA, freq=FALSE,
    col=rgb(1, 0, 0, .5), # red marks the distribution of a statistic, not data
    xlab=expression(hat(m)),
    main=NA)
title('Sampling Distribution of the mean', font.main=1)

Now draw three samples from the pretend population. Each sample mean (red) lands near, but not on, the population mean (blue dashed).

## MAIN EXAMPLE
# Population mean
population_wage_mean <- mean(population_wages)

# Three samples of workers, each with its own mean
par(mfrow=c(1, 3))
for(b in 1:3){
    x_sim <- sample(population_wages, 50, replace=FALSE)
    m_hat <- mean(x_sim)
    hist(x_sim,
        breaks=seq(0, 40, by=2), #for comparability
        freq=FALSE, main=NA, border=NA,
        xlab='Hourly wage (1980 dollars)')
    rug(x_sim, col=grey(0, .5))
    abline(v=m_hat, col=rgb(1, 0, 0, .8), lwd=2)
    abline(v=population_wage_mean, col=rgb(0, 0, 1, .8), lty=2)
    title(paste0('mean= ', round(m_hat, 2)), font.main=1)
}

Repeating this for many samples gives the sampling distribution of the mean. We color this histogram red to make clear we are looking at the distribution of a statistic, while the grey histograms above show data. In what follows, we typically use the default grey for both sampling distributions and data.

## MAIN EXAMPLE
# Many samples of workers
wage_means <- rep(NA, 500)
for(i in seq_along(wage_means)){
    x_sim <- sample(population_wages, 50, replace=FALSE)
    m_hat <- mean(x_sim)
    wage_means[i] <- m_hat
}
hist(wage_means,
    breaks=30,
    border=NA, freq=FALSE,
    col=rgb(1, 0, 0, .5), # red marks the distribution of a statistic, not data
    xlab=expression(hat(m)),
    main=NA)
abline(v=population_wage_mean, col=rgb(0, 0, 1, .8), lty=2)
title('Sampling Distribution of the mean wage', font.main=1)

In this figure, you see one of the most profound results known in statistics: the sampling distribution of the mean is approximately Normal around the true mean.

Central Limit Theorem (CLT)

ImportantKey Definition

The Central Limit Theorem (CLT) says that as \(n\) grows, the sampling distribution of \(\hat{m}\) is approximately Normal, regardless of the shape of the underlying population (under mild conditions on variance).

The CLT is useful because it lets us treat sample means as approximately Normal, even when we know little about the population distribution. For this reason, many statistical procedures, such as intervals and hypothesis tests, are built on Normal calculations. There are different variants of the theorem, but all say some version of “the sampling distribution of the mean is approximately Normal”; the histogram shown above is one such example.

NoteMust Know

The means of uniform samples are approximately Normal.

sample_means <- rep(NA, 500)
for(i in seq_along(sample_means)){
    x_sim <- runif(1000)
    m_hat <- mean(x_sim)
    sample_means[i] <- m_hat
}
hist(sample_means,
    breaks=seq(0.45, 0.55, by=.001),
    border=NA, freq=FALSE,
    xlab=expression(hat(m)),
    main=NA)
title('Sampling Distribution of the mean', font.main=1)

# Approximately normal?
m_sampling_hat <- mean(sample_means)
s_sampling_hat <- sd(sample_means)
x <- seq(0.1, 0.9, by=0.001)
fx <- dnorm(x, m_sampling_hat, s_sampling_hat)
lines(x, fx, col=rgb(0, 0, 1, .8), lty=2)

A strongly skewed population is a harder case. The figure below compares the population with the means of many samples drawn from it.

## MAIN EXAMPLE
par(mfrow=c(1, 2))

# Population: skewed
hist(population_wages,
    breaks=seq(0, 40, by=1),
    border=NA, freq=FALSE,
    xlab='Hourly wage (1980 dollars)',
    main=NA)
rug(population_wages, col=grey(0, .5))
abline(v=population_wage_mean, col=rgb(0, 0, 1, .8), lty=2)
title('Population', font.main=1)
# Wages are strongly skewed.
# Most workers earn a moderate wage and a few earn much more.

# Sampling distribution: approximately Normal?
hist(wage_means,
    breaks=30,
    border=NA, freq=FALSE,
    xlab=expression(hat(m)),
    main=NA)
m_sampling_hat <- mean(wage_means)
s_sampling_hat <- sd(wage_means)
x <- seq(min(wage_means), max(wage_means), length.out=200)
fx <- dnorm(x, m_sampling_hat, s_sampling_hat)
lines(x, fx, col=rgb(0, 0, 1, .8), lty=2)
title('Sample means, n=50', font.main=1)

# The means of samples of 50 workers are not skewed in the same way.

Many statistics have an approximately Normal sampling distribution.

TipTest Yourself

For an example with another statistic, let’s examine the sampling distribution of the standard deviation.

# CLT example of the 'sd' statistic
sample_sds <- rep(NA, 1000)
for(i in seq_along(sample_sds)){
    x_sim <- runif(100) # same distribution
    s_hat <- sqrt(mean((x_sim - mean(x_sim))^2)) # different statistic, which divides by n, unlike sd()
    sample_sds[i] <- s_hat
}
hist(sample_sds,
    breaks=seq(0.2, 0.4, by=.01),
    border=NA, freq=FALSE,
    xlab=expression(hat(s)),
    main=NA)
title('Sampling Distribution of the sd', font.main=1)

# Approximately normal?
m_sampling_hat <- mean(sample_sds)
s_sampling_hat <- sd(sample_sds)
x <- seq(0.1, 0.9, by=0.001)
fx <- dnorm(x, m_sampling_hat, s_sampling_hat)
lines(x, fx, col=rgb(0, 0, 1, .8), lty=2)

# try another random variable, such as rexp(100) instead of runif(100)

It is beyond this book to prove this result mathematically, but you should know that not all sampling distributions are Normal. The CLT approximation is better for “large \(n\)” datasets with “well behaved” variances. The CLT also does not apply to “extreme” statistics.

NoteMust Know

For an example of “extreme” statistics, examine the sampling distribution of min and max statistics.

# Create 300 samples, each with 1000 random uniform variables
X_samples <- matrix(nrow=300, ncol=1000)
for(i in seq(1, nrow(X_samples))){
    X_samples[i, ] <- runif(1000)
}
# Each row is a new sample
length(X_samples[1, ])
## [1] 1000

# Compute min and max for each sample
x_mins <- apply(X_samples, 1, quantile, probs=0)
x_maxs <- apply(X_samples, 1, quantile, probs=1)

# Plot the sampling distributions of the min and max
# Neither one looks normal.
par(mfrow=c(1, 2))
hist(x_mins, breaks=100, main=NA,
    xlab=expression(hat(q)[0]), border=NA, freq=FALSE)
title('Min', font.main=1)
hist(x_maxs, breaks=100, main=NA,
    xlab=expression(hat(q)[1]), border=NA, freq=FALSE)
title('Max', font.main=1)
title('Sampling Distributions', outer=TRUE, line=-1, adj=0, font.main=1)

TipTest Yourself

Explore the sampling distribution of another statistic, such as the range.

my_function <- function(x){ diff(range(x)) }
x_ranges <- apply(X_samples, 1, my_function)
hist(x_ranges, breaks=100, main=NA,
    xlab='Sample Range', border=NA, freq=FALSE)

TipTest Yourself

Here is an example where variance is not “well behaved” .

sample_means <- rep(NA, 999)
for(i in seq_along(sample_means)){
    x_sim <- rcauchy(1000, 0, 10)
    m_hat <- mean(x_sim)
    sample_means[i] <- m_hat
}
hist(sample_means, breaks=100,
    main='',
    border=NA, freq=FALSE) # Tails look too 'fat'

Law of Large Numbers (LLN)

The CLT tells us the sample mean is Normally distributed around the true mean. Another profound result is that the sampling distribution tightens with more observations.

ImportantKey Definition

The Law of Large Numbers (LLN) says the sample mean concentrates around the population mean as \(n\) grows: with more data, the sampling distribution of \(\hat{m}\) tightens around \(\mathbb{E}[X_{i}]\).

The LLN is useful as a guarantee that averaging works: even when individual observations are noisy, the sample mean of a large sample is close to the population mean. There are different variants of the theorem, such as weak vs. strong and with vs. without finite variance. They all say some version of “the sample mean becomes more tightly centered around the true mean as we get more data”.

NoteMust Know

Uniform draws on \([0,1]\) have a population mean of \(0.5\). Notice where the sampling distribution is centered

sample_means <- rep(NA, 500)
for(i in seq_along(sample_means)){
    x_sim <- runif(1000)
    sample_means[i] <- mean(x_sim)
}
m_sampling_hat <- mean(sample_means)
round(m_sampling_hat, 3)
## [1] 0.5

and more tightly centered with more data

par(mfrow=c(1, 3))
for(n in c(5, 50, 500)){
    sample_means_n <- rep(NA, 299)
    for(i in seq_along(sample_means_n)){
        x_sim <- runif(n)
        m_hat <- mean(x_sim)
        sample_means_n[i] <- m_hat
    }
    hist(sample_means_n,
        breaks=seq(0, 1, by=.01),
        border=NA, freq=FALSE,
        xlab=expression(hat(m)),
        main=NA)
    title(paste0('n=', n), font.main=1)
}

In the pretend population, the sample means are centered on the population mean

round(mean(wage_means), 3) # center of the sampling distribution
## [1] 5.746
round(population_wage_mean, 3) # population mean
## [1] 5.758

and are more tightly centered with larger samples.

## MAIN EXAMPLE
par(mfrow=c(1, 3))
for(n in c(5, 50, 500)){
    wage_means_n <- rep(NA, 299)
    for(i in seq_along(wage_means_n)){
        x_sim <- sample(population_wages, n, replace=FALSE)
        m_hat <- mean(x_sim)
        wage_means_n[i] <- m_hat
    }
    hist(wage_means_n,
        breaks=seq(0, 20, by=.25),
        border=NA, freq=FALSE,
        xlab=expression(hat(m)),
        main=NA)
    abline(v=population_wage_mean, col=rgb(0, 0, 1, .8), lty=2)
    title(paste0('n=', n), font.main=1)
}

TipTest Yourself

Plot the variability of the sample mean as a function of sample size

n_seq <- seq(1, 40)
sd_seq <- rep(NA, length(n_seq))
for(n in seq_along(sd_seq)){
    sample_means_n <- rep(NA, 499)
    for(i in seq_along(sample_means_n)){
        x_sim <- runif(n)
        m_hat <- mean(x_sim)
        sample_means_n[i] <- m_hat
    }
    sd_seq[n] <- sd(sample_means_n)
}
plot(n_seq, sd_seq, pch=16, col=grey(0, 0.5),
    xlab='n', ylab='sd of sample means', main=NA)

Here is the intuition for estimating the mean weight of an apple:

• With \(n=1\) apple, your estimate depends entirely on that one draw. If it happens to be unusually large or small, your estimate can be far off.
• With \(n=2\) apples, the estimate averages out their idiosyncrasies. An unusually heavy apple can be balanced by a lighter one, lowering how far off you can be. You are less likely to get two extreme values than just one.
• With \(n=100\) apples, individual apples barely move the needle. The average becomes stable.

NoteMust Know

Here is an intuitive example, from a small discrete population. Notice the extreme values

population_values <- c(18, 20, 22, 24) #student ages (population values)
mean(population_values) # Population mean
## [1] 21

# six possible samples of size 2
m1_hat <- mean( population_values[c(1, 2)] ) #{1, 2}
m2_hat <- mean( population_values[c(1, 3)] ) #{1, 3}
m3_hat <- mean( population_values[c(1, 4)] ) #{1, 4}
m4_hat <- mean( population_values[c(2, 3)] ) #{2, 3}
m5_hat <- mean( population_values[c(2, 4)] ) #{2, 4}
m6_hat <- mean( population_values[c(3, 4)] ) #{3, 4}
means_2 <- c(m1_hat, m2_hat, m3_hat, m4_hat, m5_hat, m6_hat)
sort(means_2)
## [1] 19 20 21 21 22 23

# four possible samples of size 3
m1_hat <- mean( population_values[c(1, 2, 3)] )
m2_hat <- mean( population_values[c(1, 2, 4)] )
m3_hat <- mean( population_values[c(1, 3, 4)] )
m4_hat <- mean( population_values[c(2, 3, 4)] )
means_3 <- c(m1_hat, m2_hat, m3_hat, m4_hat)
sort(means_3)
## [1] 20.00000 20.66667 21.33333 22.00000

The Fundamental Theorem of Statistics

The Law of Large Numbers generalizes to many other statistics, like the median or sd from Summary Statistics.¹ Intuitively, recall that the proportion of a sample smaller than some value \(x\) can be computed as a type of mean: \(\hat{F}(x)=\frac{1}{n}\sum_{i} \mathbf{1}(\hat{x}_{i} \leq x)\), and so the LLN applies here too. As such, the LLN also applies to the associated quantile \(\hat{F}^{-1}\).

NoteMust Know

Consider a standard uniform random variable. Here is a sampling distribution for the proportion \(\leq 3/4\), for three different sample sizes.

par(mfrow=c(1, 3))
for(n in c(5, 50, 500)){
    sample_props_n <- rep(NA, 299)
    for(i in seq_along(sample_props_n)){
        x_sim <- runif(n)
        p3_hat <- mean(x_sim<=0.75)
        sample_props_n[i] <- p3_hat
    }
    hist(sample_props_n,
        breaks=seq(0, 1, by=.01),
        border=NA, freq=FALSE,
        xlab='Sample Proportion',
        main=NA)
    title(paste0('n=', n), font.main=1)
}

Here is a sampling distribution for the upper quartile, for three different sample sizes

par(mfrow=c(1, 3))
for(n in c(5, 50, 500)){
    sample_quants_n <- rep(NA, 299)
    for(i in seq_along(sample_quants_n)){
        x_sim <- runif(n)
        q_hat <- quantile(x_sim, probs=0.75)
        sample_quants_n[i] <- q_hat
    }
    hist(sample_quants_n,
        breaks=seq(0, 1, by=.01),
        border=NA, freq=FALSE,
        xlab='Sample Quantile',
        main=NA)
    title(paste0('n=', n), font.main=1)
}

In fact, the entire empirical distribution converges, not just individual statistics like the mean, median, or a proportion.

ImportantKey Definition

The Fundamental Theorem of Statistics says that the ECDF converges to the true CDF as \(n\) grows: the two curves get uniformly close across the entire range of values.

This theorem is useful because many statistics, such as means, quantiles, and proportions, are just functions of the ECDF. When the ECDF converges to the CDF, these statistics converge too. The LLN is one example.

NoteMust Know

par(mfrow = c(1, 3))
for (n in c(50, 500, 5000)) {
  x_sim <- runif(n)
  F_hat <- ecdf(x_sim)
  plot(F_hat, main=NA)
  title(paste0('n=', n), font.main=1)
}

In the pretend population, the ECDF of a sample (red) gets closer to the population CDF (black) as the sample grows.

## MAIN EXAMPLE
F_population <- ecdf(population_wages)
par(mfrow=c(1, 3))
for(n in c(5, 50, 500)){
    x_sim <- sample(population_wages, n, replace=FALSE)
    F_hat <- ecdf(x_sim)
    plot(F_population, do.points=FALSE,
        xlim=c(0, 20), xlab='Hourly wage (1980 dollars)', main=NA)
    plot(F_hat, do.points=FALSE, add=TRUE, col=rgb(1, 0, 0, .8))
    title(paste0('n=', n), font.main=1)
}

TipTest Yourself

Show the ECDF converging to the CDF for a small discrete population of five student ages \(\{18, 20, 22, 24, 26\}\). List all \(\binom{5}{n}\) possible samples of size \(n=2, 3, 4, 5\), and compare each sample’s \(\hat{F}\) (red) to the population \(F\) (black).

population_values <- c(18, 20, 22, 24, 26) #student ages (population values)
F_population <- ecdf(population_values)
F_population(population_values) # each age adds 1/5
## [1] 0.2 0.4 0.6 0.8 1.0

# All possible samples of each size, one sample per column
for(n in 2:5){
    print( combn(population_values, n) )
}
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,]   18   18   18   18   20   20   20   22   22    24
## [2,]   20   22   24   26   22   24   26   24   26    26
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,]   18   18   18   18   18   18   20   20   20    22
## [2,]   20   20   20   22   22   24   22   22   24    24
## [3,]   22   24   26   24   26   26   24   26   26    26
##      [,1] [,2] [,3] [,4] [,5]
## [1,]   18   18   18   18   20
## [2,]   20   20   20   22   22
## [3,]   22   22   24   24   24
## [4,]   24   26   26   26   26
##      [,1]
## [1,]   18
## [2,]   20
## [3,]   22
## [4,]   24
## [5,]   26

# Population CDF nudged up by 0.01, so it stays visible above the red lines
F_population_up <- stepfun(population_values, c(0, F_population(population_values)) + 0.01)

par(mfrow=c(1, 4))
for(n in 2:5){
    samples_n <- combn(population_values, n)
    plot(F_population_up, do.points=FALSE, verticals=FALSE,
        xlim=c(16, 28), xlab='age', ylab='Fn(x)', main=NA)
    max_gaps <- rep(NA, ncol(samples_n))
    for(s in seq(ncol(samples_n))){
        F_hat <- ecdf(samples_n[, s])
        plot(F_hat, do.points=FALSE, add=TRUE, col=rgb(1, 0, 0, .3), lwd=4)
        # Both curves only jump at the five ages, so the largest gap is at one of them
        max_gaps[s] <- max(abs(F_hat(population_values) - F_population(population_values)))
    }
    # Largest gap over all possible samples of size n
    title(paste0('n=', n, ', largest gap=', round(max(max_gaps), 2)), font.main=1)
}

With \(n=5\) the only possible sample is the whole population, so \(\hat{F}=F\) exactly.

5.3 Resampling Distributions

With the pretend population, we could draw as many samples as we liked. Often, we only have one sample, and no population to draw more samples from. To mimic this, we draw one sample as our lab data, and pretend it is all we observed. How then can we estimate the sampling distribution of a statistic?

## MAIN EXAMPLE
# Lab data: one random sample of 50 workers, as if it were all we observed
x_hat <- sample(population_wages, 50, replace=FALSE)
m_hat <- mean(x_hat)
m_hat
## [1] 5.634566
# The sample mean is about 5.63, a little below the population mean of 5.758.

We can resample our data. Hesterberg (2015) provides a nice illustration of the idea. The two most basic versions are the jackknife and the bootstrap, which are discussed below. The red histogram in the figure again shows the distribution of a statistic.

Note that we do not use the mean of the resampled statistics as a replacement for the original estimate. This is because the resampled distributions are centered at the observed statistic, not the population parameter. (The bootstrapped mean is centered at the sample mean, for example, not the population mean.) This means that we do not use resampling to improve on \(\hat{m}\). We use resampling to estimate sampling variability.

Bootstrap Distribution

A flexible approach builds many resamples of full size \(n\) by allowing duplicates.

ImportantKey Definition

The bootstrap resamples \(n\) observations from the original data with replacement, computes the statistic on each resample, and uses the spread of those values to approximate the sampling distribution (Efron 1979; Efron and Tibshirani 1993; Horowitz 2001, 2019).

The bootstrap is useful when there is no formula for the sampling distribution, which is common for medians, quantiles, ratios, or any custom statistic. The recipe is the same no matter which statistic we compute on each resample. Each bootstrap sample \(b=1\dots B\) produces one value of the statistic, and we typically repeat the resampling many times (e.g., \(B=9999\)) to estimate the sampling distribution.

NoteMust Know

Given the sample \(\{1,6,7,22\}\), compute the bootstrap estimate of the median with \(B=8\) resamples. Each resample draws \(n=4\) values from \(\{1,6,7,22\}\) with replacement, then takes the median. For instance, one resample might be \(\{6,6,7,22\}\): sorted, the two middle values are \(6\) and \(7\), so the median is \((6+7)/2=6.5\). Another might be \(\{1,1,7,7\}\), with median \((1+7)/2=4\). Repeating this \(B=8\) times gives the bootstrap distribution of the median.

x0_hat <- c(1, 6, 7, 22)
boot_medians <- rep(NA, 8)
for(b in seq_along(boot_medians)){
    x_boot <- sample(x0_hat, replace=TRUE)
    boot_medians[b] <- median(x_boot)
}
boot_medians
## [1]  7.0  3.5  7.0  7.0  6.5 14.5  6.5  1.0

## MAIN EXAMPLE
# Bootstrap estimates
boot_means <- rep(NA, 9999)
for(b in seq_along(boot_means)){
    x_boot <- sample(x_hat, replace=TRUE) # c.f. jackknife
    m_b_hat <- mean(x_boot)
    boot_means[b] <-m_b_hat
}

hist(boot_means, breaks=25,
    border=NA, freq=FALSE,
    main=NA, xlab=expression(hat(m)^group('(', b, ')')))

# The bootstrap means center on the sample mean of about 5.63.

TipTest Yourself

Repeat the bootstrap with the murder arrests in USArrests.

# These 50 US states are one sample, with no known population.
x2_hat <- USArrests[, 'Murder'] # murder arrests per 100,000
m2_hat <- mean(x2_hat)
m2_hat
## [1] 7.788

# Bootstrap estimates
boot_means2 <- rep(NA, 9999)
for(b in seq_along(boot_means2)){
    x_boot <- sample(x2_hat, replace=TRUE) # c.f. jackknife
    m_b_hat <- mean(x_boot)
    boot_means2[b] <- m_b_hat
}

hist(boot_means2, breaks=25,
    border=NA, freq=FALSE,
    main=NA, xlab=expression(hat(m)^group('(', b, ')')))

Why does this work? The sample: \(\{\hat{x}_{1}, \hat{x}_{2}, ... \hat{x}_{n}\}\) is drawn from a CDF \(F\). Each bootstrap sample: \(\{\hat{x}_{1}^{(b)}, \hat{x}_{2}^{(b)}, ... \hat{x}_{n}^{(b)}\}\) is drawn from the ECDF \(\hat{F}\). With \(\hat{F} \approx F\) by the Fundamental-Theorem above, each bootstrap sample is approximately a random sample. So when we compute a statistic on each bootstrap sample, we approximate the sampling distribution of the statistic.

TipTest Yourself

Here is an intuitive example with Bernoulli random variables (unfair coin flips)

# theoretical probabilities
x <- c(0, 1)
x_probs <- c(1/4, 3/4)
# sample draws
coin_hat <- sample(x, prob=x_probs, 1000, replace=TRUE)
F_hat <- ecdf(coin_hat)

x_probs_boot <- c(F_hat(0), 1-F_hat(0))
x_probs_boot # approximately the theoretical value
## [1] 0.24 0.76
coin_boot <- sample(x, prob=x_probs_boot, 999, replace=TRUE)
# any draw from here is almost the same as the original process

Jackknife Distribution

Here, we compute all “leave-one-out” estimates. Specifically, for a dataset with \(n\) observations, the jackknife uses \(n-1\) observations other than \(i\) for each unique subsample.

NoteMust Know

Given the sample \(\{1,6,7,22\}\), compute the jackknife estimate of the median. Leave out each observation in turn and take the median of the remaining three.

• Leave out \(1\): median of \(\{6,7,22\}\) is \(7\).
• Leave out \(6\): median of \(\{1,7,22\}\) is \(7\).
• Leave out \(7\): median of \(\{1,6,22\}\) is \(6\).
• Leave out \(22\): median of \(\{1,6,7\}\) is \(6\).

The jackknife distribution of the median is \(\{7,7,6,6\}\), centered at \((7+7+6+6)/4 = 6.5\).

x0_hat <- c(1, 6, 7, 22)
jackknife_meds <- rep(NA, length(x0_hat))
for(i in seq_along(jackknife_meds)){
    jackknife_meds[i] <- median(x0_hat[-i])
}
jackknife_meds
## [1] 7 7 6 6
mean(jackknife_meds)
## [1] 6.5

## MAIN EXAMPLE
# Jackknife Estimates
n <- length(x_hat)
jack_means <- rep(NA, n)
for(i in seq_along(jack_means)){
    x_jack <- x_hat[-i]
    m_i_hat <- mean(x_jack)
    jack_means[i] <- m_i_hat
}
hist(jack_means, breaks=25,
    border=NA, freq=FALSE,
    main=NA, xlab=expression(hat(m)^group('(', -i, ')')))

# The estimates stay close together, since each drops only one of 50 workers.
# The lowest one, about 5.36, drops the highest wage of about 19 dollars.

TipTest Yourself

Repeat the jackknife with the murder arrests in USArrests.

x2_hat <- USArrests[, 'Murder'] # murder arrests per 100,000

# Jackknife Estimates
n2 <- length(x2_hat)
jack_means2 <- rep(NA, n2)
for(i in seq_along(jack_means2)){
    x_jack <- x2_hat[-i]
    m_i_hat <- mean(x_jack)
    jack_means2[i] <- m_i_hat
}
hist(jack_means2, breaks=25,
    border=NA, freq=FALSE,
    main=NA, xlab=expression(hat(m)^group('(', -i, ')')))

Note that both Jackknife and Bootstrap resampling methods provide imperfect estimates, and can give different numbers.

• Jackknife resamples are often less variable than they should be and sample \(n-1\) instead of \(n\).
• Bootstrap resamples have the right \(n\) but often have duplicated data.

5.4 Standard Errors

Having an entire sampling distribution is conceptually clear but unwieldy in practice; we usually compress it to a single number describing its spread.

ImportantKey Definition

The standard error of a statistic is the standard deviation of its sampling distribution.

The standard error is useful for reporting how precise an estimate is with a single number. A small standard error means the statistic would change little across samples. A large standard error means it could swing widely.

For a bootstrap example

• bootstrapped estimate: \(\hat{m}^{(b)}= \frac{1}{n} \sum_{i=1}^{n} \hat{x}_{i}^{(b)}\), for resampled data \(\hat{x}_{i}^{(b)}\)
• mean of the bootstrap distribution: \(\hat{m}^{\text{boot}}= \frac{1}{B} \sum_{b} \hat{m}^{(b)}\).
• standard deviation of the bootstrap distribution: \(\hat{\text{SE}}^{\text{boot}}= \sqrt{ \frac{1}{B} \sum_{b=1}^{B} \left[\hat{m}^{(b)} - \hat{m}^{\text{boot}} \right]^2}\).

Here, the superscript \((b)\) identifies one bootstrap replicate, while \(\text{boot}\) labels summaries computed across the replicates.

NoteMust Know

Using the bootstrap distribution above, with many replicates, we can compute this easily as

sd(boot_means) # standard error estimate for the mean
## [1] 0.4602069
sqrt(mean((x_hat - mean(x_hat))^2)) # standard deviation, which divides by n, unlike sd()
## [1] 3.252514
# The sample sd, about 3.25, is close to the population sd of 3.269.

We can also calculate this explicitly for a small number of replicates. Suppose we have five bootstrap estimates of the sample mean: \(\{3,5,6,7,9\}\). Then

\[\begin{aligned} \hat{m}^{\text{boot}} &= \frac{3 + 5 + 6 + 7 + 9}{5} = 6 \\ \hat{\text{SE}}^{\text{boot}} &= \sqrt{\frac{1}{5}\left[(3-6)^2 + (5-6)^2 + (6-6)^2 + (7-6)^2 + (9-6)^2\right]} \\ &= \sqrt{\frac{1}{5}\left[9 + 1 + 0 + 1 + 9\right]} = \sqrt{\frac{20}{5}} = 2 \end{aligned}\]

TipTest Yourself

For a jackknife example, we have

• jackknifed estimates: \(\hat{m}^{(-i)}=\frac{1}{n-1} \sum_{j \neq i}^{n} \hat{x}_{j}\).
• mean of the jackknife distribution: \(\hat{m}^{\text{jack}}=\frac{1}{n} \sum_{i}^{n} \hat{m}^{(-i)}\).
• standard deviation of the jackknife distribution: \(\hat{\text{SE}}^{\text{jack}} = \sqrt{ \frac{n-1}{n} \sum_{i}^{n} \left[\hat{m}^{(-i)} - \hat{m}^{\text{jack}} \right]^2 }\), approximately \(\sqrt{n}\) times the standard deviation of the \(\hat{m}^{(-i)}\).

The superscript \((-i)\) marks the sample with observation \(i\) omitted.

sd(jack_means)*sqrt(n) # standard error estimate for the mean
## [1] 0.4741274
sd(boot_means) # standard error estimate for the mean
## [1] 0.4602069
# The two estimates nearly agree, at about 0.47 and 0.46.

We modify the sd calculation with a correction, \(\sqrt{n}\), because leave-one-out estimates barely move (only one point changes out of \(n\)). Their spread is therefore smaller than the real sampling spread by a factor of about \(\sqrt{n}\), and the correction puts us back on the right scale.

The pretend population lets us check whether the bootstrap works. The lab data are one sample from it, which we pretend is all we observed. Its bootstrap standard error can then be compared with the true standard error: the standard deviation of the \(500\) sample means wage_means drawn from the population in Sampling Distributions.

## MAIN EXAMPLE
# Bootstrap standard error, from the one lab sample
round(sd(boot_means), 3)
## [1] 0.46

# True standard error, from many samples of the population
round(sd(wage_means), 3)
## [1] 0.436
# The bootstrap estimate lands in the same range as the true standard error.
# It is about 0.46 against about 0.44, so not exactly on it.

The bootstrap used only one sample, so its answer is itself noisy. With skewed data, one sample may include a few very large values or none at all, which moves the bootstrap estimate. To see how much the bootstrap estimate itself changes from one sample to the next, draw a new sample into x_hat and repeat the bootstrap.

TipTest Yourself

Compute the same standard errors for the murder arrests, using the callouts in Bootstrap Distribution and Jackknife Distribution. With no known population, there is no true standard error to compare with.

sd(boot_means2) # bootstrap standard error estimate for the mean
## [1] 0.6218279
sd(jack_means2)*sqrt(n2) # jackknife standard error estimate for the mean
## [1] 0.6285328
sqrt(mean((x2_hat - mean(x2_hat))^2)) # standard deviation, which divides by n, unlike sd()
## [1] 4.311735

Generalization

The above procedure works for many different statistics

med_hat <- quantile(x_hat, prob=0.5)
# The sample median is about 4.84, below the population median of 5.206.

# Bootstrap estimates
boot_medians <- rep(NA, 9999)
for(b in seq_along(boot_medians)){
    x_boot <- sample(x_hat, replace=TRUE) # c.f. jackknife
    med_b_hat <- quantile(x_boot, prob=0.5)
    boot_medians[b] <- med_b_hat
}
SE_med <- sd(boot_medians)

hist(boot_medians, breaks=25,
    border=NA, freq=FALSE,
    main=paste0('SE=', round(SE_med,2)),
    font.main=1,
    xlab=expression(Med[b]))

# The bootstrap medians bunch near a few values in the middle of the sample.

This means we can compute bootstrap SEs for other statistics, too.

TipTest Yourself

Taking the bootstrap distribution of the median, for example, we have

• bootstrapped estimate: \(\tilde{m}_{b}^{\text{boot}}= \text{Med}(\hat{x}_{i}^{(b)})\), for resampled data \(\hat{x}_{i}^{(b)}\)
• mean of the bootstrap distribution: \(\hat{m}^{\text{boot}}= \frac{1}{B} \sum_{b} \tilde{m}_{b}^{\text{boot}}\).
• standard deviation of the bootstrap distribution: \(\hat{\text{SE}}^{\text{boot}}= \sqrt{ \frac{1}{B} \sum_{b=1}^{B} \left[\tilde{m}_{b}^{\text{boot}} - \hat{m}^{\text{boot}} \right]^2}\).

# Standard error estimate for the median
sd(boot_medians)
## [1] 0.3879503
# In this sample it is about 0.39, below the mean's standard error, about 0.46.

TipTest Yourself

Repeat the bootstrap of the median with the murder arrests in USArrests.

x2_hat <- USArrests[, 'Murder'] # murder arrests per 100,000
med2_hat <- quantile(x2_hat, prob=0.5)

# Bootstrap estimates
boot_medians2 <- rep(NA, 9999)
for(b in seq_along(boot_medians2)){
    x_boot <- sample(x2_hat, replace=TRUE) # c.f. jackknife
    med_b_hat <- quantile(x_boot, prob=0.5)
    boot_medians2[b] <- med_b_hat
}
SE_med2 <- sd(boot_medians2)

hist(boot_medians2, breaks=25,
    border=NA, freq=FALSE,
    main=paste0('SE=', round(SE_med2,2)),
    font.main=1,
    xlab=expression(Med[b]))

5.5 Exercises

1. Comment the script you wrote for this chapter, then restart R and check that it runs from a clean session, then check the script with AI as explained in Working with AI. Write three sentences from memory on the main statistical idea of this chapter, and ask the assistant what is wrong, vague, or missing. Finish with your own questions about whatever you found hardest.

2. In your own words, explain the difference between the standard deviation of a sample and the standard error of the sample mean. Why does the standard error decrease as \(n\) grows, while the standard deviation does not necessarily change?

3. Consider a population of four students with ages \(\{18, 20, 22, 24\}\). List all possible simple random samples of size \(n=2\) (without replacement). Compute the sample mean for each. Then compute the mean and standard deviation of those sample means. Verify that the mean of the sample means equals the population mean.

4. Load USArrests in R and extract the UrbanPop column. Write a bootstrap loop with \(B=5000\) resamples to estimate the standard error of the median. That is, for each resample draw n observations with replacement, compute the median, store it, and then compute sd of the stored medians. Compare this bootstrap standard error to the bootstrap standard error of the mean from the same data.

Further Reading

• https://www.khanacademy.org/math/statistics-probability/sampling-distributions-library – Video lessons on sampling distributions, the law of large numbers, and the central limit theorem.
• https://learningstatisticswithr.com/book/estimation.html – Sampling, estimation, and standard errors explained with R code.
• https://www.econometrics-with-r.org/2.2-an-application-to-the-demand-for-economic-journals.html – Sampling and estimation concepts applied to real economic data.

Recall

A sampling distribution shows how a statistic varies across samples. The Law of Large Numbers and Central Limit Theorem describe the sampling distribution of the mean. The Fundamental Theorem of Statistics says the ECDF converges to the true CDF as \(n\) grows, so many other statistics converge too. The jackknife and bootstrap approximate sampling variability from a single dataset, and the standard error summarizes it. The four-student example tied these ideas to concrete numbers: from ages \(\{18, 20, 22, 24\}\) there are six possible samples of size \(n=2\) giving six sample means \(\{19, 20, 21, 21, 22, 23\}\) centered on the population mean \(21\), and the histogram of those six values is the sampling distribution of \(\hat{m}\).

Efron, Bradley. 1979. “Bootstrap Methods: Another Look at the Jackknife.” The Annals of Statistics 7 (1): 1–26. https://doi.org/10.1214/aos/1176344552.

Efron, Bradley, and Robert J. Tibshirani. 1993. An Introduction to the Bootstrap. Vol. 57. Monographs on Statistics and Applied Probability. Chapman & Hall. https://doi.org/10.1007/978-1-4899-4541-9.

Hesterberg, Tim C. 2015. “What Teachers Should Know about the Bootstrap: Resampling in the Undergraduate Statistics Curriculum.” The American Statistician 69 (4): 371–86. https://doi.org/10.1080/00031305.2015.1089789.

Horowitz, Joel L. 2001. “The Bootstrap.” In Handbook of Econometrics, vol. 5. Elsevier. https://doi.org/10.1016/S1573-4412(01)05005-X.

Horowitz, Joel L. 2019. “Bootstrap Methods in Econometrics.” Annual Review of Economics 11 (1): 193–224. https://doi.org/10.1146/annurev-economics-080218-025651.

---

1. A statistic that converges in probability to its target is called consistent.