Empirical Distribution Function

Before introducing the bootstrap method, we first introduce the empirical distribution function (EDF), an estimator of the cumulative distribution function (CDF).

Let us first examine the CDF \(F(x)\) more closely. Given a value \(x_0\), \[F(x_0) = P(X_i\leq x_0)\] for every \(i=1,\cdots, n\). Namely, \(F(x_0)\) is the probability of the event \(\{X_i\leq x_0\}\).

A natural estimator for the probability of an event is the relative frequency of that event within our sample. Thus, we use \[\begin{equation} \widehat{F}_n(x_0)= \frac{\mbox{number of $X_i$ $\leq x_0$}}{\mbox{total number of observations}} = \frac{\sum^n_{i=1}I(X_i\leq x_0)}{n} = \frac{1}{n}\sum^n_{i=1}I(X_i\leq x_0) \end{equation}\] as the estimator of \(F(x_0)\).

For every \(x_0\), we can use such a quantity as an estimator, so the estimator of the CDF, \(F(x)\), is \(\widehat{F}_n(x)\). This estimator, \(\widehat{F}_n(x)\), is called the empirical distribution function (EDF).

Example. Here is an example of the EDF of 5 observations of \(1, 1.2, 1.5,2,2.5\):

There are 5 jumps, each located at the position of an observation. Moreover, the height of each jump is the same: \(\frac{1}{5}\).

Example. While the previous example might not visually resemble an idealized CDF, the following provides a comparison of the EDF versus the true CDF when generating \(n=100\) and \(1000\) random samples from a standard normal distribution \(N(0,1)\):

The red curve indicates the true CDF of the standard normal. Here you can see that when the sample size is large, the EDF is pretty close to the true CDF.

Property of EDF

Because EDF is the average of \(I(X_i\leq x)\), we now study the property of \(I(X_i\leq x)\) first. For simplicity, let \(Y_i = I(X_i\leq x)\). What is the random variable \(Y_i\)?

Here is the breakdown of \(Y_i\): \[Y_i = \begin{cases} 1,\quad \mbox{if $X_i\leq x$}\\ 0,\quad \mbox{if $X_i> x$} \end{cases}.\] So \(Y_i\) only takes value \(0\) and \(1\)–so it is actually a Bernoulli random variable! We know that a Bernoulli random variable has a parameter \(p\) that determines the probability of outputing \(1\). What is the parameter \(p\) for \(Y_i\)? \[p = P(Y_i=1) = P(X_i\leq x) = F(x).\] Therefore, for a given \(x\), \[Y_i \sim {\sf Ber} (F(x)).\] This implies \[\begin{align*} \mathbb{E}(I(X_i\leq x)) &= \mathbb{E}(Y_i) = F(x)\\ {\sf Var}(I(X_i\leq x))& = {\sf Var}(Y_i) = F(x)(1-F(x)) \end{align*}\] for a given \(x\).

Now what about \(\widehat{F}_n(x)\)? Recall that \(\widehat{F}_n(x) = \frac{1}{n}\sum^n_{i=1}I(X_i\leq x) = \frac{1}{n}\sum_{i=1}^n Y_i\). Then \[\begin{align*} \mathbb{E}\left(\widehat{F}_n(x)\right) &= \mathbb{E}(I(X_1\leq x)) = F(x)\\ {\sf Var}\left(\widehat{F}_n(x)\right)& = \frac{\sum_{i=1}^n{\sf Var}(Y_i)}{n^2} = \frac{F(x)(1-F(x))}{n}. \end{align*}\] What does this tell us about using \(\widehat{F}_n(x)\) as an estimator of \(F(x)\)?

First, at each \(x\), \(\widehat{F}_n(x)\) is an unbiased estimator of \(F(x)\): \[{\bf bias}\left(\widehat{F}_n(x)\right) = \mathbb{E}\left(\widehat{F}_n(x)\right)-F(x) = 0.\] Second, the variance converges to \(0\) when \(n\rightarrow \infty\). By Lemma 0.3, this implies that for a given \(x\), \[\widehat{F}_n(x) \overset{P}{\rightarrow} F(x).\] i.e., \(\widehat{F}_n(x)\) is a consistent estimator of \(F(x)\).

