# 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")