# Monte Carlo study to estimate bootstrap CI

I am asked to conduct a Monte Carlo study to estimate the coverage probabilities of the standard normal bootstrap confidence interval, the basic bootstrap confidence interval, and the percentile confidence interval.

I am using the following code:

------------  START ---------------------

B    <- 100                 # number of replicates
mu   <- 1                   # for generating data: true mean
sd   <- 1                   # for generating data: true sd
Data <- rnorm(N, mu, sd)    # simulated data: original sample

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(), R=50 bootstrap replicates
bOot = boot(Data, statistic=getM, R=50)
boot.ci(bOot, conf=0.95, type=c("basic", "perc"))

---------------- END -----------------


I need to "wrap" this into a function that generates a vector with bootstrap confidence intervals so that I can further determine the proportion of times that the confidence intervals miss on the left, and the proportion of times that the confidence intervals miss on the right.

I am a novice to R programming and played around with some "for loops" without success. Can you direct me to where I need to go.

You may get better answers on stackoverflow, as this question relates more to programming than statistics.

It's not completely clear what you would like the code to produce, since a vector is 1 dimensional, so you could hold 1 item in the vector for each bootstrap confidence interval - but the information you want includes more than 1 thing - the type, the confidence level, the upper bound and the lower bound etc. So one simple way to do this is to use a list and enclose your code in a for loop, and add the object (which is in fact another list) returned by boot.ci like this:

require(boot)                   # for boot(), boot.ci()

N <-100
B    <- 100                 # number of replicates
mu   <- 1                   # for generating data: true mean
sd   <- 1                   # for generating data: true sd

getM <- function(orgDV, idx)
{bsM  <- mean(orgDV[idx])                       # M*
bsS2M <- (((N-1) / N) * var(orgDV[idx])) / N    # S^2*(M)
c(bsM, bsS2M)
}

N.CI <- 10             # number of botstrap confidence intervals
list.CI <- list(N.CI)  # a list to hold the objects returned by boot.ci

for (i in 1:N.CI ) {

Data <- rnorm(N, mu, sd)    # simulated data: original sample
bOot = boot(Data, statistic=getM, R=50)
list.CI[[i]] <- boot.ci(bOot, conf=0.95, type=c("basic", "perc"))
}


list.CI now contains 10 boostrap confidence intervals. You can access each one by list.CI[[1]], list.CI[[2]] etc, and then you can access the components of each CI like this:

> list.CI[[i]]$perc conf [1,] 0.95 1.28 49.73 0.7794841 1.318904 > list.CI[[i]]$basic
conf
[1,] 0.95 49.73 1.28 0.7980327 1.337453


Perhaps a better way, which will help your subsequent analysis is to store the bootstrap CIs in a matrix or a dataframe instead of a list, but I'll leave that for you to do.