In addition to the above properties, the EDF also have the following interesting feature: for a given \(x\), \[\sqrt{n}\left(\widehat{F}_n(x)-F(x)\right) \overset{D}{\rightarrow} N(0, F(x)(1-F(x))).\] Namely, \(\widehat{F}_n(x)\) is asymptotically normally distributed around \(F(x)\) with variance \(F(x)(1-F(x))\).

In fact, the asymptotic behavior extends beyond simple pointwise asymptotic normality. The scaled difference converges weakly to a Brownian bridge process, i.e., \[\sqrt{n}\left(\widehat{F}_n-F\right) \overset{D}{\rightarrow} \mathbb{B},\] where \(\mathbb{B}\) is a stochastic process called a Brownian bridge and \(\overset{D}{\rightarrow}\) here stands for weak convergence of a stochastic process.

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\) that takes a cumulative distribution function (CDF) as its argument and returns a real number.

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

• Mean of a distribution. The mean of a distribution is a statistical functional \[\mu = T_{\sf mean}(F) = \int x dF(x).\] When \(F\) has a PDF \(p(x)\), \(dF(x) = p(x ) dx\) so the mean functional reduces to the form that we are familiar with: \[\mu = T_{\sf mean}(F) = \int x dF(x) = \int x p(x) dx.\] When \(F\) is a distribution of discrete random variables, we define \[\int x dF(x) = \sum_{x} x P(x)\Longrightarrow \mu = T_{\sf mean}(F)= \sum_{x} x P(x),\] where \(P(x)\) is the PMF of the distribution \(F\).

  You may have noticed that if a random variable \(X\) has a CDF \(F\), then \[\mathbb{E}(X) = \int x dF(x) = T_{\sf mean}(F).\] Therefore, for any function \(g\), \[\mathbb{E}(g(X)) = \int g(x)dF(x).\] Using the function \(g\), we introduce another functional \(T_\omega\) such that \[T_\omega(F) = \int \omega(x) dF(x).\] Such a functional, \(T_\omega\), is called a linear functional.

• Variance of a distribution. The variance of a distribution is also a statistical functional. Let \(X\) be a random variable with CDF \(F\). Then \[\sigma^2 = T_{\sf var}(F) = {\sf Var}(X) = \mathbb{E}(X^2) - \mathbb{E}^2(X) = \int x^2 dF(x) - \left(\int xdF(x)\right)^2.\]

• Median of a distribution. Using the concept of a statistical functional, median and any quantile can be easily defined. The median of a distribution \(F\) is a point \(\theta_{\sf med}\) such that \(F(\theta_{\sf med})= 0.5\). Thus, \[T_{\sf med}(F) = F^{-1}\left(0.5\right).\] Note that when \(F\) is a CDF of a discrete random variable, \(F^{-1}\) may have multiple values. In this case, we define \[F^{-1}(q) = \inf\{x: F(x)\geq q\}.\] Any quantile of a distribution can be represented in a similar way. For instance, the \(q\)-quantile (\(0< q < 1\)) will be \[T_{\sf q}(F) = F^{-1}\left(q\right).\] As a result, the interquartile range (IQR) is \[T_{\sf IQR}(F) = F^{-1}(0.75)-F^{-1}(0.25).\]

What is the motivation for employing statistical functionals? Primarily, they elegantly formalize population parameters of interest. 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 aforementioned advantages, statistical functionals possess a powerful feature: they naturally motivate a simple class of estimators for 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 sample mean! Note that in the preceding calculation, we exploit the fact that \(\widehat{F}_n(x)\) corresponds to a discrete distribution placing equal probability mass (\(1/n\)) at each observation \(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 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, the statistical functional provides an elegant way to define a 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 because we know that \(\widehat{F}_n \rightarrow F\) in various ways so that the smoothness of \(T\) with respect to its input implies that \(T(\widehat{F}_n)\rightarrow T(F)\)¹.

Delta Method

In this section, we will talk about a very useful technique in handling the convergence–the delta method.

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 utilized when we derived the MLE of the rate parameter for 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}.\]

