Bootstrap and Resampling Methods


In-depth Articles

Introduction to Bootstrap

The bootstrap, introduced by Efron (1979), is a computationally intensive method for statistical inference. The fundamental idea is simple yet powerful: if the sample is a good representation of the population, then resampling from the sample mimics drawing new samples from the population.

The core problem in statistics is estimating the sampling distribution of a statistic \(\hat{\theta} = T(\mathbf{X})\). Traditionally, we derive this distribution analytically (often relying on asymptotic approximations). The bootstrap replaces analytical derivation with computational simulation: we resample with replacement from the observed data to generate an empirical approximation of the sampling distribution.

The key insight is the plug-in principle: we estimate the unknown population distribution \(F\) with the empirical distribution function \(\hat{F}_n\), which places mass \(1/n\) on each observed data point:

\[ \hat{F}_n(x) = \frac{1}{n} \sum_{i=1}^{n} I(X_i \leq x) \]

Any quantity that depends on \(F\) can then be estimated by substituting \(\hat{F}_n\). The bootstrap automates this substitution for complex statistics where analytical derivations are intractable.




Nonparametric Bootstrap

Given an observed sample \(\mathbf{x} = (x_1, x_2, \ldots, x_n)\) from an unknown distribution \(F\), and a statistic of interest \(\hat{\theta} = T(\mathbf{x})\), the nonparametric bootstrap algorithm proceeds as follows:

Algorithm:

  1. For \(b = 1, 2, \ldots, B\):
    1. Draw a bootstrap sample \(\mathbf{x}^*_b = (x^*_1, x^*_2, \ldots, x^*_n)\) by sampling \(n\) observations with replacement from \(\hat{F}_n\) (i.e., from the original sample)
    2. Compute the bootstrap replicate \(\hat{\theta}^*_b = T(\mathbf{x}^*_b)\)
  2. Use the empirical distribution of \(\{\hat{\theta}^*_1, \hat{\theta}^*_2, \ldots, \hat{\theta}^*_B\}\) to approximate the sampling distribution of \(\hat{\theta}\)

The bootstrap estimate of the standard error of \(\hat{\theta}\) is:

\[ \widehat{SE}_{boot} = \sqrt{\frac{1}{B-1} \sum_{b=1}^{B} \left(\hat{\theta}^*_b - \bar{\theta}^*\right)^2} \]

where \(\bar{\theta}^* = \frac{1}{B} \sum_{b=1}^{B} \hat{\theta}^*_b\) is the mean of the bootstrap replicates.

The bootstrap estimate of bias is:

\[ \widehat{Bias}_{boot} = \bar{\theta}^* - \hat{\theta} \]

In practice, \(B = 200\) is often sufficient for standard error estimation, while \(B \geq 1000\) is recommended for confidence intervals.

Conditions for Validity

The nonparametric bootstrap is consistent (i.e., the bootstrap distribution converges to the true sampling distribution) under relatively mild conditions. However, it can fail when:

  • The statistic is not smooth (e.g., extremes like \(\max(X_1, \ldots, X_n)\))
  • The parameter lies on the boundary of the parameter space
  • The rate of convergence is not \(\sqrt{n}\)
  • The sample is not i.i.d. (time series, clustered data)



Parametric Bootstrap

When the nonparametric bootstrap fails — for instance when the parameter of interest influences the support of the distribution — or when a parametric model is believed to be correct, we can use the parametric bootstrap.

Instead of resampling from \(\hat{F}_n\), we:

  1. Assume a parametric family \(F_\theta\) for the data-generating process
  2. Estimate the parameters: \(\hat{\theta} = \hat{\theta}(\mathbf{x})\)
  3. Generate bootstrap samples from the fitted model \(F_{\hat{\theta}}\)

Algorithm:

  1. Fit the parametric model to the data: obtain \(\hat{\theta}\)
  2. For \(b = 1, 2, \ldots, B\):
    1. Generate \(\mathbf{x}^*_b = (x^*_1, \ldots, x^*_n)\) where \(x^*_i \stackrel{i.i.d.}{\sim} F_{\hat{\theta}}\)
    2. Compute \(\hat{\theta}^*_b = T(\mathbf{x}^*_b)\)
  3. Use \(\{\hat{\theta}^*_1, \ldots, \hat{\theta}^*_B\}\) to approximate the sampling distribution

