# This program uses the `bootstrap' package to bootstrap
# the correlation coefficient for the Law data.

library(bootstrap)    # loads library

B<-3200               # Sets replications

par(mfrow=c(2,2))     # Sets up multiple plots 2x2
plot(law, main="Law data - Original sample of 15",
          xlim=c(450,750), xlab="LSAT", ylab="GPA")

# With bivariate data we define the basic data set as the case indices
# ie 1:15 in this example.  The bootstrap command will then
# generate bootstrap samples from 1:15 and input these samples
# to the statistic.  So the function theta defining the
# statistic is written as a function of `ind' - the case indices.

theta <- function(ind) cor(law[ind,1], law[ind,2])
cat('sample estimate of correlation', theta(1:15), '\n') # sample estimate
law.boot <- bootstrap(1:15, B, theta)
cat('bootstrap est of sd(theta)', sd(law.boot$thetastar), '\n') # bootstrap standard error
replicates <- law.boot$thetastar[is.finite(law.boot$thetastar)]
hist(replicates,xlim=c(-1.05,1.05),breaks=seq(from=-1.05,to=1.05,by=0.05))

# Use Format>Block>Comment/Uncomment to activate this section
# which compares the bootstrap distribution for
# the correlation coefficient with the
# distribution for repeated samples of size 15
# from the population of 82 values.

plot(law82[,2], law82[,3], main="Law data - Population size 82",
                           xlim=c(450,750), xlab="LSAT", ylab="GPA")

poprep<-c()
theta_82 <- function(ind) cor(law82[ind,2], law82[ind,3])
cat('population correlation', theta_82(1:82), '\n') #`population' correlation
for(i in 1:B){
  popsample <- sample(1:82,15,replace=T)
  poprep[i] <- theta_82(popsample)
}
cat('monte carlo est of sd(theta)', sd(poprep), '\n') # monte carlo standard error
hist(poprep[is.finite(poprep)], xlim=c(-1.05,1.05),
                                breaks=seq(-1.05,to=1.05,by=0.05))