Skip to content

Instantly share code, notes, and snippets.

Created August 3, 2015 07:15
Show Gist options
  • Star 3 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save anonymous/73f64e1b8ff7e972fc3b to your computer and use it in GitHub Desktop.
Save anonymous/73f64e1b8ff7e972fc3b to your computer and use it in GitHub Desktop.
getCI <- function(B, muH0, sdH0, N) {
getM <- function(orgDV, idx) {
bsM <- mean(orgDV[idx]) # M*
bsS2M <- (((N-1) / N) * var(orgDV[idx])) / N # S^2*(M)
c(bsM, bsS2M)
}
DV <- rnorm(N, muH0, sdH0) # simulated data: original sample
M <- mean(DV) # M from original sample
S2M <- (((N-1)/N) * var(DV)) / N # S^2(M) from original sample
# bootstrap
boots <- t(replicate(B, getM(DV, sample(seq(along=DV), replace=TRUE))))
Mstar <- boots[ , 1] # M* for each replicate
S2Mstar <- boots[ , 2] # S^2*(M) for each replicate
biasM <- mean(Mstar) - M # bias of estimator M
# indices for sorted vector of estimates
idx <- trunc((B + 1) * c(0.05/2, 1 - 0.05/2))
zCrit <- qnorm(c(1 - 0.05/2, 0.05/2)) # z-quantiles from std-normal distribution
tStar <- (Mstar-M) / sqrt(S2Mstar) # t*
tCrit <- sort(tStar)[rev(idx)] # t-quantiles from empirical t* distribution
ciBasic <- 2*M - sort(Mstar)[rev(idx)] # basic CI
ciPerc <- sort(Mstar)[idx] # percentile CI
ciNorm <- M-biasM - zCrit*sd(Mstar) # normal CI
ciT <- M - tCrit * sqrt(S2M) # studentized t-CI
c(basic=ciBasic, percentile=ciPerc, normal=ciNorm, t=ciT)
}
## 1000 bootstraps - this will take a while
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
Nrep <- 1000 # number of bootstraps
CIs <- t(replicate(Nrep, getCI(B=B, muH0=muH0, sdH0=sdH0, N=N)))
## coverage probabilities
sum((CIs[ , "basic1"] < muH0) & (CIs[ , "basic2"] > muH0)) / Nrep
sum((CIs[ , "percentile1"] < muH0) & (CIs[ , "percentile2"] > muH0)) / Nrep
sum((CIs[ , "normal1"] < muH0) & (CIs[ , "normal2"] > muH0)) / Nrep
sum((CIs[ , "t1"] < muH0) & (CIs[ , "t2"] > muH0)) / Nrep
@buckeye-netizen
Copy link

I just found this to be a very helpful function for my own research. So thank you! Just a quick question: what is the difference between B and Nrep? Is B just Nrep minus 1?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment