#374 Resampling Slide 43

# bootstrap comparison of medians

B<-10000
treatmentsample <- c(94, 197, 16, 38, 99, 141, 23)
controlsample <- c(52,104,146,10,50,31,40,27,46)

brep <- rep(0,B)
for(i in 1:B){                             
  btsample <- sample(treatmentsample, length(treatmentsample), replace=T)
  bcsample <- sample(controlsample, length(controlsample), replace=T)
  brep[i] <- median(btsample) - median(bcsample)
}

hist(brep, breaks=(-150:200), freq=FALSE, col="red",
        main="Bootstrap distribution for difference between medians",
        xlab="Difference between medians", ylab="Probability")
cat('Difference between medians', 
    median(treatmentsample) - median(controlsample), '\n')
cat('bootstrap estimate of std dev', sd(brep), '\n')
 

# bootstrap comparison of means

brep <- rep(0,B)
for(i in 1:B){                             
  btsample <- sample(treatmentsample, length(treatmentsample), replace=T)
  bcsample <- sample(controlsample, length(controlsample), replace=T)
  brep[i] <- mean(btsample) - mean(bcsample)
}

hist(brep, breaks=(-150:200), freq=FALSE, col="red",
        main="Bootstrap distribution for difference between means",
        xlab="Difference between means", ylab="Probability")
cat('Difference between means', 
    mean(treatmentsample) - mean(controlsample), '\n')
cat('bootstrap estimate of std dev', sd(brep), '\n')
cat('algebraic estimate of std dev', sqrt(var(treatmentsample)/length(treatmentsample) + 
                                          var(controlsample)/length(controlsample)), '\n')