UW STAT 535: Statistical Machine Learning (2026 Autumn) Lecture 12: The Bootstrap

Question 1: error of sample median? We start with a simple example: what is the error of sample median? Like sample mean is an estimate of the mean of population, the sample median is an estimate of the median of population. Because it is an estimator, we can define the bias, variance, and mean squared error (MSE) of sample median. But what are these quantities?

Question 2: confidence interval of sample median? Moreover, how can we construct a confidence interval for the population median? We know that given a random sample $X_1,\cdots, X_n\sim F$, a $1-\alpha$ confidence interval of population mean is \[\bar{X}_n \pm z_{1-\alpha/2} \cdot \frac{\widehat{\sigma}_n}{\sqrt{n}},\] where $\bar{X}_n$ and $\widehat{\sigma}_n$ are the sample mean and sample standard deviation. Can we do the same thing (construct a confidence interval) for the median?

In this lecture, we will address these problems for median and many other statistics using the well-known approach: the bootstrap.

The Bootstrap

Here is how we can estimate the error of sample median and construct the corresponding confidence interval. Assume we are given data points $X_1,\cdots,X_n$. Let $\widehat\theta_n$ be our estimator for an unknown parameter of interest $\theta$.

The bootstrap procedure works as follows. First, we sample with replacement from these $n$ points, leading to a set of new observations denoted as $X_1^{*(1)},\cdots,X_n^{*(1)}.$ Again, we repeat the sampling procedure, generating a new sample from the original dataset $X_1,\cdots,X_n$ by sampling with replacement, leading to another new set of observations $X_1^{*(2)},\cdots,X_n^{*(2)}.$ Now we keep repeating the same process of generating new sets of observations; after $B$ rounds, we will obtain \[\begin{align*} X_1^{*(1)},&\cdots,X_n^{*(1)}\\ X_1^{*(2)},&\cdots,X_n^{*(2)}\\ \vdots\qquad&\quad\vdots\quad\quad\vdots\\ X_1^{*(B)},&\cdots,X_n^{*(B)}. \end{align*}\] So in total, we will have $B$ sets of data points. Each set of the data points, say $X_1^{*(1)},\cdots,X_n^{*(1)}$, is called a bootstrap sample. This sampling approach – sample with replacement from the original dataset – is called the empirical bootstrap, invented by Bradley Efron (sometimes this approach is also called Efron’s bootstrap or nonparametric bootstrap)¹. For each bootstrap sample, we can compute the estimator again, leading to $B$ bootstrap estimates $\widehat{\theta}^{*(1)}_n,\cdots, \widehat{\theta}^{*(B)}_n$.

With the bootstrap estimators $\widehat{\theta}^{*(1)}_n,\cdots, \widehat{\theta}^{*(B)}_n$, we can access the uncertainties via various measures.

Example: inference of median. Suppose our parameter of interest $\theta$ is the median of the distribution. Let $\widehat{\theta}_n = {\sf median} \{X_1,\cdots, X_n\}$ be the sample median. Recall that the bootstrap generates $B$ samples from sampling with replacement of the original data: \[\begin{align*} X_1^{*(1)},&\cdots,X_n^{*(1)}\\ X_1^{*(2)},&\cdots,X_n^{*(2)}\\ \vdots\qquad&\quad\vdots\quad\quad\vdots\\ X_1^{*(B)},&\cdots,X_n^{*(B)}. \end{align*}\] We then compute the median of each bootstrap sample: \[\begin{align*} \widehat{\theta}^{*(1)}_n &= {\sf median} \{X^{*(1)}_1,\cdots, X^{*(1)}_n\}\\ \widehat{\theta}^{*(2)}_n &= {\sf median} \{X^{*(2)}_1,\cdots, X^{*(2)}_n\}\\ &\vdots\\ \widehat{\theta}^{*(B)}_n &= {\sf median} \{X^{*(B)}_1,\cdots, X^{*(B)}_n\}. \end{align*}\] A $1-\alpha$ confidence interval of $\theta$ is then $\left[t_{\alpha/2}, t_{1-\alpha/2}\right],$ where $t_{\beta}$ is the $\beta$-percentile of the bootstrap medians $\widehat{\theta}^{*(1)}_n, \cdots, \widehat{\theta}^{*(B)}_n$. As you can see, the bootstrap procedure bypass the needs of computing the standard errors of the original median estimator (sample median).

Other variants of bootstrap confidence intervals

There are many variants of the bootstrap confidence interval. The above percentile method is the most popular approach. Here are three other popular bootstrap confidence intervals.

High level idea of how bootstrap works

Let $X_1,\cdots, X_n\sim F$. Recall that a statistic $S(X_1,\cdots,X_n)$ is a function of random variables so its distribution will depend on the CDF $F$ and the sample size $n$. Thus, the distribution of the estimator $\widehat{\theta}_n$, denoted as $F_{\widehat{\theta}_n}$, will also be determined by the CDF $F$ and sample size $n$. Namely, we may write the CDF of the estimator as \[\begin{equation} F_{\widehat{\theta}_n}(x) = \Psi(x; F,n), \label{eq::BT} \end{equation}\] where $\Psi$ is some complicated function that depends on the CDF of each observation $F$ and the sample size $n$.

When we sample with replacement from $X_1,\cdots, X_n$, what is the distribution we are sampling from? Let $\widehat{F}_n(x) = \frac{1}{n} \sum_{i=1}^n I(X_i\leq x)$ be the EDF of these data points. The EDF is a step function that jumps at each data point. We know that for a discrete random variable, each jump point in its CDF corresponds to the possible value of this random variable and the size of the jump is the probability of selecting that value.

Therefore, if we generate a random variable $Z$ from $\widehat{F}_n$, then $Z$ has the following probability distribution: \[P(Z = X_i) = \frac{1}{n}, \quad \mbox{ for each }i=1,2,\cdots, n.\] If we generated IID $Z_1,\cdots, Z_n\sim \widehat{F}_n$, then the distribution of each $Z_\ell$ is \[P(Z_\ell = X_i) = \frac{1}{n}, \quad \mbox{ for each }i=1,2,\cdots, n, \mbox{ and for all }\ell =1,\cdots, n.\] What is this sample $Z_1,\cdots, Z_n$? This sample is a sample generated by sampling with replacement from $X_1,\cdots, X_n$. Namely, a bootstrap sample is an IID random sample from the EDF $\widehat F_n$ of the original data.

Recall that each set of the bootstrap sample, say $X_1^{*(1)},\cdots,X_n^{*(1)}$, is obtained via sampling with replacement from $X_1,\cdots,X_n$. Thus, each set of the bootstrap sample is an IID sample from $\widehat{F}_n$. Namely, \[\begin{align*} X_1^{*(1)},\cdots,X_n^{*(1)}&\sim \widehat{F}_n\\ X_1^{*(2)},\cdots,X_n^{*(2)}&\sim \widehat{F}_n\\ &\vdots\\ X_1^{*(B)},\cdots,X_n^{*(B)}&\sim \widehat{F}_n. \end{align*}\]

Because a bootstrap estimator, say $\widehat{\theta}^{*(1)}_n$, is the estimator based on $X_1^{*(1)},\cdots,X_n^{*(1)}$. Its CDF, by equation $\eqref{eq::BT}$, is \[F_{\widehat{\theta}^{*(1)}_n}(x) = \Psi(x; \widehat{F}_n,n).\] And because each bootstrap sample is from the distribution $\widehat{F}_n$, we will have \[\Psi(x; \widehat{F}_n,n) = F_{\widehat{\theta}^{*(1)}_n}(x)=F_{\widehat{\theta}^{*(2)}_n}(x)=\cdots=F_{\widehat{\theta}^{*(B)}_n}(x).\] We know that $\widehat{F}_n$ is very similar to $F$ when the sample size is large. Thus, as long as $\Psi$ is smooth (smoothly changing) with respect to $F$, $\Psi(x; \widehat{F}_n,n)$ will also be very similar to $\Psi(x; F,n)$, i.e., \[\widehat{F}_n\approx F \Longrightarrow F_{\widehat{\theta}^{*(\ell)}_n}(x)=\Psi(x; \widehat{F}_n,n)\approx \Psi(x; F,n) = F_{\widehat{\theta}_n}(x).\] This means:

This has many implications. For example, when two CDFs are similar, their variances will be similar as well, i.e., \[{\sf Var}\left(\widehat{\theta}^{*(\ell)}_n|X_1,\cdots,X_n\right)\approx {\sf Var} (\widehat{\theta}_n).\footnote{ The reason why on the left-hand side the variance is conditioned on $X_1,\cdots,X_n$ is because when we compute the bootstrap estimate, the original observations $X_1,\cdots,X_n$ are fixed.}\] Now the bootstrap variance estimate $\widehat{{\sf Var}}_B(\widehat{\theta}_n)$ is just a sample variance of $\widehat{\theta}^{*(\ell)}_n$. When $B$ is large, the sample variance is about the same as the population variance, implying \[\widehat{{\sf Var}}_B(\widehat{\theta}_n) = \frac{1}{B-1}\sum_{\ell=1}^B\left(\widehat{\theta}^{*(\ell)}_n- \bar{\theta}^*_B\right)^2\approx {\sf Var}\left(\widehat{\theta}^{*(\ell)}_n|X_1,\cdots,X_n\right).\] Therefore, \[\widehat{{\sf Var}}_B(\widehat{\theta}_n)\approx{\sf Var}\left(\widehat{\theta}^{*(\ell)}_n|X_1,\cdots,X_n\right)\approx {\sf Var} (\widehat{\theta}_n),\] which explains why the bootstrap variance is a good estimate of the true variance of the estimator.

Generalization to other statistics. The bootstrap can be applied to many other statistics such as sample quantiles, interquartile range, skewness (related to $\mathbb{E}(X^3)$), kurtosis (related to $\mathbb{E}(X^4)$), etc. The theory basically follows from the same idea.

Other variants of bootstrap

From the above argument, we see that as long as we are sampling from a distribution $G_n$ that is converging toward the true underlying distribution $F$, the bootstrap method is expected to work. Therefore, there are two common variants of the bootstrap.

In modern times, we may have many powerful generative models approximating our data generating process. We can also use these generative models for creating bootstrap samples as well.

Failure of the bootstrap

However, the bootstrap may fail for some statistics. One example is the minimum value of a distribution. Here is an illustration of why the bootstrap fails. Let $X_1,\cdots,X_n\sim {\sf Uni}[0,1]$ and $\widehat{\theta}_n = \min \{X_1,\cdots, X_n\}$ be the minimum value of the sample. Then it is known that \[n \cdot \widehat{\theta}_n \overset{D}{\rightarrow} {\sf Exp}(1).\] $\spadesuit:$ Think about why it converges to an exponential distribution.

Thus, $\widehat{\theta}_n$ has a continuous distribution. Assume we generate a bootstrap sample $X^*_1,\cdots,X^*_n$ from the original observations. Now let $\widehat{\theta}_n^*=\min \{X^*_1,\cdots, X^*_n\}$ be the minimum value of a bootstrap sample. Because each $X^*_\ell$ has an equal probability ($\frac{1}{n}$) of selecting each of $X_1,\cdots,X_n$, this implies \[P(X^*_\ell = \widehat{\theta}_n ) = \frac{1}{n}.\] Namely, for each observation in the bootstrap sample, we have a probability of $1/n$ selecting the minimum value of the original sample. Thus, the probability that we do not select $\widehat{\theta}_n$ in the bootstrap sample is \[\begin{align*} P(\mbox{none of }X^*_1,\cdots,X^*_n\mbox{ select }\widehat{\theta}_n) = \left(1-\frac{1}{n}\right)^n \approx e^{-1}. \end{align*}\] This implies that with probability $1-e^{-1}$, one of the observations in the bootstrap sample will select the minimum value of the original sample $\widehat{\theta}_n$. Namely, \[P(\widehat\theta^*_n = \widehat{\theta}_n) = 1-e^{-1}.\] Thus, $\widehat\theta^*_n$ has a huge probability mass at the value $\widehat{\theta}_n$, meaning that the distribution of $\widehat\theta^*_n$ will not be close to an exponential distribution.

Statistical Functionals

To study how the bootstrap works, we first introduce the concept of statistical functionals.

What is a functional? A functional is just a function of a function. Namely, it is a ‘function’ such that the input is another function and the output is a number. Formally speaking, a functional is a mapping $T: \mathcal{F}\mapsto \mathbb{R}$, where $\mathcal{F}$ is a collection of functions. A statistical functional is a mapping $T$ such that you input a distribution (CDF) and it returns a number.

This sounds very complicated but actually, we have encountered numerous statistical functionals. Here are some examples.

Why do we want to use the form of statistical functionals? One answer is: it elegantly describes a population quantity that we may be interested in. Recall that the statistical model about how the data is generated is that we observe a random sample $X_1,\cdots,X_n$ IID from an unknown distribution $F$. Thus, the distribution $F$ is our model for the population. Because the statistical functionals map $F$ into some real numbers, they can be viewed as quantities describing the features of the population. The mean, variance, median, quantiles of $F$ are numbers characterizing the population. Thus, using statistical functionals, we have a more rigorous way to define the concepts of population parameters.

In addition to the above advantage, there is a very powerful feature of statistical functionals – they provide a simple estimator of these population quantities. Recall that the EDF $\widehat{F}_n(x) = \frac{1}{n}\sum_{i=1}^n I(X_i\leq x)$ is a good estimator of $F$. Thus, if we want to estimate a population quantity $\theta=T_{\sf target}(F)$, we can use $T_{\sf target}(\widehat{F}_n)=\widehat{\theta}_n$ as our estimator. Actually, many estimators do follow this form. For instance, in the case of estimating the mean $\mu=T_{\sf mean}(F)$, we often use the sample mean $\bar{X}_n$ as our estimator. However, if you plug in $\widehat{F}_n$ into the statistical functional: \[T_{\sf mean}(\widehat{F}_n) = \int x d\widehat{F}_n(x)= \sum_{i=1}^nX_i\frac{1}{n} = \sum_{i=1}^n\frac{X_i}{n} = \bar{X}_n.\] This implies that the estimator from the statistical functional is the same as the sample mean! Note that in the above calculation, we use the fact that $\widehat{F}_n(x)$ is a distribution whose PMF puts equal probability ($1/n$) at $X_1,\cdots,X_n$. The estimator formed via replacing $F$ by $\widehat{F}_n$ is called a plug-in estimator.

Similarly, we may estimate the variance $\sigma^2 = T_{\sf var} (F)$ via \[T_{\sf var}(\widehat{F}_n) = \int x^2 d\widehat{F}_n(x) - \left( \int x d\widehat{F}_n(x)\right)^2 = \frac{1}{n}\sum_{i=1}^n X_i^2 - \bar{X}_n^2 = \frac{1}{n}\sum_{i=1}^n (X_i-\bar{X}_n)^2 =\frac{n-1}{n}S_n^2.\] This estimator is very similar to the sample variance $S_n^2$ (they are asymptotically the same).

Using how we define the inverse of a CDF of a discrete random variable, we can define the estimator of the median \[T_{\sf med}(\widehat{F}_n) = \widehat{F}_n^{-1}(0.5)\] and other quantiles of a distribution. And it turns out that this estimator is the sample median (and the corresponding sample quantiles)!

Therefore, statistical functionals provide an elegant way to define population quantities as well as an estimator. And the plug-in estimator will be a good estimator if the statistical functional $T(\cdot)$ is ‘smooth’ with respect to the input function so that the smoothness of $T$ with respect to the input will imply $T(\widehat{F}_n)\rightarrow T(F)$².

Bootstrap and Statistical Functionals

We have learned that the (empirical) bootstrap sample is a new random sample from the EDF $\widehat{F}_n$. The bootstrap sample forms another EDF called the bootstrap EDF, denoted as $\widehat{F}_n^*$. Namely, let $X^*_1,\cdots,X^*_n$ be a bootstrap sample. Then the bootstrap EDF is \[\widehat{F}_n^*(x) = \frac{1}{n} \sum_{i=1}^n I(X^*_i\leq x).\]

Here is how the statistical functionals and the bootstrap is connected. In estimating the parameter $\theta = T_{\sf target}(F)$, we often use a plug-in estimate from the EDF $\widehat{\theta}_n = T_{\sf target}(\widehat{F}_n)$ (just think of how we estimate the sample mean). In this case, the bootstrap estimator, the estimator using the bootstrap sample, will be \[\widehat{\theta}_n^* = T_{\sf target}(\widehat{F}^*_n),\] another plug-in estimator but now we are plugging in the bootstrap EDF $\widehat{F}^*_n$.

