# S-PLUS code for bootstrap confidence intervals # for medians: Problem 4 on assignment 9 in 2002 # Problem 4.(a) # Enter the data heart1 <- scan("heart1.dat") heart1 # Set the number of bootstrap samples nboot <- 5000 # Compute the sample median median(heart1) # Generate the bootstrap samples m.boot1 <- bootstrap(heart1, statistic=median, B=nboot) summary(m.boot1) # Draw a histogram of m.boot1. First # find the extremes of the bootstrap # samples. Unix users may use the # truehist plotting function library(MASS) mm <- range(m.boot1\$rep) min.int <- floor(mm[1]) max.int <- ceiling(mm[2]) truehist(m.boot1\$rep,xlim=c(min.int,max.int),density=.001) abline(v=median(heart1),lty=2) title("5000 bootstrap samples: median(heart1.dat)") # Problem 4.(b) # Enter the data heart2 <- scan("heart2.dat") # Compute the sample median median(heart2) # Generate the bootstrap samples m.boot2 <- bootstrap(heart2, statistic=median, B=nboot) summary(m.boot2) # Draw a histogram of m.boot2. First # find the extremes of the bootstrap # samples. Unix users may use the # truehist plotting function library(MASS) mm <- range(m.boot2\$rep) min.int <- floor(mm[1]) max.int <- ceiling(mm[2]) truehist(m.boot2\$rep,xlim=c(min.int,max.int),density=.001) abline(v=median(heart2),lty=2) title("5000 bootstrap samples: median(heart2.dat)") # Problem 4.(c) # construct 95% percentile confidence limits # for the difference in population medians # Compute differences in bootstrapped replicates diff.boot <- m.boot1\$rep- m.boot2\$rep diff.boot <- sort(diff.boot) # Compute upper and lower limits kl<-floor((nboot+1)*.025) ku<-nboot+1-kl diff.lower <- round(diff.boot[kl],digits=4) diff.upper <- round(diff.boot[ku],digits=4) c("Lower Limit"=diff.lower, "Upper Limit"=diff.upper) # Create a histogram of differences mm <- range(diff.boot) min.int <- floor(mm[1]) max.int <- ceiling(mm[2]) truehist(diff.boot,xlim=c(min.int,max.int),density=.001) title("Differences of 5000 bootstrapped medians")