# Bootstrap, Monte Carlo

I have been set the following question as part of homework:

Design and implement a simulation study to examine the performance of the bootstrap for obtaining 95% confidence intervals on the mean of a univariate sample of data. Your implementation can be in R or SAS.

Aspects of performance you may want to look at are confidence interval coverage (i.e., what proportion of times the confidence interval contains the true mean) and Monte Carlo variation (i.e., how much do the upper and lower confidence limits vary between simulations)'

Does anyone know how to go about the Monte Carlo variation aspect of this? I can't seem to even work out an algorithm or anything. Is it to do with Monte Carlo integration? Thanks!

This confusion between bootstrap procedures and Monte Carlo procedures keeps recurring, so perhaps this is as good a place as any to address it. (The examples of R code may help with the homework, too.)

Consider this implementation of the bootstrap in R:

boot <- function(x, t) { # Exact bootstrap of procedure t on data x
n <- length(x)       # Must lie between 2 and 7 inclusive.
if (n > 7) {
stop("Sample size exceeds 7; use an approximate method instead.")
}
p <- c(n, 1:(n-1))
a <- rep(x, n^(n-1))
dim(a) <- rep(n, n)
y <- as.vector(a)
while (n > 1) {
n <- n-1
a <- aperm(a, p)
y <- cbind(as.vector(a), y)
}
apply(y, 1, t)
}


A quick look will confirm that this is a deterministic calculation: no random values are generated or used. (I will leave the details of its inner workings for interested readers to figure out for themselves.)

The arguments to boot are a batch of numeric data in the array x and a reference t to a function (that can be applied to arrays exactly like x) to return a single numeric value; in other words, t is a statistic. It generates all possible samples with replacement from x and applies t to each of them, thereby producing one number for each such sample: that's the bootstrap in a nutshell. The return value is an array representing the exact bootstrap distribution of t for the sample x.

As a tiny example, let's bootstrap the mean for a sample x = c(1,3):

> boot(c(1,3), mean)
> [1] 1 2 2 3


There are indeed four possible samples of size $2$ with replacement from $(1,3)$; namely, $(1,1)$, $(1,3)$, $(3,1)$, and $(3,3)$. boot generates them all (in the order just listed) and applies t to each of them. In this case t computes the mean and those turn out to be $1$, $2$, $2$, and $3$, respectively, as shown in the output.

Where you go from here depends on how you want to use the bootstrap. The full information about the bootstrap is contained in this output array, so it's usually a good idea to display it. Here is an example where the standard deviation is bootstrapped from the sample $(1,3,3,4,7)$:

hist(boot(c(1,3,3,4,7), sd))


Now we are prepared to talk about Monte Carlo simulation. Suppose, say, we were going to bootstrap a 95% upper confidence limit on the SD from a sample of $5$ by using the upper 95th percentile of its bootstrap distribution. What properties would this procedure have? One way to find out is to suppose the sample were obtained randomly from, say, a uniform distribution. (The application will often indicate what a reasonable distributional assumption may be; here, I arbitrarily chose one that is simple for computation but not easy to deal with analytically.) We can simulate what happens by taking such a sample and computing the UCL:

> set.seed(17)
> quantile(boot(runif(5, min=0, max=10), sd), .95)[1]
95%
3.835870


The result for this particular random sample is 3.83587. This is definite: were you to call boot again with the same set of data, the answer would be exactly the same. But how might the answer change with different random samples? Find out by repeating this process a few times and drawing a histogram of the results:

> boot.sd <- replicate(100, quantile(boot(runif(5, min=0, max=10), sd), .95)[1])
> hist(boot.sd)


Were we to do another set of simulations, the random draws would come out different, creating a (slightly) different histogram--but not very different from this one. We can use it with some confidence to understand how the bootstrap UCL of the SD is working. For reference, notice that the standard deviation of a uniform distribution (spanning the range from $0$ to $10$ as specified here) equals $10/\sqrt{12} \approx 2.887$. As one would hope for any UCL worth its salt, the majority (three-quarters, or 0.75) of the values in the histogram exceed this:

> length(boot.sd[boot.sd >= 10/sqrt(12)]) / length(boot.sd)
[1] 0.75


But that's nowhere near the nominal 95% we specified (and were hoping for)! This is one value of simulation: it compares our hopes to what is really going on. (Why the discrepancy? I believe it is because bootstrapping an SD does not work well with really small samples.)

### Review

• Bootstrap statistics are conceptually just the same as any other statistic like a mean or standard deviation; they just tend to take a long time to compute. (See the warning message in the boot code!)

• Monte-Carlo simulation can be useful for studying how a bootstrap statistic varies due to randomness in obtaining samples. The variation observed in such simulation is due to variation in the samples, not variation in the bootstrap.

