By the end of this lecture, you will be able to:
We spent much of the first half of the course on linear models, in part due to their ubiquity but also due to our ability to analyze them statistically.
However, these methods for uncertaity quantification, variance reduction, and bias reduction all relied on linear structure. Other predictive/statistical models do not have such convenient analytic properties.
In this module, we will replace mathematical analysis with computational simulation to achieve similar goals.
In this lecture, we will introduce the bootstrap, a computational method for estimating predictive uncertainty in arbitrary predictive models (with almost no assumptions about the statistical or predictive model)!
In the following lectures, we will introduce ensemble methods as generic tools to either provably reduce variance (bagging, random forests) or reduce bias (boosting).
In STAT 306, you produced confidence intervals for the response \(Y\) of a (fixed) test covariate \(x\) under the assumptions of a linear statistical model:
library(Stat406)
library(dplyr)
data(prostate, package = "ElemStatLearn")
prostate <- prostate |> filter(train == TRUE) |> select(-train)
lm_fit <- lm(lpsa ~ . , data = prostate)
test_x <- data.frame(
lcavol = 0, lweight = 3, age = 65,
lbph = 0, svi = 0, lcp = 0, gleason = 7, pgg45 = 0
)
predict(lm_fit, newdata = test_x, interval = "prediction", level = 0.95) fit lwr upr
1 0.8296435 -0.7034454 2.362732These confidence intervals tell us that, under our linearity assumptions, the test response \(Y\) will fall within the interval with 95% probability.
We derive these confidence intervals analytically using our statistical model, properties of our OLS estimator:
\[ \begin{align} \mathbb E[ \hat \beta_\mathrm{OLS}] &= \beta \\ \mathrm{Cov}[ \hat \beta_\mathrm{OLS} \mid \boldsymbol X] &= \sigma^2 (\boldsymbol X^\top \boldsymbol X)^{-1}, \end{align} \]
and by using the fact that \(\hat \beta_\mathrm{OLS} \mid \boldsymbol X\) is Gaussian under our statistical model.
However, consider what happens if we use a predictive model that’s based on a more complicated (or implicit) statistical model, like a decision tree or k-nearest neighbors.
The key idea in this lecture (and module) is to use brute-force computation when closed-form math becomes too difficult.
Imagine that we want to construct a confidence interval around the response \(Y\) for some (fixed) test covariates \(x\) based on our prediction \(\hat Y = \hat f_{\mathcal D}(x)\)
In order to compute a confidence interval for \(Y\) based on the central limit theorem, we need to estimate \(\mathrm{Var}[ \hat Y \mid X=x ] = \mathrm{Var}[ \hat f_{\mathcal D}(x) ]\)
If we had access to \(100\) different training samples, we could train \(100\) different models \(\hat f_{\mathcal D^{(1)}}, \ldots, \hat f_{\mathcal D^{(100)}}\), compute their predictions at \(x\), and estimate this variance empirically.
However, we only have access to one training dataset \(\mathcal D\). We thus need a way to simulate different training datasets \(\mathcal D^{(1)}, \ldots, \mathcal D^{(100)}\) using only the data in \(\mathcal D\).
This idea of (1) simulating new datasets and (2) using those datasets to estimate uncertainty is the core idea behind simulation-based statistics. We use computation to approximate mathematical expectations/variances that we cannot compute analytically.
Imagine that I have \(X_1, \ldots, X_n\) drawn i.i.d. from some distribution \(P(X)\).
The empirical distribution \(\hat P(X)\) is the discrete distribution that places mass \(1/n\) at each observed data point \(X_i\):
\[ \hat P(X=x) = \frac{1}{n} \sum_{i=1}^n \delta_{X_i}(x) \]
where \(\delta_{X_i}(x)\) is a point mass at \(X_i\).
If \(n\) is large enough, then \(\hat P(X) \approx P(X)\):
set.seed(1)
n <- 100
x <- rnorm(n)
hist(x, probability = TRUE, breaks = 10, main = "Empirical Distribution vs True Distribution")
curve(dnorm(x), col = "red", add = TRUE, lwd = 2)
legend("topright", legend = c("True Distribution", "Empirical Distribution"), col = c("red", "black"), lwd = 2)
This idea gives us a recipe to sample many datasets \(\mathcal D^{(1)}, \ldots, \mathcal D^{(B)}\) that are (approximately) drawn from the same distribution as our original dataset \(\mathcal D\), which is often referred to as the bootstrap:
Note that we sample with replacement from the empirical distribution because we want each data point in the new dataset to be an independent draw from \(\hat P(X, Y) \approx P(X, Y)\).
Once we have sampled these new datasets, we can train models \(\hat f_{\mathcal D^{(1)}}, \ldots, \hat f_{\mathcal D^{(B)}}\) on each dataset and use the predictions at \(x\) to estimate the variance of \(\hat Y = \hat f_{\mathcal D}(x)\).
Amazingly, this simple procedure can be used to estimate the variance of the predictions from any model, or just about any statistic that we could ever care about! We don’t need to derive complicated analytic expressions, nor do we even need to make any (significant) assumptions about the statistical model underlying our data!
While we motivated the bootstrap as a way to obtain confidence intervals around responses in learning problems, it’s really a general way to form confidence intervals around any statistical quantity.
Here we will walk through a non-predictive example, since you will look at a predictive example for homework!
Imagine that we have \(X_1, \ldots, X_n\) drawn i.i.d. from \(\mathrm{Exponential}(1/5)\) with pdf \(f(x) = \frac{1}{5}e^{-x/5}\).
Imagine I want to estimate a confidence interval for the sample mean \(\bar X\).
The CLT (if \(n\) is sufficiently large) gives \[ \frac{\sqrt{n}(\bar X - \mathbb E[X])}{{\mathrm{se}}} \approx \mathcal{N}(0, 1), \] where \(\mathrm{se}\) is the standard error of the sample mean, giving us a 95% confidence interval \[ \bar X \pm 2 \hat{\mathrm{se}}/\sqrt{n}. \]
However, what if we wanted to estimate a confidence interval for the sample median instead?
The median is different from the mean for the exponential distribution:
library(tidyverse)
ggplot(data.frame(x = c(0, 12)), aes(x)) +
stat_function(fun = function(x) dexp(x, 1 / 5), color = "orange") +
geom_vline(xintercept = 5, color = "blue") + # mean
geom_vline(xintercept = qexp(.5, 1 / 5), color = "red") + # median
annotate("label",
x = c(2.5, 5.5, 10), y = c(.15, .15, .05),
label = c("median", "mean", "pdf"), parse = TRUE,
color = c("red", "blue", "orange"), size = 6
)
We can use the bootstrap along with the pivotal interval to construct a confidence interval for the sample median!
Assume I have some dataset \(\mathcal D = \{ X_1, \ldots, X_n \}\) of size \(n=500\) drawn i.i.d. from \(\mathrm{Exponential}(1/5)\).
set.seed(406406406)
n <- 500
x <- rexp(n, 1 / 5)First I will draw \(100\) bootstrap samples \(\mathcal D^{(1)}, \ldots, \mathcal D^{(100)}\) by sampling \(n\) points with replacement from my original dataset \(\mathcal D\).
set.seed(406406406)
B <- 100
bootstrap_samples <- map(1:B, ~ sample(x, replace = TRUE))Next, I will compute the sample median on each bootstrap sample:
median_bootstrap <- map_dbl(bootstrap_samples, ~ median(.x))Then I will compute a 95% confidence interval for the sample median using the bootstrap sample medians (using a CI formula that we are about to discuss):
alpha <- 0.05
med <- median(x) # Sample median
# Use the pivotal confidence interval formula for CI
CI <- 2 * med - quantile(median_bootstrap, probs = c(1 - alpha / 2, alpha / 2))The plot below shows the distribution of the bootstrap sample medians:
ggplot(data.frame(median_bootstrap), aes(median_bootstrap)) +
geom_density(color = "orange") +
geom_vline(xintercept = CI, color = "orange", linetype = 2) +
geom_vline(xintercept = med, col = "blue") +
geom_vline(xintercept = qexp(.5, 1 / 5), col = "red") +
annotate("label",
x = c(3.15, 3.5, 3.75), y = c(.5, .5, 1),
color = c("orange", "red", "blue"),
label = c("widehat(F)", "true~median", "widehat(median)"),
parse = TRUE
) +
xlab("x") +
geom_rug(aes(2 * med - median_bootstrap))
This confidence interval is quite good because it contains the true median (the red line) and is fairly tight around the sample median!
Recall that we want to generate a confidence interval for our true response \(Y \mid X=x\) based on our prediction \(\hat Y = \hat f_{\mathcal D}(x)\) and bootstrapped predictions \(\hat Y^{(1)} = \hat f_{\mathcal D^{(1)}}(x), \ldots, \hat Y^{(B)} = \hat f_{\mathcal D^{(B)}}(x)\).
There are two primary techniques used to construct confidence intervals from bootstrap samples.
The percentile interval is given by the empirical quantiles of the bootstrap sample statistics:
\[ [\hat Y^{(b)}_{\alpha/2},\ \hat Y^{(b)}_{1-\alpha/2}] \]
where \(\hat Y^{(b)}_q\) is the \(q\) quantile of \(\hat Y^{(b)}\).
We teach this one in DSCI100 because it’s easy and “looks right”, but it’s not theoretically justified in general. It does have the right coverage (asymptotically) in some specific cases.
When is the Percentile Interval Valid?
If there exists a monotonic function \(m: \mathbb R \to \mathbb R\) so that \(m(\hat Y)\) is Gaussian distributed, i.e. \(m(\hat Y) \sim N(m(\hat Y), c^2)\), then the percentile interval has asymptotically correct coverage.
For most complicated statistics (e.g. the predictions of a decision tree), it’s hard to know whether such a function \(m\) exists.
If we want (asymptotically) correct confidence intervals around the true response, we should use the pivotal interval.
The pivotal interval is given by
\[ [2\hat Y - \hat Y^{(b)}_{1-\alpha/2},\ 2\hat Y - \hat Y^{(b)}_{\alpha/2}] \]
where again \(\hat Y\) is the statistic computed on the original dataset, and \(\hat Y^{(b)}_q\) is the \(q\) quantile of the bootstrap sample statistics.
This interval is justified under fairly general conditions, and thus is the preferred method for constructing bootstrap confidence intervals of the true response \(Y \mid X=x\)
For details on why this interval works, see All of Statistics, Ch. 8.3, or visit me in office hours!
The bootstrap is, at the end of the day, just an approximation. There are two sources of error that could make our variance estimates or confidence intervals incorrect:
Statistical error from using only \(n\) data points to estimate the true data distribution. We don’t have the whole population so we make an error by using a sample. (Note: this part is what always happens with data, and what the science of statistics analyzes.)
Simulation error from using only \(B\) samples to estimate \(\mathrm{Var}[\hat Y]\).
In the homework, you’ll apply the bootstrap to compute uncertainty estimates of machine learning models. In the next lecture, we’ll show how this bootstrap principle can be applied to make predictions more accurate!