library(bootstrap) # for law data
library(boot) # for freq.array

B <- 1000               # Sets replications

theta <- function(ind) cor(law[ind,1], law[ind,2])
cat('ratio estimate', theta(1:15), '\n')
boot.samples <- matrix(0, nrow=B, ncol=15)
boot.replicates <- rep(0, B)
for (i in 1:B) {
  boot.samples[i,] <- sample(1:15, 15, T)
  boot.replicates[i] <- theta(boot.samples[i,])
}
cat('bootstrap est of sd', sd(boot.replicates), '\n')
cat('bootstrap est of bias', mean(boot.replicates) - theta(1:15), '\n')

theta_f <- function(f) {
  xbar <- sum(law[,1]*f)
  ybar <- sum(law[,2]*f)
  xycov <- sum((law[,1] - xbar)*(law[,2] - ybar)*f)
  xvar <- sum((law[,1] - xbar)^2*f)
  yvar <- sum((law[,2] - ybar)^2*f)
  xycor <- xycov/sqrt(xvar*yvar)
  return(xycor)
}
cat('ratio estimate', theta_f(rep(1/15, 15)), '\n')
boot.samples_f <- freq.array(boot.samples)/15
boot.replicates_f <- rep(0, B)
for (i in 1:B) {
  boot.replicates_f[i] <- theta_f(boot.samples_f[i,])
}
cat('bootstrap est of sd', sd(boot.replicates_f), '\n')
cat('bootstrap est of bias', mean(boot.replicates_f) - theta_f(rep(1/15, 15)), '\n')

cat('controlled est of bias', 
     mean(boot.replicates_f) - theta_f(apply(boot.samples_f, 2, mean)), '\n')