• (Not illustrated here) Because bootstrap statistics can take a lot of computation (apparently, up to $n^n$ calculations for samples of size $n$), it is convenient to approximate the bootstrap distribution. This is usually done by creating a "black box" program to obtain one value at random from the true bootstrap distribution and calling that program repeatedly. The collective output approximates the exact distribution. The approximation can vary due to randomness in the black box--but that variation is an artifact of the approximation procedure. It is not (conceptually) inherent in the bootstrap procedure itself.

• +1 This is very nice. I don't think I could have done it better! I would just add that I have some papers on the performance of the bootstrap in error rate estimation and for confidence intervals when estimating variances. These results can be found in the second edition of my bootstrap methods book and my book on bootstrapping with applications to R. – Michael Chernick Jul 24 '12 at 16:15
• Amazon links to the books are here: amazon.com/… amazon.com/Introduction-Bootstrap-Methods-Applications/dp/… – Michael Chernick Jul 24 '12 at 16:17
• I would just say that the idea is to simulate data from known distributions apply the bootstrap procedure of interest (e.g. confidence intervals for a parameter or estimation of a classifier's error rate) and you can evaluate how well it does because you know the truth. – Michael Chernick Jul 24 '12 at 16:20

The bootstrap is a Monte Carlo technique, in that it involves some type of random sampling. If you run the bootstrap twice on the same set of data, you'll get different answers. The more samples you use in your bootstrap, the less variation you'll get.

The coverage refers to the variation of the results on different data sets from the same hypothetical sampling distribution.

• thanks, I'll assume it just means to come up with code to see if the confidence limits vary between the simulations then! – Sarah Apr 15 '12 at 16:37
• To be precise (pedantic?), we should clarify that the bootstrap is not per se a Monte Carlo technique. However, it is usually approximated by using Monte Carlo simulation since it is not often tractable to calculate the particular quantities of interest. I find the wording of the problem very prone to misinterpretation as well! One must be careful to distinguish sampling error from Monte Carlo approximation error. The latter is only sensibly estimated by repeating the MC bootstrap approximation many times on the same sample, as you point out. :) – cardinal Apr 15 '12 at 16:46

I'm not sure what is meant exactly by "Monte Carlo variation" per se, either. Certainly it should be possible to look at how much variation there is amongst iterations in things like the value of the upper (or lower) limit, for example (hint). Perhaps they only want you to do this for the Monte Carlo, and not the bootstrap? That isn't a requirement I would have for an exercise, though. It may be best to just ask what is meant by that phrase.

• Thanks for your help! Perhaps I will just look at the variation between upper and lower limits for now and ask about Monte Carlo later. Now just to work out how I go about that.. – Sarah Apr 15 '12 at 16:11
• There are several possibilities here. One is upper-lower (ie, the difference). My interpretation of the phrasing in the question is variation in upper, & variation in lower (ie, 2 different variations), but I suppose I can't be sure. On a different note, is part of your question how to conduct MC simulation in R? – gung Apr 15 '12 at 16:28
• I came to the conclusion that it was looking for 2 different variations aswell, so think I will try that! – Sarah Apr 15 '12 at 16:32
• Fair enough, what I'm wondering now is, are you good to go at this point, or do you need extra hints about how to approach what we've guessed the question is asking for? – gung Apr 15 '12 at 17:27
• So, on each iteration, are you storing the upper & lower limits (hint)? – gung Apr 15 '12 at 18:34

My understanding of this homework assignment is that it is asking you to do a Monte Carlo simulation for a certain statistical technique. That is, you simulate a bunch of random data sets, apply this technique to these data set, and store the numbers to summarize later in a convenient way (means, simulated probabilities, etc.)

Now, the technique in question is the bootstrap, which involves Monte Carlo simulation within it (unless, as whuber demonstrated, you are asked to perform the exact bootstrap, which is highly unlikely). So as outcomes of your simulations, you may be storing the sample mean, the sample standard deviation, and the limits of confidence interval for the mean obtained by the bootstrap.

I think to make this study interesting, it is worth taking a small size sample ($n=10$) from a skewed distribution (exponential). That way, the standard CLT-based confidence interval will likely have inadequate performance: say, for 90% confidence interval that would need to miss 5% at either end, you will see something like 10% and 2%. The bootstrap CI may show better performance (say, 6% and 4%).

• One thing I have learned from working with the exact bootstrap is that the standard normal-theory CI will often work much better with small samples even when they are skewed. The bootstrap comes into its own with large samples: the theory is an asymptotic one, not an exact one. – whuber Apr 17 '12 at 21:50
• @whuber, that's an interesting observation. I usually stress the asymptotic nature of the bootstrap (see Peter Hall's The Bootstrap and Edgeworth Expansions book to understand why), but failed to do it this time. – StasK Apr 19 '12 at 18:36