# Bootstrap vs numerical integration

My understanding of the bootstrap approach is based on Wasserman's framework (almost verbatim):

Let $T_n = g(X_1, ..., X_n)$ be a statistic ($X_i$ is the iid sample drawn from distribution $F$). Suppose we want to estimate $V_F(T_n)$ - the variance of $T_n$ given $F$.

The bootstrap approach follows these two steps:

1. Estimate $V_F(T_n)$ with $V_{\hat{F}}(T_n)$, where $\hat{F}$ is the empirical distribution function.

2. Approximate $V_{\hat{F}}(T_n)$ using simulation.

Do I understand correctly that the simulation in step 2 could be replaced with a precise calculation, except that it is infeasible for practically useful values of $n$? Here's my thinking: $V_{\hat{F}}$ precisely equals an integral of $T_n(X_1, ..., X_n)d\hat{F}(X_1)d\hat{F}(X_2)...d\hat{F}(X_n)$. $\hat{F}$ is a step-function, with a finite number $n$ steps; so we can ignore all points except the $n$ points where $d\hat{F}(x)$ has non-zero mass. So the integral is precisely equal to a sum of $n^n$ terms. Once $n$ exceeds 14, a simple direct calculation is impossible.

But all we're trying to do is calculate an integral. Why not replace the brute-force bootstrap simulation with any of the traditional numerical algorithms for taking integrals? Wouldn't it result in much higher precision for the same computational time?

Even something as simple as splitting the sample space in sections (perhaps with smaller volumes where sample statistic varies faster), and estimating the value of the statistic in each section by using the middle point, seems to be better than blind bootstrap.

What am I missing?

Perhaps bootstrap works so well and so fast that there's no need to do anything more complicated? (For example, if the loss of precision in step 1 is so much greater than in step 2, then improvements to step 2 are rather useless.)

The bootstrap works remarkably well. If you want to estimate the mean, variance and some not-too-extreme quantiles of the distribution of some low-dimensional $\hat\theta(Y)$, a few hundred to a few thousand re-samples will make the Monte Carlo error negligible, for many realistic problems. As a happy by-product, it also gives you a sample of $\hat\theta(Y^*)$, that can be used for diagnostic procedures, if desired, and it's not too difficult to get acceptably-good measures of how big the Monte Carlo errors actually are.

Fitting a regression model e.g. a thousand times over is (today) not a big deal, either in terms of CPU time or coding effort.

In contrast, numerical integration (excluding Monte Carlo methods) can be difficult to code - you'd have to decide how to split the sample space, for example, which is a non-trivial task. These methods also don't give the diagnostics, and the accuracy with which they estimate the true integral is notoriously hard to assess.

To do most of what the bootstrap does, but quicker, take a look at Generalized Method of Moments - for inference based on regression models (and much else) you can think of it as a quick, accurate approximation to what the non-parametric bootstrap would give.

• Thanks. Since step 2 is handled pretty well, I'm curious, can GMM or any other technique address the imprecision in step 1 (where we estimate the variance of the true distribution with the variance of the empirical distribution)? – max Mar 26 '12 at 17:11
• "Plain vanilla" GMM uses pretty straightforward approximations to the true covariance. Use of higher-order approximations (saddlepoint approximations and the like) can be used, but you'd have to code them yourself, and possibly make slightly-stronger assumptions than typical GMM to ensure you're getting the "best" approximation. – guest Mar 27 '12 at 6:46

The simulation most often used in bootstrapping for the numerical computation of the variance could in principle be replaced by an exact computation or an alternative approximation of the integral. One should, however, be aware that a "brute-force" simulation as an alternative to other numerical integration techniques is actually a good idea. The answer to the question "Wouldn't it result in much higher precision for the same computational time?" is no.

But why is that so? The thing is that standard numerical integration in high-dimensions scales badly with the dimension. If you are to divide the space into regular grid points, say, with $r$ grid points in each coordinate, you end up with $r^n$ grid points in total. The approximation achieved by simulation (know as Monte Carlo integration) can be viewed as a clever choice of function evaluations. Instead of time consuming grid evaluations we only evaluate the function we integrate at selected points. The error is, due to the random nature of the selected points, random, but can usually be controlled by the central limit theorem.

There are other methods such as quasi-Monte Carlo integration, that I know virtually nothing about, that make clever function evaluations based on quasi-random numbers instead of the pseudo-random numbers that we use for ordinary Monte Carlo integration.