# mouse : bootstrap mean and se of the mean # Efron and Tibshirani, 1994, An Introduction to the Bootstrap, chapter 2 # ============================================= treatment <- c(94, 197, 16, 38, 99, 141, 23) # ============================================= # bootstrap estimates of standard error of the mean of treatment group, mouse data # Chihara-Hesterberg style treatment <- c(94, 197, 16, 38, 99, 141, 23) n <- length(treatment) N <- 50 mean.hat.star <- numeric(N) for (i in 1:N){ samp <- sample(1:n, n, replace=TRUE) mean.hat.star[i] <- mean(treatment[samp]) } hist(mean.hat.star, col="dark red") sd(mean.hat.star) # => 22.54373 Ns <- c(50, 100, 250, 500, 1000) sd.mean.hat.star <- c(26.67048, 24.3416, 23.60291, 23.2366, 23.29813) # ======================================== # bootstrap - Efron-Tibshirani style bootstrap <- function(x, nboot, theta, ...){ data <- matrix(sample(x, size=length(x) * nboot, replace=TRUE), nrow=nboot) return(apply(data, 1, theta)) } # arguments # x a vector containing the data # nboot the number of bootstrap samples desired # theta function to be bootstrapped; takes x as an argument and may take additional arguments # ======================================== # Efron-Tibshirani style treatment <- c(94, 197, 16, 38, 99, 141, 23) n <- length(treatment) N <- 50 nboot <- N theta <- function(x){ mean(treatment[x]) } theta.hat.star <- bootstrap(1:n, nboot, theta) hist(theta.hat.star, col="dark red") mean(theta.hat.star) # => 88.12 sd(theta.hat.star) # => 23.41702 alpha <- .05 quantile(theta.hat.star, c(alpha / 2, 1 - alpha / 2)) # => 2.5% 97.5% # => 45.01786 134.56071 Ns <- c(50, 100, 250, 500, 1000) sd.mean.hat.star <- c(26.67048, 24.3416, 23.60291, 23.2366, 23.29813) mean(bootstrap(1:n, 50, theta)) mean(bootstrap(1:n, 100, theta)) mean(bootstrap(1:n, 250, theta)) mean(bootstrap(1:n, 500, theta)) mean(bootstrap(1:n, 1000, theta)) # ======================================== # mean.star : apply bootstrap # abstract bootstrap over N == nboot bootstrapET <- function(N){ mean(bootstrap(1:n, N, theta)) } bootstrapET(50) # and apply it to a list of Ns # => list of mean.stars for those Ns Ns <- c(50, 100, 250, 500, 1000) apply(cbind(Ns), 1, bootstrapET) # => 86.21714 85.52429 85.17086 87.15229 86.89257 # ======================================== # se.star : apply bootstrap # abstract bootstrap over N == nboot bootstrapET <- function(N){ sd(bootstrap(1:n, N, theta)) } bootstrapET(50) # and apply it to a list of Ns # => list of se.stars for those Ns Ns <- c(50, 100, 250, 500, 1000) apply(cbind(Ns), 1, bootstrapET) # => 24.03611 22.28045 23.98496 23.18574 23.37262 # ========================================