Consistency of bootstrap variance estimator. How do we use the bootstrap to estimate the variance and construct a confidence interval? We keep generating bootstrap samples from the EDF $\widehat{F}_n$ and obtain several realizations of $\widehat{\theta}_n^*$’s. Namely, we generate \[\widehat{\theta}_n^{*(1)},\cdots, \widehat{\theta}_n^{*(B)}\] and use their sample variance, $\widehat{{\sf Var }}_B(\widehat{\theta}_n^*)$, as an estimator of ${\sf Var}(\widehat{\theta}_n)$. Note that $\widehat{{\sf Var }}_B(\widehat{\theta}_n^*)$ is \[\widehat{{\sf Var }}_B(\widehat{\theta}_n^*) = \frac{1}{B-1}\sum_{\ell=1}^N\left(\widehat{\theta}_n^{*(\ell)}-\bar{\widehat{\theta}}_{n,B}^{*}\right), \quad \bar{\widehat{\theta}}_{n,B}^{*} = \frac{1}{B}\sum_{\ell=1}^B\widehat{\theta}_{n}^{*(\ell)}.\]

When $B$ is large, the sample variance of the bootstrap estimators \[\begin{equation} \widehat{{\sf Var }}_B(\widehat{\theta}_n^*) \approx{\sf Var }(\widehat{\theta}_n^*|\widehat{F}_n). \label{eq::BT1} \end{equation}\] Note that $\cdot|\widehat{F}_n$ means conditioned on $\widehat{F}_n$ being fixed. The reason why here it converges to this conditioned variance is because when we generate bootstrap samples, the original EDF $\widehat{F}_n$ is fixed (and we are generating from it). Thus, the variance is conditioned on $\widehat{F}_n$ being fixed.

To argue that the bootstrap variance $\widehat{{\sf Var }}_B(\widehat{\theta}_n^*)$ is a good estimate of the original variance, we need to argue \[\begin{equation*} \widehat{{\sf Var }}_B(\widehat{\theta}_n^*) \approx{\sf Var }(\widehat{\theta}_n^*|\widehat{F}_n)\approx {\sf Var}(\widehat{\theta}_n). \end{equation*}\] However, because of equation $\ref{eq::BT1}$ and we can select $B$ as large as we wish, so what really matters is \[\begin{equation*} {\sf Var }(\widehat{\theta}_n^*|\widehat{F}_n)\approx {\sf Var}(\widehat{\theta}_n). %\label{eq::BT2} \end{equation*}\] Or more formally, \[\begin{equation} \frac{{\sf Var }(\widehat{\theta}_n^*|\widehat{F}_n)}{{\sf Var }(\widehat{\theta}_n)} \approx 1 \label{eq::BT3} \end{equation}\] (people generally use the ratio expression because both quantities often converge to $0$ when the sample size $n\rightarrow \infty$).

Because $\widehat{\theta}_n = T_{\sf target}(\widehat{F}_n)$ is a statistic (a function of our random sample $X_1,\cdots,X_n$), its distribution is completely determined by the distribution $X_1,\cdots,X_n$ are sampling from, which is $F$, and the sample size $n$. This implies that the variance of $\widehat{\theta}_n$ is determined by $F$ and $n$ as well. Therefore, we can write \[{\sf Var}(\widehat{\theta}_n) = {\sf Var}(T_{\sf target}(\widehat{F}_n)) = \mathbb{V}_{n,\sf target}(F).\] And it turns out that we often have \[\mathbb{V}_{n,\sf target}(F) \approx \frac{1}{n}\mathbb{V}_{1,\sf target}(F)\equiv \frac{1}{n}\mathbb{V}_{\sf target}(F).\] Note that here $\mathbb{V}_{n,\sf target}(\cdot), \mathbb{V}_{\sf targe}(\cdot)$ are both again statistical functionals!

Because the bootstrap estimator $\widehat{\theta}^*_n = T_{\sf target}(\widehat{F}^*_n)$, its conditional variance will be \[{\sf Var}(\widehat{\theta}^*_n|\widehat{F}_n) = {\sf Var}(T_{\sf target}(\widehat{F}^*_n)|\widehat{F}_n) = \mathbb{V}_{n,\sf target}(\widehat{F}_n) \approx \frac{1}{n}\mathbb{V}_{\sf target}(\widehat{F}_n).\] Thus, as long as \[\begin{equation} \mathbb{V}_{\sf target}(\widehat{F}_n)\approx\mathbb{V}_{\sf target}(F), \label{eq::BT4} \end{equation}\] equation $\ref{eq::BT3}$ holds. Namely, the bootstrap variance estimate will be a good estimator of the variance of the true estimator³.

Validity of bootstrap confidence interval. How about the validity of the bootstrap confidence interval? Here is a derivation showing that the consistency of bootstrap variance estimator implies the validity of bootstrap confidence interval.

For the bootstrap confidence interval, a simple way is first show that \[\begin{equation} \sqrt{n}(\widehat{\theta}_n - \theta) = \sqrt{n}\left(T_{\sf target}(\widehat{F}_n) - T_{\sf target}(F)\right) \approx N(0, \mathbb{V}_{\sf target}(F)) \label{eq::BT5} \end{equation}\] which implies \[\sqrt{n}(\widehat{\theta}^*_n - \widehat{\theta}_n) = \sqrt{n}\left(T_{\sf target}(\widehat{F^*}_n) - T_{\sf target}(\widehat{F}_n)\right) \approx N(0, \mathbb{V}_{\sf target}(\widehat{F}_n)).\] Thus, as long as the bootstrap variance converges, we also have the convergence of the entire distribution, implying the validity of a bootstrap confidence interval. Note that to formally prove this, we need to show the convergence in terms of CDF of the difference. In more details, let $Z_n = \sqrt{n}(\widehat{\theta}_n - \theta)$ and $Z_n^* = \sqrt{n}(\widehat{\theta}^*_n - \widehat{\theta}_n)$. We need to prove \[\sup_t\left|P(Z_n^*\leq t|\widehat{F}_n) - P(Z_n\leq t)\right|\overset{P}{\rightarrow } 0.\] Later we will show examples about this using the sample mean as a starting point.

Example: mean. We now consider a simple example: the mean of a distribution $T_{\sf target} = T_{\sf mean}$. The mean of a distribution has the form \[\mu = T_{\sf mean}(F) = \int x dF(x).\] The plug-in estimator is \[\widehat{\mu}_n = T_{\sf mean}(\widehat{F}_n) = \int x d\widehat{F}_n(x) = \bar{X}_n\] and the bootstrap estimator is \[\widehat{\mu}^*_n = T_{\sf mean}(\widehat{F}^*_n) = \int x d\widehat{F}^*_n(x) = \bar{X}^*_n.\] It is clear from the Central Limit Theorem that \[\sqrt{n}(\widehat{\mu}_n-\mu) \approx N(0, n{\sf Var}( T_{\sf mean}(\widehat{F}_n)))\] so equation $\ref{eq::BT5}$ holds and \[\sqrt{n}(\widehat{\mu}^*_n-\widehat{\mu}_n) \approx N(0, n{\sf Var}( T_{\sf mean}(\widehat{F}^*_n)|\widehat{F}_n)).\]

In this case, we know that \[{\sf Var}( T_{\sf mean}(\widehat{F}_n)) = {\sf Var}(\bar{X}_n) = \frac{1}{n}{\sf Var}(X_i)\Longrightarrow \mathbb{V}_{\sf mean}(F) = {\sf Var}(X_i) = \mathbb{E}(X_i^2) - \mathbb{E}^2(X_i) = \int x^2 dF(x) - \left(\int xdF(x) \right)^2.\] Therefore, the bootstrap variance is \[{\sf Var}( T_{\sf mean}(\widehat{F}^*_n)|\widehat{F}_n) = \frac{1}{n} \mathbb{V}_{\sf mean}(\widehat{F}_n) = \int x^2 d\widehat{F}_n(x) - \left(\int xd\widehat{F}_n(x) \right)^2.\] Because of the Law of Large Numbers, \[\begin{align*} \int x^2 d\widehat{F}_n(x) = \frac{1}{n}\sum_{i=1}^nX_i^2&\overset{P}{\rightarrow} \mathbb{E}(X_i^2) = \int x^2 dP(x)\\ \int x d\widehat{F}_n(x) = \frac{1}{n}\sum_{i=1}^nX_i&\overset{P}{\rightarrow} \mathbb{E}(X_i) = \int x dP(x). \end{align*}\] Thus, ⁴ \[\mathbb{V}_{\sf mean}(\widehat{F}_n) \overset{P}{\rightarrow} \mathbb{V}_{\sf mean}(F),\] which shows that equation $\ref{eq::BT4}$ holds and so does equation $\ref{eq::BT3}$. Thus, the bootstrap variance estimator converges to the true variance estimator and we conclude that \[\frac{{\sf Var}( T_{\sf mean}(\widehat{F}^*_n)|\widehat{F}_n)}{{\sf Var}( T_{\sf mean}(\widehat{F}_n))} \overset{P}{\rightarrow} 1.\] As a result, the bootstrap variance estimator is consistent and the bootstrap confidence interval is also valid.

