# Coverage probabilities of the basic bootstrap confidence Interval

I have the following question for a course I'm working on:

Conduct a Monte Carlo study to estimate the coverage probabilities of the standard normal bootstrap confidence interval and the basic bootstrap confidence interval. Sample from a normal population and check the empirical coverage rates for the sample mean.

Coverage probabilities for the standard normal bootstrap CI are easy:

n = 1000;
alpha = c(0.025, 0.975);
x = rnorm(n, 0, 1);
mu = mean(x);
sqrt.n = sqrt(n);

LNorm = numeric(B);
UNorm = numeric(B);

for(j in 1:B)
{
smpl = x[sample(1:n, size = n, replace = TRUE)];
xbar = mean(smpl);
s = sd(smpl);

LNorm[j] = xbar + qnorm(alpha) * (s / sqrt.n);
UNorm[j] = xbar + qnorm(alpha) * (s / sqrt.n);
}

mean(LNorm < 0 & UNorm > 0); # Approximates to 0.95
# NOTE: it is not good enough to look at overall coverage
# Must compute separately for each tail


From what I've been taught for this course, the basic bootstrap confidence interval can be calculated like this:

# Using x from previous...
R = boot(data = x, R=1000, statistic = function(x, i){ mean(x[i]); });
result = 2 * mu - quantile(R$t, alpha, type=1);  That makes sense. What I don't understand is how to calculate coverage probabilities for the basic bootstrap CI. I understand that the coverage probability would represent the number of times that the CI contains the true value (in this case mu). Do I simply run the boot function many times? How can I approach this question differently? • Is your size=100 a typo? I don't believe you're getting the right upper and lower bounds since the implicit sample size appears to be 1000 when you calculate your CIs in the loop (since you use sqrt.n in the calculation). Also, why are you comparing to mu and not 0 directly (the latter being the true mean)? Apr 3, 2011 at 19:24 • Also, smpl = x[sample(1:n, size = 100, replace = TRUE)];  can be simplified to smpl = sample(x, size=100, replace=TRUE). Apr 3, 2011 at 19:25 • @cardinal - Yes it was a typo and same with mu being 0. The normal CI works fine, it's the basic bootstrap CI that I've having difficulty with. Apr 3, 2011 at 21:47 ## 1 Answer The terminology is probably not used consistently, so the following is only how I understand the original question. From my understanding, the normal CIs you computed are not what was asked for. Each set of bootstrap replicates gives you one confidence interval, not many. The way to compute different CI-types from the results of a set of bootstrap replicates is as follows: B <- 999 # number of replicates muH0 <- 100 # for generating data: true mean sdH0 <- 40 # for generating data: true sd N <- 200 # sample size DV <- rnorm(N, muH0, sdH0) # simulated data: original sample  Since I want to compare the calculations against the results from package boot, I first define a function that will be called for each replicate. Its arguments are the original sample, and an index vector specifying the cases for a single replicate. It returns$M^{\star}$, the plug-in estimate for$\mu$, as well as$S_{M}^{2\star}$, the plug-in estimate for the variance of the mean$\sigma_{M}^{2}$. The latter will be required only for the bootstrap$t$-CI. > getM <- function(orgDV, idx) { + bsM <- mean(orgDV[idx]) # M* + bsS2M <- (((N-1) / N) * var(orgDV[idx])) / N # S^2*(M) + c(bsM, bsS2M) + } > library(boot) # for boot(), boot.ci() > bOut <- boot(DV, statistic=getM, R=B) > boot.ci(bOut, conf=0.95, type=c("basic", "perc", "norm", "stud")) BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS Based on 999 bootstrap replicates CALL : boot.ci(boot.out = bOut, conf = 0.95, type = c("basic", "perc", "norm", "stud")) Intervals : Level Normal Basic Studentized Percentile 95% ( 95.6, 106.0 ) ( 95.7, 106.2 ) ( 95.4, 106.2 ) ( 95.4, 106.0 ) Calculations and Intervals on Original Scale  Without using package boot you can simply use replicate() to get a set of bootstrap replicates. boots <- t(replicate(B, getM(DV, sample(seq(along=DV), replace=TRUE))))  But let's stick with the results from boot.ci() to have a reference. boots <- bOut$t                     # estimates from all replicates
M       <- mean(DV)                   # M from original sample
S2M     <- (((N-1)/N) * var(DV)) / N  # S^2(M) from original sample
Mstar   <- boots[ , 1]                # M* for each replicate
S2Mstar <- boots[ , 2]                # S^2*(M) for each replicate
biasM   <- mean(Mstar) - M            # bias of estimator M


The basic, percentile, and $t$-CI rely on the empirical distribution of bootstrap estimates. To get the $\alpha/2$ and $1 - \alpha/2$ quantiles, we find the corresponding indices to the sorted vector of bootstrap estimates (note that boot.ci() will do a more complicated interpolation to find the empirical quantiles when the indices are not natural numbers).

(idx <- trunc((B + 1) * c(0.05/2, 1 - 0.05/2)) # indices for sorted vector of estimates
 25 975

> (ciBasic <- 2*M - sort(Mstar)[idx])          # basic CI
 106.21826  95.65911

> (ciPerc <- sort(Mstar)[idx])                 # percentile CI
 95.42188 105.98103


For the $t$-CI, we need the bootstrap $t^{\star}$ estimates to calculate the critical $t$-values. For the standard normal CI, the critical value will just be the $z$-value from the standard normal distribution.

# standard normal CI with bias correction
> zCrit   <- qnorm(c(0.025, 0.975))   # z-quantiles from std-normal distribution
> (ciNorm <- M - biasM + zCrit * sqrt(var(Mstar)))
 95.5566 106.0043

> tStar <- (Mstar-M) / sqrt(S2Mstar)  # t*
> tCrit <- sort(tStar)[idx]           # t-quantiles from empirical t* distribution
> (ciT  <- M - tCrit * sqrt(S2M))     # studentized t-CI
 106.20690  95.44878


In order to estimate the coverage probabilities of these CI-types, you will have to run this simulation many times. Just wrap the code into a function, return a list with the CI-results and run it with replicate() like demonstrated in this gist.

• Wow! - Awesome explanation on what I was doing wrong. Also - thanks for the code tips! This works perfectly! Apr 4, 2011 at 15:36
• Ok one last question: when I try to replicate this information, I created a function computeCIs and called results = replicate(500, computeCIs());. At the end of computeCIs it returns c(ciBasic, ciPerc). To test the coverage probabilities, shouldn't I then test for mean(results[1, ] < 0 & results[2, ] > 0) to test for all Basic CI's containing the true mean (the coverage probability)? When I run this, I get 1 when I think I should get 0.95. Apr 4, 2011 at 16:02
• @TheCloudlessSky For the complete function and full simulation with the expected results in terms of coverage frequencies, see pastebin.com/qKpNKK0D Apr 4, 2011 at 20:33
• Yup, I'm an idiot :)... I made a typo when copying the code in R... thanks for all your help! :) Apr 4, 2011 at 22:28
• Thanks @caracal for nice answer. The link pastebin.com/qKpNKK0D is broken. Would appreciate if you update it and provide the complete function and full simulation. Thanks Aug 1, 2015 at 5:24