Now using equation \(\eqref{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}.\]

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 \omega(X_i). %T_\omega(\hat{F}^*_n) &= \int \omega(x) d\hat{F}^*_n(x)= \frac{1}{n}\sum_{i=1}^n \omega(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\).

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)\]

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 Number (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 the equation \(\eqref{eq::BT4}\) holds.

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}\] Equation \(\eqref{eq::IF1}\) formally represents a directional derivative. Indeed, it is the derivative of a functional with respect to a probability measure. This type of derivative is called Gâteaux derivative², a type of derivative of functionals. You can check that applying equation\(\eqref{eq::IF1}\) to a linear functional leads to an influence function as we defined previously.

A powerful feature of this generalized influence function is that when the statistical functional \(T_{\sf target}\) is sufficiently smooth³’, equation \(\eqref{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}\]

Example: median. Why is the sample median asymptotically normally distributed? 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). Note that \(F^{-1}(0.5) = T_{\sf med}(F)\) is the median of \(F\). So this shows not only the asymptotic normality of 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.

Empirical Bootstrap

The bootstrap is a powerful tool for assessing the uncertainty of an estimate. It can be used in many complex scenarios such as assessing the uncertainty of sample ‘median’. Here is how we can estimate the error of sample median and construct the corresponding confidence interval. Assume we are given the data points \(X_1,\cdots,X_n\). Let \(M_n = {\sf median} \{X_1,\cdots, X_n\}\). 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 sample procedure again, generating a new sample from the original dataset \(X_1,\cdots,X_n\) by sampling with replacement, leading to another new sets 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 totally, 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)⁶.

Now for each set of data, we then compute their sample median. This leads to \(B\) sample medians, called bootstrap medians: \[\begin{align*} M^{*(1)}_n &= {\sf median} \{X^{*(1)}_1,\cdots, X^{*(1)}_n\}\\ M^{*(2)}_n &= {\sf median} \{X^{*(2)}_1,\cdots, X^{*(2)}_n\}\\ &\vdots\\ M^{*(B)}_n &= {\sf median} \{X^{*(B)}_1,\cdots, X^{*(B)}_n\}. \end{align*}\]

This empirical bootstrap procedure facilitates several key statistical techniques.

• Bootstrap estimate of the variance. We will use the sample variance of \(M^{*(1)}_n, \cdots, M^{*(B)}_n\) as an estimate of the variance of sample median \(M_n\). Namely, we will use \[\widehat{{\sf Var}}_B(M_n) = \frac{1}{B-1}\sum_{\ell=1}^B\left(M^{*(\ell)}_n- \bar{M}^*_B\right)^2, \quad \bar{M}^*_B = \frac{1}{B }\sum_{\ell=1}^B M^{*(\ell)}_n,\] as an estimate of \({\sf Var}(M_n)\).

• Bootstrap estimate of the MSE. Moreover, we can estimate the MSE by \[\widehat{{\sf MSE}(M_n)} = \frac{1}{B}\sum_{\ell=1}^B \left(M^{*(\ell)}_n-M_n\right)^2.\]

• Bootstrap confidence interval. In addition, we can construct a \(1-\alpha\) confidence interval of the population median via \[M_n \pm z_{1-\alpha/2}\cdot \sqrt{\widehat{{\sf Var}}_B(M_n)}.\]

At first glance, this procedure may appear circular—generating new data by resampling from the existing dataset. However, under some conditions, this approach does work! And here is a brief explanation on why this approach works.

Let \(X_1,\cdots, X_n\sim F\). Recall from Lecture 1, 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 median \(M_n\), denoted as \(F_{M_n}\), will also be determined by the CDF \(F\) and sample size \(n\). Namely, we may write the CDF of median as \[\begin{equation} F_{M_n}(x) = \Psi(x; F,n), \label{eq::BT} \end{equation}\] where \(\Psi\) is some complicated function that depends on CDF of each observation \(F\) and the sample size \(n\).