Consistency of Bootstrap

The analysis in the previous section is a high-level sketch of the consistency of bootstrap. In this section, we provide a formal way to derive the consistency of bootstrap.

Let $Z_n = \sqrt{n}(\widehat{\theta}_n - \theta)$ and $Z_n^* = \sqrt{n}(\widehat{\theta}^*_n - \widehat{\theta}_n)$. To formally prove the validity bootstrap, we need to prove that \[\begin{equation} \sup_t\left|P(Z_n^*\leq t|\widehat{F}_n) - P(Z_n\leq t)\right|\overset{P}{\rightarrow } 0. \label{eq::uniform} \end{equation}\] The above bound is also known as the Kolomogrov distance between two random variables.

Although this seems to be hard to prove, there are two popular approaches to derive equation $\ref{eq::uniform}$. The first approach is to show that $Z_n$ has an asymptotic linear form and then apply Lindeberg-Feller central limit theorem (triangular arrays) to $Z_n^*$ since $Z_n^*$ is sampled from a ‘random distribution function’ $\widehat F_n$. The second approach is via the Berry-Esseen bound of the sample mean, which is our preferred route.

Lindeberg-Feller’s central limit theorem

In the plug-in estimate of statistical functionals, we see that many estimators can be written as a sample mean problem. Therefore, here we consider a simple scenario that we observe univariate $X_1,\cdots, X_n$ and we are interested in estimating the population mean, i.e., $\theta = \mathbb{E}(X_1)$, using the sample mean $\widehat\theta_n =\bar X_n = \frac{1}{n}\sum_{i=1}^n$.

The conventional central limit theorem (CLT) shows that \[\sqrt{n}(\widehat\theta_n - \theta) \overset{d}{\rightarrow}N(0,\sigma^2),\] where $\sigma^2 = {\sf Var}(X_1)$. This result is NOT enough for bootstrap consistency because it assumes that observations $X_1,\cdots, X_n$ are sampled from a fixed CDF $F$, not a distribution function that can change with respect to the sample size $n$. In the bootstrap case, the bootstrap sample $X^*_1,\cdots, X^*_n$ are IID from the EDF $\widehat F_n$ given $X_1,\cdots, X_n$. Thus, the ‘population’ of the bootstrap sample changes with respect to $n$, so conventional CLT is not applicable.

To resolve this issue, we use the Lindeberg-Feller’s CLT (also known as triangular array CLT).

Theorem 1 (Lindeberg-Feller). For each $n=1,2,3,\cdots$, let $W_n = (W_{n,1},\cdots, W_{n,k_n})$ be a vector of independent elements with finite variance, i.e., $W_{n,1},\cdots, W_{n,k_n}$ are independent from each other. Assume that

Uniform integrability. $\sum_{i=1}^{k_n}\mathbb{E}[W^2_{n,i}I(|W_{n,i}|>\epsilon)]\rightarrow0$ for every $\epsilon>0$.
Finite variance. $\sum_{i=1}^{k_n} {\sf Var}(W_{n,i}) \rightarrow \sigma^2$.

Then \[\sum_{i=1}^{k_n} (W_{n,i} - \mathbb{E}(W_{n,i}))\overset{d}{\rightarrow}N(0, \sigma^2).\]

We consider the scenario where the original data $X_1,\cdots, X_n$ are fixed so that the bootstrap sample $X^*_1,\cdots, X^*_n$ are IID from $\widehat F_n$.

To use Theorem 1 in the bootstrap setting, each $W_{n, i} = \frac{1}{\sqrt{n}}X^*_i$ is sampled from $\widehat F_n$, so $\mathbb{E}_{\widehat F_n}[W_{n, i}] = \frac{1}{\sqrt{n}} \widehat\theta_n$. Note that the expectation $\mathbb{E}_{\widehat F_n}(\cdot)$ is with respect to the distribution $\widehat F_n$. Also, $k_n = n$. Under this setting, $\sum_{i=1}^{k_n} W_{n,i} = \frac{1}{\sqrt{n}}\sum_{i=1}^n X^*_{i} = \sqrt{n}\widehat\theta_n^*$ and \[\sum_{i=1}^{k_n} (W_{n,i} - \mathbb{E}_{\widehat F_n}(W_{n,i})) = \sqrt{n}(\widehat\theta_n^* - \widehat\theta_n).\] Thus, the conclusion of Theorem 1 is applicable to the setting of the bootstrap.

Now we investigate the two conditions in Theorem 1. The first uniform integrability condition \[\sum_{i=1}^{k_n}\mathbb{E}_{\widehat F_n}[W^2_{n,i}I(|W_{n,i}|>\epsilon)]\rightarrow0\] becomes \[\frac{1}{n}\sum_{i=1}^n X^{2}_{i}I(|X_{i}|>\sqrt{n}\epsilon) \rightarrow0,\] When the true distribution $F$ has a finite second moment, i.e., $\mathbb{E}(X_i^2)<\infty$, strong Law of Large Numbers implies $\frac{1}{n}\sum_{i=1}^n X^{2}_{i}\overset{a.s.}{\rightarrow} \mathbb{E}(X_i^2)<\infty$, so $\frac{1}{n}\sum_{i=1}^n X^{2}_{i}I(|X_{i}|>\sqrt{n}\epsilon) \overset{a.s.}{\rightarrow}0$.

The finite variance condition becomes \[\sum_{i=1}^{k_n} {\sf Var}(W_{n,i}) = \frac{1}{n}\sum_{i=1}^n {\sf Var}_{\widehat F_n}(X^*_i) = \widehat\sigma^*_n \overset{P}{\rightarrow} \sigma^2,\] where $\widehat\sigma^*_n = \frac{1}{n}\sum_{i=1}^n (X_i-\bar X_n)^2$ and $\sigma^2 = {\sf Var}(X_1)$.

As a result, we conclude that \[\sqrt{n}(\widehat\theta_n^* - \widehat\theta_n) \overset{d}{\rightarrow} N(0,\sigma^2);\] namely, it converges to the same limit as the original estimator $\sqrt{n}(\widehat\theta_n - \theta)\overset{d}{\rightarrow} N(0,\sigma^2)$. Thus, the equation $\ref{eq::uniform}$ holds for bootstrapping $\widehat\theta_n$.

A more general form of this approach can be found in Theorem 23.4 of [Van der Vaart (2000)].

Berry-Esseen bound

Consider again the simple scenario that we observe univariate $X_1,\cdots, X_n$ and we are interested in estimating the population mean, i.e., $\theta = \mathbb{E}(X_1)$.

Theorem 2 (Berry-Esseen bound). Assume that $\mathbb{E}(|X_1|^3)<\infty$. Let $Z\sim N(0,1)$ and $\theta = \mathbb{E}(X_1)$ and $\sigma^2 = {\sf Var}(X_1)$. Then for any $n$, we have \[\sup_t\left|P\left(\sqrt{n}\left(\frac{\bar{X}_n - \theta}{\sigma}\right)<t\right) - P(Z<t)\right|\leq C\frac{\mathbb{E}|X_1|^3}{\sigma^3\sqrt{n}},\] for a constant $C\geq \frac{\sqrt{10}+3}{6\sqrt{2\pi}}$.

The Berry-Esseen bound quantifies how fast the limiting distribution converges to a Gaussian and the result is uniform across different quantiles.

Finite Sample bound and the bootstrap. It is important to note that the Berry-Esseen bound is a finite sample bound, meaning that its result holds for any $n$ (some finite sample bound holds when $n$ is larger than some constant). So it is a much stronger result than the conventional central limit theorem. The finite sample bound is important in deriving the validity of the bootstrap (see the proof below).

The Berry-Esseen bound can be used to derive bounds like equation $\ref{eq::uniform}$. Now consider very simple scenario that we are interested in estimating the population mean $\theta = \mathbb{E}(X_1)$ and we use the sample mean as the estimator $\widehat\theta_n$.

Theorem 3. Suppose that we are considering the sample mean problem, i.e., $\theta = \mathbb{E}(X_1)$ and $\widehat\theta_n = \bar X_n$ is the original sample mean estimator and $\widehat\theta_n^* = \bar X_n^*$ is the sample mean of the bootstrap sample. Assume that $\mathbb{E}(|X_1|^3)<\infty$. Let \[Z_n = \sqrt{n}(\widehat\theta_n - \theta),\qquad Z_n^* = \sqrt{n} (\widehat\theta_n^* - \widehat\theta_n).\] Then \[\sup_t\left|P(Z_n^*\leq t|\widehat{F}_n) - P(Z_n\leq t)\right|= O_P\left(\frac{1}{\sqrt{n}}\right).\]

