library(bootstrap)

# part 1
# how much variation in our bootstrap estimate is due to B?

theta <- function(ind) cor(law[ind,1], law[ind,2])

for (B in c(25, 100, 400)) {
  se_boot <- rep(0, 100)
  for (n in 1:100) {
    se_boot[n] <- sd(bootstrap(1:15, B, theta)$thetastar)
  }
  cat('B =', B, 'mean se_boot =', mean(se_boot), 'sd of se_boot =', sd(se_boot), '\n') 
}

# part 2
# ideally instead of bootstrap samples we would like to use repeated samples
# from the original population

theta_82 <- function(ind) cor(law82[ind,2], law82[ind,3])

for (B in c(25, 100, 400)) {
  se_repeat <- rep(0, 100)
  for (n in 1:100) {
    cor_repeat <- rep(0, B)
    for (i in 1:B) {
      newsample <- sample(1:82, 15)
      cor_repeat[i] <- theta_82(newsample)
    }
    se_repeat[n] <- sd(cor_repeat)
  }
  cat('B =', B, 'mean se_repeat =', mean(se_repeat), 'sd of se_repeat =', sd(se_repeat), '\n')
}