When we sample with replacement from \(X_1,\cdots, X_n\), what is the underlying distribution from which we are sampling? 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 functions 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\).

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 median, say \(M^{*(1)}_n\), is the sample median of \(X_1^{*(1)},\cdots,X_n^{*(1)}\). Its CDF, by equation \(\eqref{eq::BT}\), is \[F_{M^{*(1)}_n}(x) = \Psi(x; \widehat{F}_n,n).\] And because each of the bootstrap sample are all from the distribution \(\widehat{F}_n\), we will have \[\Psi(x; \widehat{F}_n,n) = F_{M^{*(1)}_n}(x)=F_{M^{*(2)}_n}(x)=\cdots=F_{M^{*(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_{M^{*(\ell)}_n}(x)=\Psi(x; \widehat{F}_n,n)\approx \Psi(x; F,n) = F_{M_n}(x).\] This means:

The CDF of a bootstrap median, \(F_{M^{*(\ell)}_n}(x)\), is approximating the CDF of the true median, \(F_{M_n}(x)\).

This has many implications. For an example, when two CDFs are similar, their variances will be similar as well, i.e., \[{\sf Var}\left(M^{*(\ell)}_n|X_1,\cdots,X_n\right)\approx {\sf Var} (M_n).\footnote{ The reason why in 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(M_n)\) is just a sample variance of \(M^{*(\ell)}\). When \(B\) is large, the sample variance is about the same as the population variance, implying \[\widehat{{\sf Var}}_B(M_n) = \frac{1}{B-1}\sum_{\ell=1}^B\left(M^{*(\ell)}_n- \bar{M}^*_B\right)^2\approx {\sf Var}\left(M^{*(\ell)}_n|X_1,\cdots,X_n\right).\] Therefore, \[\widehat{{\sf Var}}_B(M_n)\approx{\sf Var}\left(M^{*(\ell)}_n|X_1,\cdots,X_n\right)\approx {\sf Var} (M_n),\] which explains why the bootstrap variance is a good estimate of the true variance of the median.

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.

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 why the bootstrap fails. Let \(X_1,\cdots,X_n\sim {\sf Uni}[0,1]\) and \(M_n = \min \{X_1,\cdots, X_n\}\) be the minimum value of the sample. Then it is known that \[n \cdot M_n \overset{D}{\rightarrow} {\sf Exp}(1).\] \(\spadesuit:\) Think about why it converges to exponential distribution.

Thus, \(M_n\) has a continuous distribution. Assume we generate a bootstrap sample \(X^*_1,\cdots,X^*_n\) from the original observations. Now let \(M_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 = M_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 \(M_n\) in the bootstrap sample is \[\begin{align*} P(\mbox{none of }X^*_1,\cdots,X^*_n\mbox{ select }M_n) = \left(1-\frac{1}{n}\right)^n \approx e^{-1}. \end{align*}\] This implies that with a probability \(1-e^{-1}\), one of the observation in the bootstrap sample will select the minimum value of the original sample \(M_n\). Namely, \[P(M^*_n = M_n) = 1-e^{-1}.\] Thus, \(M^*_n\) has a huge probability mass at the value \(M_n\), meaning that the distribution of \(M^*_n\) will not be close to an exponential distribution.

Bootstrap and Statistical Functionals

So far, we have not yet talked about the bootstrap. However, 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 \(\eqref{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\)).

Therefore, we conclude that

as long as we can show that equation \(\eqref{eq::BT3}\) holds, the bootstrap variance is a good estimate of the variance of the estimator \(\widehat{\theta}_n\).

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 \(\eqref{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 \(\eqref{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 Number, \[\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 \(\eqref{eq::BT4}\) holds and so is equation \(\eqref{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.

Berry-Esseen Bound (optional)

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 exist rigorous analytical routes to establish convergence of this form. One famous result is the Berry-Esseen bound of the sample mean.

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

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

The Berry-Esseen bound can be used to derive bounds like equation \(\eqref{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\).

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}}.\] By strong law of large number, 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.