Let $\Psi_{\sigma}(t)$ be the CDF of $N(0,\sigma^2)$ and $\widehat\sigma^2 = \frac{1}{n}\sum_{i=1}^n (X_i-\bar X_n)^2$. We bound the difference using \[\begin{align*} \sup_t\left|P(Z_n^*\leq t|\widehat{F}_n) - P(Z_n\leq t)\right| & \leq \sup_t\left|P(Z_n^*\leq t|\widehat{F}_n) - \Psi_{\widehat\sigma}(t)\right| + \sup_t\left|\Psi_{\widehat\sigma}(t)) - \Psi_{ \sigma}(t)\right|+ \sup_t\left|P(Z_n\leq t) - \Psi_{ \sigma}(t)\right|. \end{align*}\]

The Berry Esseen theorem implies that \[\sup_t |P(Z_n\leq t) - \Psi_{\sigma}(t)| = O_P\left(\frac{1}{\sqrt{n}}\right)\] so the third quantity is bounded. Similarly, we can apply the Berry-Esseen bound to the first quantity by replacing $\mathbb{E}(\cdot)$ with the empirical version of it (sample average operation), which implies \[\sup_t\left|P(Z_n^*\leq t|\widehat{F}_n) - P(Z_n\leq t)\right| \leq C\frac{\frac{1}{n}\sum_{i=1}^nX_i^3}{\sigma^3\sqrt{n}}.\] Note that we can apply the Berry-Esseen theory to the bootstrap because this theory holds in finite sample! In the bootstrap world, the EDF is the population distribution generating our data, and that is why we replace the expectation $\mathbb{E}$ by the empirical version of it.

By strong Law of Large Numbers, the probability that the right hand side is less than $2C\frac{\mathbb{E}|X_1|^3}{\sigma^3\sqrt{n}}$ is $1$. Thus, we conclude that \[\sup_t\left|P(Z_n^*\leq t|\widehat{F}_n) - P(Z_n\leq t)\right| = O_P\left(\frac{1}{\sqrt{n}}\right).\]

For the second term, $\sup_t\left|\Psi_{\widehat\sigma}(t)) - \Psi_{ \sigma}(t)\right|$, because $|\widehat\sigma - \sigma| = O_P\left(\frac{1}{\sqrt{n}}\right)$ so differentiating the CDF with respect to $\sigma$ and take a uniform bound leads to \[\sup_t\left|\Psi_{\widehat\sigma}(t)) - \Psi_{ \sigma}(t)\right| = O_P\left(\frac{1}{\sqrt{n}}\right),\] which completes the proof.

The Lindeberg-Feller central limit theorem approach requires a slightly less condition than the Berry-Esseen bound (we do not need third-moment but just need a bounded second moment). However, the Lindeberg-Feller approach will not give us a convergence rate while the Berry-Esseen approach gives us a convergence rate.

Delta Method

In this section, we will talk about a very useful technique in handling the convergence–the delta method. We start with an example of proving consistency theorem of some bootstrap estimates.

Example: inverse of mean. Assume we are interested in the inverse of the population mean. Namely, the statistical functional we will be using is \[T_{\sf inv}(F) = \frac{1}{\int x dF(x)} = \lambda.\] This statistical functional was implicitly used when we compute the MLE of the rate parameter of an exponential distribution. The plug-in estimator (as well as the MLE of estimating an exponential model) is \[\widehat{\lambda}_n = T_{\sf inv}(\widehat{F}_n) = \frac{1}{\int x d\widehat{F}_n(x)} = \frac{1}{\bar{X}_n}.\] The bootstrap estimator is \[\widehat{\lambda}^*_n = T_{\sf inv}(\widehat{F}^*_n) = \frac{1}{\int x d\widehat{F}^*_n(x)} = \frac{1}{\bar{X}^*_n}.\]

In the lab session, we have shown that this estimator follows asymptotically a normal distribution. But how do we show this, and how do we compute the variance of the estimator $\widehat{\lambda}_n$? Here is how the delta method will help us.

Now using equation $\ref{eq::delta02}$ and identifying $Y_n$ as $\bar{X}_n$ and $f(x)$ as $\frac{1}{x}$, we obtain \[%\sqrt{n}(\hat{\lambda}_n -\lambda) = \sqrt{n}\left(\frac{1}{\bar{X}_n} - \frac{1}{\E(X_i)}\right) \approx -\frac{1}{\E^2(X_i)} \sqrt{n} \left(\bar{X}_n-\frac{1}{\E^2(X_i)}\right) \sqrt{n}(\widehat{\lambda}_n -\lambda) = \sqrt{n}\left(\frac{1}{\bar{X}_n} - \frac{1}{\mathbb{E}(X_i)}\right) \approx -\frac{1}{\mathbb{E}^2(X_i)} \sqrt{n} \left(\bar{X}_n-{\mathbb{E}(X_i)}\right) \approx N\left(0, \underbrace{\frac{1}{\mathbb{E}^4(X_i)}{\sf Var}(X_i)}_{=\mathbb{V}_{\sf inv}(F)}\right).\] Using the fact that $\mathbb{E}(X_i) = \int x dF(x)$ and ${\sf Var}(X_i) = \int x^2 dF(x) - \left(\int xdF(x) \right)^2$, we obtain \[\sqrt{n}(\widehat{\lambda}_n -\lambda) \approx N(0, \mathbb{V}_{\sf inv}(F)),\] where \[\mathbb{V}_{\sf inv}(F) = \frac{ \int x^2 dF(x) - \left(\int xdF(x) \right)^2}{\left(\int x dF(x)\right)^4}.\] So equation $\ref{eq::BT5}$ holds and \[\sqrt{n}(\widehat{\lambda}^*_n -\widehat{\lambda}_n) \approx N(0, \mathbb{V}_{\sf inv}(\widehat{F}_n)),\] where \[\mathbb{V}_{\sf inv}(\widehat{F}_n) = \frac{ \int x^2 d\widehat{F}_n(x) - \left(\int xd\widehat{F}_n(x) \right)^2}{\left(\int x d\widehat{F}_n(x)\right)^4}\] is the corresponding bootstrap variance component.

As long as $\int x dF(x) \neq 0$, each component in $\mathbb{V}_{\sf inv}(\widehat{F}_n)$ is a natural estimator of the corresponding component in $\mathbb{V}_{\sf inv}(F)$. Therefore, we conclude \[\mathbb{V}_{\sf inv}(\widehat{F}_n)\overset{P}{\rightarrow}\mathbb{V}_{\sf inv}(F),\] which shows that equation $\ref{eq::BT4}$ holds, implying that the bootstrap variance estimator is consistent: \[\frac{{\sf Var}( T_{\sf inv}(\widehat{F}^*_n)|\widehat{F}_n)}{{\sf Var}( T_{\sf inv}(\widehat{F}_n))} \overset{P}{\rightarrow} 1\] and moreover, the bootstrap confidence interval is also valid.

Influence Function

Linear Functional

In the above derivations, we see many examples of statistical functionals that are of the form \[T_\omega(F) = \int \omega(x) dF(x),\] where $g$ is a function. As we have mentioned, this type of statistical functionals are called linear functionals.

Linear functionals has a feature that the estimators \[\begin{align*} T_\omega(\widehat{F}_n) &= \int \omega(x) d\widehat{F}_n(x)= \frac{1}{n}\sum_{i=1}^n g(X_i),\\ T_\omega(\widehat{F}^*_n) &= \int \omega(x) d\widehat{F}^*_n(x)= \frac{1}{n}\sum_{i=1}^n g(X^*_i). \end{align*}\]

Moreover, a powerful feature of the linear functional is that for another CDF $G$, we always have \[\begin{align*} T_\omega(G) - T_\omega(F) &= \int \omega(x)dG(x)-T_\omega(F) \\ & = \int \omega(x)dG(x) - \int T_\omega(F) dG(x)\\ & = \int L_F(x) dG(x), \end{align*}\] where \[\begin{equation} L_F(x) = \omega(x) - T_\omega(F) \label{eq::LF1} \end{equation}\] is called the influence function of the functional $T_\omega$.

Theorem 4. Suppose that $T_\omega$ is a linear functional with an influence function $L_F(x)$ define in equation $\ref{eq::LF1}$ and $\int \omega^2(x)dF(x)<\infty$ . Then \[\sqrt{n}\left(T_\omega(\widehat{F}_n) - T_\omega(F)\right) \overset{D}{\rightarrow} N\left(0, \mathbb{V}_\omega(F) = \int L^2_F(x)dF(x)\right)\] and a consistent estimator of $\mathbb{V}_\omega(F)$ is $\mathbb{V}_\omega(\widehat F_n) = \frac{1}{n}\sum_{i=1}^n L^2_F(X_i)$.