Example: Consider estimating the endpoint \(\theta\) of a Uniform\((0, \theta)\) distribution using \(\hat{\theta} = X_{(n)} = \max(X_1, \ldots, X_n)\). The nonparametric bootstrap fails here because bootstrap samples from \(\hat{F}_n\) cannot exceed \(X_{(n)}\), so the bootstrap distribution of \(\hat{\theta}^*\) is degenerate at \(X_{(n)}\). The parametric bootstrap generates from Uniform\((0, \hat{\theta})\), correctly capturing the variability of the maximum.

The parametric bootstrap is more efficient than the nonparametric version when the model is correctly specified, but is inconsistent if the model is misspecified.




Bootstrap for Regression

Consider the linear regression model:

\[ Y_i = \mathbf{x}_i^\top \boldsymbol{\beta} + \varepsilon_i, \quad i = 1, \ldots, n \]

There are two distinct bootstrap strategies depending on whether the covariates are treated as random or fixed.

Random-X Resampling (Resample Cases)

Treat each pair \((Y_i, \mathbf{x}_i)\) as a single observation drawn from a joint distribution. The bootstrap resamples entire cases:

  1. For \(b = 1, \ldots, B\): draw \(n\) pairs \((Y^*_i, \mathbf{x}^*_i)\) with replacement from \(\{(Y_1, \mathbf{x}_1), \ldots, (Y_n, \mathbf{x}_n)\}\)
  2. Fit the regression on the bootstrap sample to obtain \(\hat{\boldsymbol{\beta}}^*_b\)

This approach:

  • Makes no assumptions about the error distribution or homoscedasticity
  • Is valid under heteroscedasticity
  • Is the default choice when \(\mathbf{x}\) is observational (not experimentally controlled)

Fixed-X Resampling (Resample Residuals)

When the design matrix \(\mathbf{X}\) is fixed (experimental design), we resample residuals:

  1. Fit the model: obtain \(\hat{\boldsymbol{\beta}}\) and residuals \(\hat{e}_i = Y_i - \mathbf{x}_i^\top \hat{\boldsymbol{\beta}}\)
  2. For \(b = 1, \ldots, B\): draw \(e^*_1, \ldots, e^*_n\) with replacement from \(\{\hat{e}_1, \ldots, \hat{e}_n\}\)
  3. Construct bootstrap responses: \(Y^*_i = \mathbf{x}_i^\top \hat{\boldsymbol{\beta}} + e^*_i\)
  4. Fit the regression on \((\mathbf{X}, \mathbf{Y}^*)\) to obtain \(\hat{\boldsymbol{\beta}}^*_b\)

This approach:

  • Preserves the design structure
  • Assumes errors are i.i.d. (homoscedastic)
  • Is more efficient than case resampling when the model is correctly specified
  • Fails under heteroscedasticity (use Wild Bootstrap instead)



The Jackknife

The jackknife, introduced by Quenouille (1949) and developed by Tukey (1958), is a resampling method based on systematically leaving out one observation at a time. It predates the bootstrap and can be seen as a linear approximation to it.

Given a sample \(\mathbf{x} = (x_1, \ldots, x_n)\) and statistic \(\hat{\theta} = T(\mathbf{x})\), the \(i\)-th jackknife replicate is computed on the sample with the \(i\)-th observation removed:

\[ \hat{\theta}_{(i)} = T(\mathbf{x}_{(-i)}) = T(x_1, \ldots, x_{i-1}, x_{i+1}, \ldots, x_n) \]

The jackknife mean is:

\[ \bar{\theta}_{(\cdot)} = \frac{1}{n} \sum_{i=1}^{n} \hat{\theta}_{(i)} \]

Bias Estimation

The jackknife estimate of bias is:

\[ \widehat{Bias}_{jack} = (n-1)\left(\bar{\theta}_{(\cdot)} - \hat{\theta}\right) \]

