Subsampling

statisticsresampling

Write \(P_F\) for probability under the data-generating distribution \(F\), \(\widehat F_n\) for the empirical distribution, and \(T_n\) for the root whose distribution is being approximated. Bootstrap probabilities condition on the observed sample. The setup and conditional convergence conventions fix the distinction between these two sources of randomness.

Subsampling approximates a sampling distribution by recomputing the statistic on smaller subsets of the observations. Its justification starts from a limit law for the statistic at the true population, so it can apply even when replacing that population by the empirical distribution is difficult to justify. For example, in the case where the data are \(n\) observations which are independent and identically distributed, a statistic \(\hat \theta_n\) is computed based on the entire data set and is recomputed over all \(\binom{n}{m}\) data sets of size \(m\). The same estimation rule must be defined at both sample sizes \(n\) and \(m\). These recomputed values of the statistic are suitably normalized to approximate the true sampling distribution. The consistency theorem below requires a limiting distribution and a separation between the full-sample and subsample scales. The bootstrap, on the other hand, requires that the distribution of the statistic is somehow locally smooth as a function of the unknown model. In contrast, no such assumption is required in the theory for subsampling. However, when both subsampling and the bootstrap are consistent, the bootstrap is typically more accurate.

Begin with independent observations and a real-valued parameter. Let \(X_1, \dots, X_n\) be a sample of \(n\) i.i.d. random variables taking values in an arbitrary sample space \(\mathcal S\). The common probability measure generating the observations is denoted by \(P\). The goal is to construct a confidence region for some parameter \(\theta(P)\). For now, assume \(\theta\) is real-valued, but this can and will be generalized to allow for the construction of confidence regions for multivariate parameters or confidence bands for functions. Let \(\hat \theta_n = \hat \theta_n (X_1, \dots, X_n)\) denote an estimator for \(\theta(P)\). Let

\[\begin{aligned} R_n(X^n, \theta(P)) = \tau_n (\hat \theta_n - \theta(P)) \end{aligned}\]

denote the root of interest, and let \(J_n(P)\) denote the sampling distribution of the root based on the sample \(X^n = (X_1, \dots, X_n)\) from \(P\). \(\tau_n\) is known, does not depend on \(P\), and is a normalizing constant. As before, the cumulative distribution function of \(J_n(P)\) is defined as

\[\begin{aligned} J_n (x,P) = P( \tau_n (\hat \theta_n - \theta(P)) \leq x ) ~. \end{aligned}\]

The first requirement is a limit for the normalized estimation error. This is the distribution that the subsample calculation will estimate when constructing confidence intervals for \(\theta(P)\).

Assumption 1. There exists a limiting distribution \(J(P)\) such that \(J_n (P) \rightsquigarrow J(P)\) as \(n \to \infty\) (for some sequence \(\tau_n\)).

The subsample statistic is centered at the full-sample estimate, since the population parameter is unknown. The scale condition below will make the error from this centering negligible.

Algorithm 2. Subsampling procedure

  1. Take all \(N_n=\binom{n}{m}\) subsets of size \(m\), and order these subsets in any way desired. Call them \(Y_1, \dots, Y_{N_n}\). Assume \(\frac{m}{n} \to 0\) as \(m, n \to \infty\).

    Define the approximation to \(J_n(x,P)\) as

    \[\begin{aligned} K_{n,m} (x) = \frac{1}{N_n} \sum_{1 \leq i \leq N_n} \mathbb 1 \{ \tau_m (\hat \theta_{n,m,i} - \hat \theta_n) \leq x \} ~, \end{aligned}\]

    where \(\hat \theta_{n,m,i}\) is the statistic \(\hat \theta_m\) evaluated using the subsample \(Y_i\).

Note that as the result of this procedure, we have, for any \(i\), \(Y_i\) is actually a random sample of \(m\) i.i.d. observations from \(P\). Hence, the exact distribution of \(\tau_m (\hat \theta_{n,m,i} - \theta (P))\) is \(J_m (P)\). We hope that the empirical distribution of the \(N_n\) values of this exact distribution serves as a good approximation to \(J_n (P)\), under some assumptions. This leads us to the following consistency result for this approximation procedure.

Theorem 3. Assume the estimator is a measurable symmetric function of its observations and that the normalizing constants are positive. Assume \(\frac{m}{n} \to 0\) and \(\frac{\tau_m}{\tau_n} \to 0\) as \(m,n \to \infty\). Furthermore, suppose Assumption 1 holds. Then,

  1. if \(x\) is a continuity point of \(J(\cdot, P)\), we have \(K_{n,m} (x) \to J(x,P)\) in probability.

  2. if \(J(\cdot, P)\) is continuous, we have \(K_{n,m} (x)\) converges uniformly in probability to \(J_n(x,P)\), i.e.

    \[\begin{aligned} \sup_x |K_{n,m} (x) - J_n (x,P)| \xrightarrow{P} 0 ~. \end{aligned}\]
  3. Define

    \[\begin{aligned} c_{n,m} (1-\alpha) &\triangleq \inf \{ x : K_{n,m}(x) \geq 1-\alpha \} ~, \\ c (1-\alpha, P) &\triangleq \inf \{ x : J(x,P) \geq 1-\alpha \} ~. \end{aligned}\]

    Suppose \(J(\cdot, P)\) is continuous (and strictly increasing) at \(c(1-\alpha,P)\). Then, the asymptotic coverage probability under \(P\) of the confidence interval \([\hat \theta_n - \tau_n^{-1} c_{n,m} (1-\alpha), \infty)\) is the nominal level \(1-\alpha\), i.e.

    \[\begin{aligned} P( \tau_n (\hat \theta_n - \theta(P)) \leq c_{n,m} (1-\alpha) ) \to 1-\alpha ~~ \text{ as } n \to \infty ~. \end{aligned}\]