As a result, the bootstrap always works for the linear functional whenever $T_{\omega^2}(F)<\infty$.

It is easy to see that \[T_\omega(\widehat{F}_n) - T_\omega(F) = \int L_F(x)d\widehat{F}_n(x) = \frac{1}{n}\sum_{i=1}^n L_F(X_i).\] Moreover, \[\mathbb{E}(L_F(X_i)) = \int L_F(x)dF(x) = \int\left(\omega(x) - T_\omega(F)\right) dF(x) = T_\omega(F)-T_\omega(F) = 0.\] Thus, by central limit theorem, \[\sqrt{n}\left(T_\omega(\widehat{F}_n) - T_\omega(F)\right) \overset{D}{\rightarrow} N\left(0, \mathbb{V}_\omega(F) = \int L^2_F(x)dF(x)\right)\] Note that for a linear functional $T_\omega$, equation $\ref{eq::BT5}$ always holds with \[\begin{equation} \mathbb{V}_\omega(F) = \int L^2_F(x)dF(x) = \int \left(\omega^2(x)-2\omega(x)T_\omega(F)-T^2_\omega(F)\right)dF(x) = \int \omega^2(x)dF(x) = T_{\omega^2}(F). \label{eq::VIF1} \end{equation}\]

Moreover, \[\begin{equation} \begin{aligned} \mathbb{V}_\omega(\widehat{F}_n) &= \int L^2_{\widehat{F}_n}(x)d\widehat{F}_n(x) = \int \left(\omega^2(x)-2\omega(x)T_\omega(\widehat{F}_n) + T^2_\omega(\widehat{F}_n)\right)d\widehat{F}_n(x)\\ & = \int \omega^2(x)d\widehat{F}_n(x) - T^2_\omega(\widehat{F}_n). \end{aligned} \label{eq::VIF2} \end{equation}\] By Law of Large Numbers (and continuous mapping theorem), \[T^2_\omega(\widehat{F}_n)\overset{P}{\rightarrow} T^2_\omega(F)\] if $\mathbb{E}(|\omega(X_i)|) = T_{|\omega|}<\infty$. And \[\int \omega^2(x)d\widehat{F}_n(x)=T_{\omega^2}(\widehat{F}_n) \overset{P}{\rightarrow} T_{\omega^2}(F) %= \V_\omega(F)\] if $\mathbb{E}(\omega(X_i)^2) = T_{\omega^2}(F)<\infty$. Therefore, we conclude that when $T_{\omega^2}(F)<\infty$, \[\mathbb{V}_\omega(\widehat{F}_n)= \int \omega^2(x)d\widehat{F}_n(x) - T^2_\omega(\widehat{F}_n) \overset{P}{\rightarrow} \mathbb{V}_\omega(F) %+0^2 = \mathbb{V}_\omega(F),\] implying that equation $\ref{eq::BT4}$ holds and so does equation $\ref{eq::BT3}$.

Non-linear Functional

Although the linear functional has so many beautiful properties, many statistical functionals are not linear. For instance, the median \[T_{\sf med}(F) = F^{-1}(0.5)\] is not a linear functional. Therefore, our results of linear functional cannot be directly applied to analyze the median.

Then how can we analyze the properties of non-linear statistical functionals? One way to proceed is to generalize the notion of influence function. And here is the formal definition of the influence function.

Let $\delta_x$ be a point mass at location $x$. The influence function of a (general) statistical function $T_{\sf target}$ is \[\begin{equation} L_{F}(x) = \lim_{\epsilon\rightarrow 0}\frac{T_{\sf target}((1-\epsilon)F + \epsilon\delta_x)-T_{\sf target}(F)}{\epsilon}. \label{eq::IF1} \end{equation}\] Some of you may find equation $\ref{eq::IF1}$ very familiar; it seems to be taking a derivative. And yes – it is a derivative of a functional with respect to a function. This type of derivative is called Gâteaux derivative⁵, a type of derivative of functionals. You can check that applying equation $\ref{eq::IF1}$ to a linear functional leads to an influence function as we defined previously.

A powerful feature of this generalized version of influence function is that when the statistical functional $T_{\sf target}$ is ‘smooth⁶’, equation $\ref{eq::VIF1}$ and $\ref{eq::VIF2}$ hold in the sense that \[\begin{equation} \begin{aligned} \mathbb{V}_{\sf target}(F) = \int L^2_F(x)dF(x),\quad \mathbb{V}_{\sf target}(\widehat{F}_n) &= \int L^2_{\widehat{F}_n}(x)d\widehat{F}_n(x) \end{aligned} \end{equation}\] and, moreover, equation $\ref{eq::BT5}$ holds. Note that $L_{\widehat{F}_n}(x)$ is defined via replacing $F$ by $\widehat{F}_n$ in equation $\ref{eq::IF1}$. That is, when the statistical functional $T_{\sf target}$ is smooth, we only need to verify \[\begin{equation} \int L^2_{\widehat{F}_n}(x)d\widehat{F}_n(x) \approx \int L^2_F(x)dF(x) \label{eq::FD} \end{equation}\] to argue the validity of bootstrap consistency.

Example: median. Why does the median follow a normal distribution? Here we will show this using the influence function. The influence function of the functional $T_{\sf med}$ is \[L_{F}(x) =\frac{1}{2p(F^{-1}(0.5))},\] where $p$ is the PDF of $F$ (you can verify it). Thus, equation $\ref{eq::BT5}$ implies \[\sqrt{n}\left(\underbrace{T_{\sf med}(\widehat{F}_n)}_{\mbox{sample median}}-\underbrace{T_{\sf med}(F)}_{\mbox{population median}}\right) \approx N\left(0, \frac{1}{4p^2(F^{-1}(0.5))}\right).\] Note that $F^{-1}(0.5) = T_{\sf med}(F)$ is the median of $F$. So this shows not only the asymptotic normality of the sample median but also its limiting variance, which is inversely related to the PDF at the median.

The influence function is also related to the robustness of an estimator⁷ and plays a key role in the semi-parametric statistics⁸. You would encounter it several times if you want to pursue a career in statistics.

Functional Delta Method

The Berry-Esseen theory shows that we can use the bootstrap to a sample mean problem and the delta method further implies that the bootstrap is applicable to any statistical functional that can be written as a smoothed function of a sample mean.

However, in many scenarios such as the sample median, the regular delta method does not work but the analysis using influence function shows that the bootstrap is still applicable. This motivates us to generalize the regular delta method. It turns out that there is a much wider classes of statistical functionals that the bootstrap method works. Here we will introduce a new technique called functional delta method that generalizes the regular delta method to a much wider class of statistical functionals.

Before we formally explain the functional delta method, we first introduce the concept of Hadamard differentiation. Let $\mathbb{D}, \mathbb{F}$ be two normed spaces (you can think of function spaces) and let $\phi:\mathbb{D}\mapsto \mathbb{F}$ be a mapping. $\phi$ is said to be Hadamard differentiable at $\omega\in\mathbb{D}$ with a differentiation $\dot{\phi}_\omega$ if for any sequence $\eta_t\rightarrow\eta$ when $t\rightarrow 0$, \[\lim_{t\rightarrow0}\left\|\frac{\phi(\omega+t\cdot\eta_t) - \phi(\omega)}{t} - \dot{\phi}_\omega(\eta)\right\|_{\mathbb{F}}=0,\] where $\|\cdot\|_{\mathbb{F}}$ is the norm of the space $\mathbb{F}$.

The followings are two informal statements about the functional delta method, which are simplified from

Theorem 5 (Functional delta method). Let $X_1,\cdots,X_n\sim F_0$ and $\widehat{F}_n$ be the EDF and $\theta(F_0)$ be the parameter of interest. If $\theta$ is Hadamard differentiable at $F_0$ with respect to the $L_\infty$ norm, then \[\sqrt{n}(\theta(\widehat F_n) - \theta(F_0)) \overset{D}{\rightarrow} \dot{\theta}_{F_0}(\mathbb{B}),\] where $\mathbb{B}$ is a Gaussian process defined over $\mathbb{R}^d$ with ${\sf Cov}(\mathbb{B}(x), \mathbb{B}(y)) = \sqrt{n}{\sf Cov}(\widehat F_n(x), \widehat F_n(y))$.

Theorem 5 shows that after rescaling, the difference between the estimator and the parameter of interest converges to a stable distribution (often is also a normal distribution). A formal statement of the above result can be found in Theorem 2.8 of Kosorok (2008). Moreover, we have the following result about the bootstrap.