The bias-corrected jackknife estimator is:

\[ \tilde{\theta}_{jack} = \hat{\theta} - \widehat{Bias}_{jack} = n\hat{\theta} - (n-1)\bar{\theta}_{(\cdot)} \]

Variance Estimation

The jackknife estimate of variance is:

\[ \widehat{Var}_{jack} = \frac{n-1}{n} \sum_{i=1}^{n} \left(\hat{\theta}_{(i)} - \bar{\theta}_{(\cdot)}\right)^2 \]

Note the factor \(\frac{n-1}{n}\) instead of the usual \(\frac{1}{n-1}\) — this compensates for the high correlation between jackknife replicates (each pair shares \(n-2\) observations).

Jackknife Pseudo-Values

The pseudo-values are defined as:

\[ \tilde{\theta}_i = n\hat{\theta} - (n-1)\hat{\theta}_{(i)} \]

These pseudo-values are approximately independent and can be used as input to standard parametric procedures (e.g., t-tests).

Limitations

The jackknife fails for non-smooth statistics such as the median or quantiles. For such statistics, the delete-\(d\) jackknife (removing \(d > 1\) observations at a time) may be used, but the bootstrap is generally preferred.




Bootstrap Confidence Intervals

Let \(\hat{\theta}\) be an estimator of \(\theta\) and \(\{\hat{\theta}^*_1, \ldots, \hat{\theta}^*_B\}\) the bootstrap replicates. We seek a \((1-\alpha)\) confidence interval for \(\theta\). Several methods exist, with increasing sophistication.

(a) Asymptotic Normal Interval

The simplest approach assumes the sampling distribution of \(\hat{\theta}\) is approximately normal:

\[ CI_{1-\alpha} = \left[\hat{\theta} - z_{\alpha/2} \cdot \widehat{SE}_{boot}, \quad \hat{\theta} + z_{\alpha/2} \cdot \widehat{SE}_{boot}\right] \]

where \(z_{\alpha/2}\) is the \((1-\alpha/2)\) quantile of the standard normal and \(\widehat{SE}_{boot}\) is the bootstrap standard error.

This interval is first-order accurate: its coverage error is \(O(n^{-1/2})\). It fails when the sampling distribution is skewed or when the normal approximation is poor.

(b) Studentized (Pivotal) Bootstrap Interval

This method bootstraps the pivotal quantity (t-statistic) rather than the estimator itself. For each bootstrap replicate, compute:

\[ t^*_b = \frac{\hat{\theta}^*_b - \hat{\theta}}{\widehat{SE}^*_b} \]

where \(\widehat{SE}^*_b\) is the estimated standard error of \(\hat{\theta}^*_b\) (which may itself require a nested bootstrap or analytical formula). Let \(t^*_{(\alpha)}\) denote the \(\alpha\)-quantile of the bootstrap t-distribution. The confidence interval is:

\[ CI_{1-\alpha} = \left[\hat{\theta} - t^*_{(1-\alpha/2)} \cdot \widehat{SE}, \quad \hat{\theta} - t^*_{(\alpha/2)} \cdot \widehat{SE}\right] \]

This interval is second-order accurate: coverage error is \(O(n^{-1})\). It respects asymmetry in the sampling distribution and is transformation-invariant in the pivotal sense. However, it requires an estimate of \(\widehat{SE}^*_b\) for each replicate, making it computationally expensive.

(c) Basic Bootstrap Interval

The basic (or reverse percentile) method uses the bootstrap distribution of \((\hat{\theta}^* - \hat{\theta})\) to approximate the distribution of \((\hat{\theta} - \theta)\). The interval is:

\[ CI_{1-\alpha} = \left[2\hat{\theta} - \hat{\theta}^*_{(1-\alpha/2)}, \quad 2\hat{\theta} - \hat{\theta}^*_{(\alpha/2)}\right] \]

where \(\hat{\theta}^*_{(\alpha)}\) is the \(\alpha\)-quantile of the bootstrap distribution. Equivalently:

\[ CI_{1-\alpha} = \left[\hat{\theta} - (\hat{\theta}^*_{(1-\alpha/2)} - \hat{\theta}), \quad \hat{\theta} - (\hat{\theta}^*_{(\alpha/2)} - \hat{\theta})\right] \]

This interval reflects the bootstrap distribution about \(\hat{\theta}\). It is first-order accurate like the normal interval but does not require normality.

(d) Percentile Interval

The percentile method directly uses quantiles of the bootstrap distribution:

\[ CI_{1-\alpha} = \left[\hat{\theta}^*_{(\alpha/2)}, \quad \hat{\theta}^*_{(1-\alpha/2)}\right] \]

This is the simplest and most intuitive method. It is transformation-respecting: if \(\phi = g(\theta)\), then the percentile interval for \(\phi\) is \([g(\hat{\theta}^*_{(\alpha/2)}), g(\hat{\theta}^*_{(1-\alpha/2)})]\).

However, it is only first-order accurate and can have poor coverage when the bootstrap distribution is biased (i.e., not centred at \(\hat{\theta}\)).

(e) Bias-Corrected (BC) Interval

The BC interval corrects the percentile interval for bias in the bootstrap distribution. Define the bias-correction constant:

\[ \hat{z}_0 = \Phi^{-1}\left(\frac{\#\{\hat{\theta}^*_b < \hat{\theta}\}}{B}\right) \]

where \(\Phi^{-1}\) is the standard normal quantile function. If the bootstrap distribution is centred at \(\hat{\theta}\), then \(\hat{z}_0 = 0\) and BC reduces to the percentile method.

The BC interval adjusts the percentile levels:

\[ CI_{1-\alpha} = \left[\hat{\theta}^*_{(\alpha_1)}, \quad \hat{\theta}^*_{(\alpha_2)}\right] \]

where:

\[ \alpha_1 = \Phi\left(2\hat{z}_0 + z_{\alpha/2}\right), \qquad \alpha_2 = \Phi\left(2\hat{z}_0 + z_{1-\alpha/2}\right) \]

(f) BCa (Bias-Corrected and Accelerated) Interval

The BCa method extends BC by adding an acceleration constant \(\hat{a}\) that accounts for the rate at which the standard error of \(\hat{\theta}\) changes with \(\theta\). The adjusted percentile levels are:

\[ \alpha_1 = \Phi\left(\hat{z}_0 + \frac{\hat{z}_0 + z_{\alpha/2}}{1 - \hat{a}(\hat{z}_0 + z_{\alpha/2})}\right) \]

\[ \alpha_2 = \Phi\left(\hat{z}_0 + \frac{\hat{z}_0 + z_{1-\alpha/2}}{1 - \hat{a}(\hat{z}_0 + z_{1-\alpha/2})}\right) \]

The acceleration constant is typically estimated using the jackknife:

\[ \hat{a} = \frac{\sum_{i=1}^{n} \left(\bar{\theta}_{(\cdot)} - \hat{\theta}_{(i)}\right)^3}{6 \left[\sum_{i=1}^{n} \left(\bar{\theta}_{(\cdot)} - \hat{\theta}_{(i)}\right)^2\right]^{3/2}} \]

The BCa interval is:

  • Second-order accurate: coverage error is \(O(n^{-1})\)
  • Transformation-respecting: invariant under monotone reparametrizations
  • Range-preserving: endpoints respect the natural range of \(\theta\)

When \(\hat{z}_0 = 0\) and \(\hat{a} = 0\), BCa reduces to the standard percentile method. In practice, BCa is considered the best general-purpose bootstrap confidence interval.

Summary Comparison

MethodAccuracyTransform-RespectingRequires
Normal\(O(n^{-1/2})\)NoNormality of \(\hat{\theta}\)
Basic\(O(n^{-1/2})\)NoSymmetry of pivot
Percentile\(O(n^{-1/2})\)YesUnbiased bootstrap
Studentized\(O(n^{-1})\)NoSE estimate per replicate
BC\(O(n^{-1/2})\)YesBias correction \(\hat{z}_0\)
BCa\(O(n^{-1})\)Yes\(\hat{z}_0\) + acceleration \(\hat{a}\)