Proof. First center every subsample at the true parameter and write

\[U_n(x)=\binom nm^{-1}\sum_{|A|=m}\mathbf1\{\tau_m(\widehat\theta_m(X_A)-\theta(P))\le x\}.\]

Let \(\Pi\) be a uniform random permutation independent of the data, put \(k=\lfloor n/m\rfloor\), and form the average over its disjoint blocks,

\[V_n(x)=\frac1k\sum_{j=1}^k\mathbf1\{\tau_m(\widehat\theta_m(X_{\Pi(m(j-1)+1)},\ldots,X_{\Pi(mj)})-\theta(P))\le x\}.\]

Each block is a uniform subset after conditioning on the data. Symmetry of the estimator therefore gives \(\mathbb{E}[V_n(x)\mid X_1,\ldots,X_n]=U_n(x)\). Unconditionally, the permuted observations remain iid, so the block indicators are independent Bernoulli variables with mean \(J_m(x,P)\). Conditional Jensen now gives the required direction:

\[\operatorname{Var}(U_n(x))\le\operatorname{Var}(V_n(x))\le\frac1{4k}.\]

The conditioning identity can also be seen by counting permutations. For each fixed block position, every subset \(A\) of size \(m\) occurs in exactly \(m!(n-m)!\) of the \(n!\) permutations. Hence its conditional probability is \(\binom nm^{-1}\), independently of the block position. Averaging over the \(k\) positions gives \(E[V_n(x)\mid X]=U_n(x)\). Therefore

\[E\{U_n(x)-J_m(x,P)\}^2=E\{E[V_n(x)-J_m(x,P)\mid X]\}^2\le E\{V_n(x)-J_m(x,P)\}^2.\]

The many overlapping subsamples need not be independent. Independence is used only for the disjoint blocks before taking their conditional average.

Since \(m/n\to0\), \(k\to\infty\). Thus \(U_n(x)-J_m(x,P)\to_P0\), and Assumption 1 gives \(U_n(x)\to_PJ(x,P)\) at every continuity point.

The difference between the two centerings is \(D_n=\tau_m(\widehat\theta_n-\theta(P))=(\tau_m/\tau_n)R_n=o_P(1)\), because the limiting law makes \(R_n\) tight. Exactly, \(K_{n,m}(x)=U_n(x+D_n)\). On \(\{|D_n|\le\epsilon\}\),

\[U_n(x-\epsilon)\le K_{n,m}(x)\le U_n(x+\epsilon).\tag{1}\]

Choose \(x\pm\epsilon\) among continuity points and let \(\epsilon\downarrow0\); this proves part (a). If \(J\) is continuous, a finite grid whose successive cdf increments are small, together with monotonicity and tail bounds, upgrades pointwise convergence in probability to uniform convergence. The deterministic cdfs \(J_n\) also converge uniformly to \(J\). The triangle inequality proves part (b).

Here is the probability estimate behind part (a). Given \(\eta>0\), choose continuity points \(x\pm\epsilon\) so close to \(x\) that both \(|J(x\pm\epsilon,P)-J(x,P)|<\eta/2\). Then

\[\begin{aligned}P(|K_{n,m}(x)-J(x,P)|>\eta)&\le P(|D_n|>\epsilon)\\&\quad+P(|U_n(x-\epsilon)-J(x-\epsilon,P)|>\eta/2)\\&\quad+P(|U_n(x+\epsilon)-J(x+\epsilon,P)|>\eta/2)\longrightarrow0.\end{aligned}\]

For part (b), choose a finite grid spanning all but \(\eta\) of each limiting tail, with successive limiting cdf increments at most \(\eta\). The union bound makes the maximum error on this finite grid tend to zero in probability; monotonicity bounds every intervening error by the grid error plus \(\eta\). Thus no uncountable intersection of probability-one events is required.

For part (c), let \(c=c(1-\alpha,P)\). Continuity and strict increase at this quantile mean that, for every \(\epsilon>0\), \(J(c-\epsilon,P)<1-\alpha<J(c+\epsilon,P)\) and \(J(c,P)=1-\alpha\). If either endpoint is a discontinuity, choose nearby continuity points still on the respective sides of the quantile. Cdf convergence at these continuity points traps the empirical quantile between them with probability tending to one, so \(c_{n,m}(1-\alpha)\to_Pc\). The limit law of \(R_n\), continuity at \(c\), and Slutsky's theorem then give the stated coverage. □

← All posts