Theorem 6 (Bootstrap Functional delta method). Under the assumption of Theorem 5, then \[%\sup_t\left|P\left(\sqrt{n}(\theta(\hat F^*_n) - \theta(\hat F_n))<t|X_1,\cdots, X_n\right)- P\left(\sqrt{n}(\theta(\hat F_n) - \theta(F_0))<t\right)\right| = o_P(1). \sqrt{n}(\theta(\widehat F^*_n) - \theta(\widehat F_n))\approx \dot{\theta}_{F_0}(\mathbb{B}),\] where $\dot{\theta}_{F_0}(\mathbb{B})$ is the same as the one in Theorem 5.

Theorem 6 shows that as long as the functional delta method works, the bootstrap method also works. The formal statement of Theorem 6 can be found in Theorem 2.9 of Kosorok (2008). Formally, we say that the bootstrap method is consistent if the random variable $\sqrt{n}(\theta(\widehat F^*_n) - \theta(\widehat F_n))$ converges in the Kolmogorov distance to $\dot{\theta}_{F_0}(\mathbb{B})$, i.e., \[\begin{equation} \sup_t\left|P\left(\sqrt{n}(\theta(\widehat F^*_n) - \theta(\widehat F_n))<t|X_1,\cdots, X_n\right)- P\left(\dot{\theta}_{F_0}(\mathbb{B})<t\right)\right| = o_P(1). \label{eq::Kdist} \end{equation}\] The Berry-Esseen bound is a common approach to establish the consistency of a bootstrap method.

Remark (higher-order accuracy). There are many ways of constructing a bootstrap confidence interval. So one may be wondering which method give the best confidence interval. Here is a common way of measuring how good a bootstrap confidence interval is using the idea of accuracy. In a bound like the one in equation $\ref{eq::Kdist}$, the $o_P(1)$ can often be explicitly written as $O_P(n^{-k})$. The Berry-Esseen bound gives an accuracy at rate $k=1/2$ and the percentile method also give the same rate. However, the bootstrap t-percentile method may lead to $k=1$. When the limit is a Gaussian distribution, a common approach of finding $k$ is via the Edgeworth expansion: we try to expand the CDF of $\sqrt{n}(\theta(\widehat F^*_n) - \theta(\widehat F_n))$ around a normal CDF using a Taylor expansion-like method. See Chapter 23.3 of Van der Vaart (2000) for more details.

Bootstrap and empirical process

In density estimator and regression analysis, we have seen that our estimator may be a function. So it will be of great interest to study if the bootstrap can be applied to a function estimator.

In this case, we need to introduce the concept of empirical processes. Let $\mathcal{F} = \{f_t: t\in \mathbb{T}\}$ be a collection of functions such that $f_t:\mathcal{X} \mapsto\mathbb{R}$, where $\mathcal{X}$ is the support of the observations and $\mathbb{T}$ is the index set. Consider the sample mean type estimator: \[\widehat{\mathbb{P}}_n(f_t) = \frac{1}{n}\sum_{i=1}^n f_t(X_i).\] It is easy to see that this is an unbiased estimator of \[\mathbb{P}(f_t) = \mathbb{E}(f_t(X_1)).\] The scaled difference is called the empirical process: \[\mathbb{G}_n (f_t) = \sqrt{n} (\widehat{\mathbb{P}}_n(f_t) - {\mathbb{P}}(f_t)).\] When we vary this quantity over $t\in \mathbb{T}$, you can see that $\mathbb{G}_n(f_t) = G(t)$ is a random function of the argument $t\in\mathbb{T}$. Although this seems to be abstract, here are some examples of empirical processes.

Example: EDF. Consider $\mathbb{T} = \mathcal{X}$ and \[f_t(x) = I(x\leq t).\] In this case, \[\widehat{\mathbb{P}}_n(f_t) = \frac{1}{n}\sum_{i=1}^n I(X_i\leq t) = \widehat F_n(t)\] is the EDF.

Example: KDE. Consider $\mathbb{T} = \mathcal{X}$ and \[f_t(x) = \frac{1}{h^d} K\left(\frac{x-t}{h}\right).\] Then it is easy to see that \[\widehat{\mathbb{P}}_n(f_t) = \frac{1}{nh^d}\sum_{i=1}^n K\left(\frac{X_i-t}{h}\right) = \widehat p_h(t),\] which is the KDE.

Example: log-likelihood function. Consider $\mathbb{T} = \Theta$ and \[f_\theta(x) = \ell(\theta|x),\] where $\ell(\theta|x) = \log p_\theta(x)$ is a log-likelihood function. Then \[\widehat{\mathbb{P}}_n(f_\theta) = \frac{1}{n}\sum_{i=1}^n \ell(\theta|X_i)\] is the joint log-likelihood function.

Example: empirical risk. Consider a prediction problem where we want to predict $Y$ based on $X$ and we assume that our predictor $m_\theta: \mathcal{X}\rightarrow \mathcal{Y}$ is a parametric model with parameter $\theta$. Let $L:\mathcal{Y}\times \mathcal{Y}\rightarrow \mathbb{R}$ be the loss function. Then the empirical risk \[R_n(\theta) = \frac{1}{n}\sum_{i=1}^n L(m_\theta(X_i), Y_i)=\widehat{\mathbb{P}}_n( f_\theta),\] where $f_\theta(x,y) = L(m_\theta(x), y)$.

When the estimator is a vector, we have central limit theorem showing that the estimator converges to a Gaussian vector after rescaling. A similar pattern occurs for the empirical process as well. We call $\mathcal{F} = \{f_t: t\in \mathbb{T}\}$ a Donsker class if \[\mathbb{G}_n(f_t)\overset{D}{\rightarrow} \mathbb{B}(f_t),\] where $\mathbb{B}(f_t)$ is a Gaussian process and the notation $\overset{D}{\rightarrow}$ stands for convergence in distribution (weak convergence) of a stochastic process under $L_{\infty}$ norm. Informally, we have the following result, known as the uniform central limit theorem.

In the formal statement, the approximation sign $\approx$ will be replaced by weak convergence of a stochastic process (notice that $\mathbb{G}_n^*(f_t)$ is a process indexed by $t$) conditioned on the sample $X_1,\cdots, X_n$. The take away message is that when the function class is Donsker, the bootstrap method works.

Roughly speaking, most common parametric models is a Donsker class as described in the following proposition.

Proposition 8 (Example 19.7 in van der Vaart (2000)). Let $\mathcal{F} = \{f_\theta: \theta\in \Theta\}$ be a collection of functions such that $\Theta$ is bounded subset of $\mathbb{R}^d$. Suppose that there exists a function $m$ such that \[|f_{\theta_1}(x) - f_{\theta_2}(x)| \leq m(x) \|\theta_1-\theta_2\|\] for any $\theta_1,\theta_2\in\Theta$ and $\mathbb{E}(|m(X)|^2)<\infty$. Then $\mathcal{F}$ is Donsker.

There are several ways to argue that a function class $\mathcal{F}$ is Donsker. Common approaches are based on bounding the uniform covering number or the bracketing number. The two textbooks mentioned in the above (Kosorok 2008 and van der Vaart & Wellner 1996) are good references. You will learn more about this in STAT 580 sequence.

The Donsker theory can be combined with the functional delta method. Roughly speaking, if the parameter of interest is a statistical functional such that $\theta (\mathbb{P}(\mathcal{F}))$, where $\mathbb{P}(\mathcal{F}) = \{\mathbb{P}(f_t): t\in \mathbb{T}\}$ is a stochastic process, then \[%\sup_t\left|P\left(\sqrt{n}(\theta(\hat{\mathbb{P}}_n^*(\mathcal{F})) - \theta(\hat{\mathbb{P}}_n(\mathcal{F})))<t|X_1,\cdots, X_n\right)- P\left(\sqrt{n}(\theta(\hat{\mathbb{P}}_n(\mathcal{F})) - \theta(\mathbb{P}(\mathcal{F})))<t\right)\right| = o_P(1) \sqrt{n}(\theta(\widehat{\mathbb{P}}_n^*(\mathcal{F})) - \theta(\widehat{\mathbb{P}}_n(\mathcal{F})))\approx \sqrt{n}(\theta(\widehat{\mathbb{P}}_n(\mathcal{F})) - \theta(\mathbb{P}(\mathcal{F}))),\] when $\theta$ is Hadamard differentiable at $\mathbb{P}(\mathcal{F})$.

Beyond Functional Delta Method and Donsker Class

The functional delta method along with the Donsker theory makes the bootstrap a widely applicable approach. However, there are cases where these methodologies cannot be applied. Here are two examples that the regular Donsker theory does not work.

High-dimensional models. Let $X_1,\cdots, X_n\in\mathbb{R}^d$ be IID. Suppose that we are interested in estimating the mean vector $\mu = \mathbb{E}(X_1)$. It is easy to see that the sample mean $\bar{X}_n$ has asymptotic normality when $d$ is fixed and the covariance matrix $\Sigma = \mathbb{E}(X_1X_1^T)$ is finite. However, when $d = d_n\rightarrow \infty$, we may not have the same asymptotic normality⁹. In particular, when $\log d \asymp n$ (regular high-dimensional model), there is no asymptotic normality of the sample mean so it is unclear how to use the bootstrap to construct a confidence set of $\mu$.

KDE with decreasing $h$. In the case of KDE, if $h$ is fixed, then the function class \[\left\{K\left(\frac{\cdot - x}{h}\right):x\in\mathcal{X}\right\}\] is a Donsker class. However, if $h=h_n\rightarrow0$, then the function class being considered \[\left\{K\left(\frac{\cdot - x}{h}\right):x\in\mathcal{X},1> h>0\right\}\] is no longer a Donsker class. Note that we use $1$ as the upper bound of $h$ since asymptotically the bandwidth will be less than $1$. It can be replaced by any finite upper bound.

Bootstrap and High-Dimensional Models

Let $X_1,\cdots, X_n\in\mathbb{R}^d$ be IID with $\mathbb{E}(X_i) = \mu$ and $\Sigma = \mathbb{E}(X_1X_1^T)$. In recent work of Chernozhukov and his collaborators, they showed that under good conditions (allowing $d =o( e^{n^{1/8}})$), \[\sqrt{\frac{n}{\log d}}\|\bar X_n- \mu\|_{\max } \approx \sqrt{\frac{n}{\log d}}\|\bar Z_n-\mu\|_{\max},\] where $\bar Z_n = \frac{1}{n}\sum_{i=1}^nZ_i$ is the mean of a standard normal vector. The $\approx$ here is actually under the Kolmogorov distance.

Moreover, the above result can be applied to the bootstrap sample mean $\bar X_n^*$, which leads to the following “rectangle confidence set". Let \[\widehat t_\alpha = \widehat G^{-1}(1-\alpha/2),\quad \widehat G(s) = P(\|\bar X_n^* - \bar X_n\|_{\max }\leq s|X_1,\cdots, X_n).\] Namely, $\widehat t_\alpha$ is the $1-\alpha$ quantile of the bootstrap maximum $\|\bar X_n^* - \bar X_n\|_{\max }$. Define a confidence set as \[C_n = \{\mu\in\mathbb{R}^d: \|\bar X_n-\mu\|_{\max}\leq \widehat t_\alpha\}.\] Note that $C_n$ looks like a rectangle in the parameter space.

Covariance matrix. Using the delta method, this result applies to other smoothed functionals of sample mean or higher moments as well. One notable result is the use of this idea to construct confidence sets of the population covariance matrix in high dimensions. Let $\widehat\Sigma_n$ be the sample covariance matrix and let $\widehat\Sigma^*_n$ be the bootstrap sample covariance matrix. Define \[\widehat\eta_\alpha = \widehat G^{-1}(1-\alpha/2),\quad \widehat G(s) = P(\|\widehat\Sigma_n^* - \widehat\Sigma_n\|_{\max }\leq s|X_1,\cdots, X_n)\] and two matrices \[U_\alpha = \widehat\Sigma_n +\widehat\eta_\alpha,\quad L_\alpha = \widehat\Sigma_n - \widehat\eta_\alpha.\] Then you can show that \[P(L_{\alpha,ij} \leq \Sigma_{ij} \leq U_{\alpha,ij},\,\, \forall i,\,j) \geq 1-\alpha -\frac{c_0 \log d}{n^{1/8}}\] for some constant. See the following paper for more details

Bootstrap and the Supremum of an Empirical Process

The convergence of a maximum of a Gaussian vector also implies that the supremum of a good empirical process can also be well-approximated by the supremum of a Gaussian process. This result works even for some non-Donsker classes. In particular, it can be applied to KDE.

Let $X_1,\cdots, X_n$ be IID from some distribution supported on $\mathcal{X}$ and $\widehat p_h(x) = \frac{1}{nh^d}\sum_{i=1}^n K\left(\frac{X_i-x}{h}\right)$ be the KDE and $\widehat p_h^*(x) = \frac{1}{nh^d}\sum_{i=1}^n K\left(\frac{X^*_i-x}{h}\right)$ be the bootstrap KDE. Also let $p_h(x) = \mathbb{E}(\widehat p_h(x))$ be the expected version of the KDE. Under the conditions in the following paper

we have the bootstrap validity of the confidence band of the KDE (actually, their result is stronger than the following theorem).

Theorem 10 (Chernozhukov, Chetverikov, and Kato (2014)). Assume conditions in Chernozhukov, Chetverikov, and Kato (2014). Then there exists a mean zero Gaussian process $\mathbb{B}$ defined over $\mathcal{X}$ such that \[\sup_t \left|P\left(\sqrt{\frac{nh^d}{\log n}} \|\widehat p_h - p_h\|_{\infty}<t\right) - P\left(\|\mathbb{B}\|_{\infty}<t\right)\right| = O_P\left(\left(\frac{\log n}{nh^{d}}\right)^{1/8}\right).\] Moreover, \[\sup_t \left|P\left(\sqrt{\frac{nh^d}{\log n}} \|\widehat p^*_h -\widehat p_h\|_{\infty}<t|X_1,\cdots, X_n\right) -P\left(\sqrt{\frac{nh^d}{\log n}} \|\widehat p_h - p_h\|_{\infty}<t\right)\right| = O_P\left(\left(\frac{\log n}{nh^{d}}\right)^{1/8}\right).\]

The conditions in Theorem 10 are quite mild – actually, they are the same conditions as the rate of uniform convergence. This Theorem also implies a construction of a simultaneous confidence band. Let \[\widehat t_\alpha = \widehat G(1-\alpha/2),\quad \widehat G(s) = P( \|\widehat p^*_h -\widehat p_h\|_{\infty}\leq s|X_1,\cdots, X_n)\] be the $1-\alpha$ quantile of the bootstrap supremum error $\|\widehat p^*_h -\widehat p_h\|_{\infty}$. We define \[L_\alpha = \widehat p_h - \widehat t_\alpha,\quad U_\alpha = \widehat p_h + \widehat t_\alpha.\] Then you can show that $L_\alpha,U_\alpha$ can be used as a simultaneous confidence band and we have \[P(L_\alpha(x)\leq p_h(x)\leq U_\alpha(x)\,\,\forall x) \geq 1-\alpha+ c_0 \left(\frac{\log n}{nh^d}\right)^{1/8}\] for some constant $c_0>0$.

Note that the above confidence band is simultaneous for $p_h$, not the true PDF $p$. So we have to undersmooth the KDE (i.e, choosing $h$ to be at a fast rate to $0$ than the optimal rate) to obtain a valid confidence band. One approach that can by pass this problem is via bootstrapping the debiased estimator. See the following paper

For more details, check wikipedia: https://en.wikipedia.org/wiki/Bootstrapping_(statistics)↩︎
Note that here we ignore lots of technical details. The smoothness of a ‘functional’ is an advanced topic in mathematics called functional analysis: https://en.wikipedia.org/wiki/Functional_analysis. There are formal ways of defining continuity of functionals and even ‘differentiation’ of functionals; see, e.g., https://en.wikipedia.org/wiki/G%C3%A2teaux_derivative.↩︎
A more formal way is to show that it converges in probability.↩︎
Note that here we use the continuous mapping theorem: if $f$ is a continuous function and random variable $A_n\overset{P}{\rightarrow} a_0$, then $f(A_n)\overset{P}{\rightarrow}f(a_0)$. Setting $f(x) = x^2$, we obtain the convergence of the second quantity.↩︎
https://en.wikipedia.org/wiki/G%C3%A2teaux_derivative.↩︎
More precisely, we need it to be Hadamard differentiable with respect to the $L_\infty$ metric $d(F,G) = \sup_x|F(x)-G(x)|$; see https://en.wikipedia.org/wiki/Hadamard_derivative ↩︎
https://en.wikipedia.org/wiki/Robust_statistics#Influence_function_and_sensitivity_curve ↩︎
https://en.wikipedia.org/wiki/Semiparametric_model ↩︎
We may still have CLT with $d\rightarrow\infty, d^{2}/n\rightarrow 0$, see Portnoy, S. (1984). Asymptotic behavior of M-estimators of p regression parameters when p 2/n is large. I. Consistency. The Annals of Statistics, 1298-1309.↩︎

Introduction