llawboot.samp<-matrix(0,ncol=1000,nrow=15)
for (i in 1:1000) {
  llawboot.samp[,i]<-sample(law[,1],15,replace=T) }

llawboot.mean<-apply(llawboot.samp,2,mean)
llawboot.med<-apply(llawboot.samp,2,median)
llawboot.stddev<-sqrt(apply(llawboot.samp,2,var))

llbmean.se <- rep(0,20)
llbmed.se <- rep(0,20)
llbsd.se <- rep(0,20)

for(i in 50:50:1000) {
  llbmean.se[i/50] <-sqrt(var(llawboot.mean[1:i]))
  llbmed.se[i/50] <- sqrt(var(llawboot.med[1:i]))
  llbsd.se[i/50] <- sqrt(var(llawboot.stddev[1:i]))
}

llawsamp.samp<-matrix(0,ncol=1000,nrow=15)
for (i in 1:1000) {
  llawsamp.samp[,i]<-sample(law82[,2],15,replace=T) }

llawsamp.mean<-apply(llawsamp.samp,2,mean)
llawsamp.med<-apply(llawsamp.samp,2,median)
llawsamp.stddev<-sqrt(apply(llawsamp.samp,2,var))

llsmean.se <- rep(0,20)
llsmed.se <- rep(0,20)
llssd.se <- rep(0,20)

for(i in 50:50:1000) {
  llsmean.se[i/50] <-sqrt(var(llawsamp.mean[1:i]))
  llsmed.se[i/50] <- sqrt(var(llawsamp.med[1:i]))
  llssd.se[i/50] <- sqrt(var(llawsamp.stddev[1:i]))
}

par(mfrow=c(1,2))
hist(llawboot.mean,prob=T,main="LSAT - Bootstrap - Mean",xlab="Mean",xlim=c(560,640),ylim=c(0,0.04))
hist(llawsamp.mean,prob=T,main="LSAT - Sample - Mean",xlab="Mean",xlim=c(560,640),ylim=c(0,0.04))

hist(llawboot.med,prob=T,main="LSAT - Bootstrap - Median",xlab="Median",xlim=c(540,660),ylim=c(0,0.06))
hist(llawsamp.med,prob=T,main="LSAT - Sample - Median",xlab="Median",xlim=c(540,660),ylim=c(0,0.06))

hist(llawboot.stddev,prob=T,main="LSAT - Bootstrap - Std Deviation",xlab="Std Dev",xlim=c(10,65),ylim=c(0,0.08))
hist(llawsamp.stddev,prob=T,main="LSAT - Sample - Std Deviation",xlab="Std Dev",xlim=c(10,65),ylim=c(0,0.08))

llbmean.bias <- rep(0,20)
llbmed.bias <- rep(0,20)
llbsd.bias <- rep(0,20)

for(i in 50:50:1000) {
  llbmean.bias[i/50] <- mean(llawboot.mean[1:i])-mean(law[,1])
  llbmed.bias[i/50] <- mean(llawboot.med[1:i])-median(law[,1])
  llbsd.bias[i/50] <- mean(llawboot.stddev[1:i])-sqrt(var(law[,1]))
}

llsmean.bias <- rep(0,20)
llsmed.bias <- rep(0,20)
llssd.bias <- rep(0,20)

for(i in 50:50:1000) {
  llsmean.bias[i/50] <- mean(llawsamp.mean[1:i])-mean(law82[,2])
  llsmed.bias[i/50] <- mean(llawsamp.med[1:i])-median(law82[,2])
  llssd.bias[i/50] <- mean(llawsamp.stddev[1:i])-sqrt(var(law82[,2])*81/82)
}

par(mfrow=c(1,2))
hist(llawboot.mean,prob=T,main="LSAT - Bootstrap - Mean",xlab="Mean",xlim=c(560,640),ylim=c(0,0.04))
points(mean(law[,1]),0)
hist(llawsamp.mean,prob=T,main="LSAT - Sample - Mean",xlab="Mean",xlim=c(560,640),ylim=c(0,0.04))
points(mean(law82[,2]),0)

hist(llawboot.med,prob=T,main="LSAT - Bootstrap - Median",xlab="Median",xlim=c(540,660),ylim=c(0,0.06))
points(median(law[,1]),0)
hist(llawsamp.med,prob=T,main="LSAT - Sample - Median",xlab="Median",xlim=c(540,660),ylim=c(0,0.06))
points(median(law82[,2]),0)

hist(llawboot.stddev,prob=T,main="LSAT - Bootstrap - Std Deviation",xlab="Std Dev",xlim=c(10,65),ylim=c(0,0.08))
points(sqrt(var(law[,1])),0)
hist(llawsamp.stddev,prob=T,main="LSAT - Sample - Std Deviation",xlab="Std Dev",xlim=c(10,65),ylim=c(0,0.08))
points(sqrt(var(law82[,2])*81/82),0)

glawboot.samp<-matrix(0,ncol=1000,nrow=15)
for (i in 1:1000) {
  glawboot.samp[,i]<-sample(law[,2],15,replace=T) }

glawboot.mean<-apply(glawboot.samp,2,mean)
glawboot.med<-apply(glawboot.samp,2,median)
glawboot.stddev<-sqrt(apply(glawboot.samp,2,var))

glbmean.se <- rep(0,20)
glbmed.se <- rep(0,20)
glbsd.se <- rep(0,20)

for(i in 50:50:1000) {
  glbmean.se[i/50] <-sqrt(var(glawboot.mean[1:i]))
  glbmed.se[i/50] <- sqrt(var(glawboot.med[1:i]))
  glbsd.se[i/50] <- sqrt(var(glawboot.stddev[1:i]))
}

glawsamp.samp<-matrix(0,ncol=1000,nrow=15)
for (i in 1:1000) {
  glawsamp.samp[,i]<-sample(law82[,3],15,replace=T) }

glawsamp.mean<-apply(glawsamp.samp,2,mean)
glawsamp.med<-apply(glawsamp.samp,2,median)
glawsamp.stddev<-sqrt(apply(glawsamp.samp,2,var))

glsmean.se <- rep(0,20)
glsmed.se <- rep(0,20)
glssd.se <- rep(0,20)

for(i in 50:50:1000) {
  glsmean.se[i/50] <-sqrt(var(glawsamp.mean[1:i]))
  glsmed.se[i/50] <- sqrt(var(glawsamp.med[1:i]))
  glssd.se[i/50] <- sqrt(var(glawsamp.stddev[1:i]))
}

par(mfrow=c(1,2))
hist(glawboot.mean,prob=T,main="GPA - Bootstrap - Mean",xlab="Mean",xlim=c(2.85,3.35),ylim=c(0,8))
hist(glawsamp.mean,prob=T,main="GPA - Sample - Mean",xlab="Mean",xlim=c(2.85,3.35),ylim=c(0,8))

hist(glawboot.med,prob=T,main="GPA - Bootstrap - Median",xlab="Median",xlim=c(2.8,3.45),ylim=c(0,7))
hist(glawsamp.med,prob=T,main="GPA - Sample - Median",xlab="Median",xlim=c(2.8,3.45),ylim=c(0,7))

hist(glawboot.stddev,prob=T,main="GPA - Bootstrap - Std Deviation",xlab="Std Dev",xlim=c(0.08,0.32),ylim=c(0,14))
hist(glawsamp.stddev,prob=T,main="GPA - Sample - Std Deviation",xlab="Std Dev",xlim=c(0.08,0.32),ylim=c(0,14))

glbmean.bias <- rep(0,20)
glbmed.bias <- rep(0,20)
glbsd.bias <- rep(0,20)

for(i in 50:50:1000) {
  glbmean.bias[i/50] <- mean(glawboot.mean[1:i])-mean(law[,2])
  glbmed.bias[i/50] <- mean(glawboot.med[1:i])-median(law[,2])
  glbsd.bias[i/50] <- mean(glawboot.stddev[1:i])-sqrt(var(law[,2]))
}

glsmean.bias <- rep(0,20)
glsmed.bias <- rep(0,20)
glssd.bias <- rep(0,20)

for(i in 50:50:1000) {
  glsmean.bias[i/50] <- mean(glawsamp.mean[1:i])-mean(law82[,3])
  glsmed.bias[i/50] <- mean(glawsamp.med[1:i])-median(law82[,3])
  glssd.bias[i/50] <- mean(glawsamp.stddev[1:i])-sqrt(var(law82[,3])*81/82)
}

par(mfrow=c(1,2))
hist(glawboot.mean,prob=T,main="GPA - Bootstrap - Mean",xlab="Mean",xlim=c(2.85,3.35),ylim=c(0,8))
points(mean(law[,2]),0)
hist(glawsamp.mean,prob=T,main="GPA - Sample - Mean",xlab="Mean",xlim=c(2.85,3.35),ylim=c(0,8))
points(mean(law82[,3]),0)

hist(glawboot.med,prob=T,main="GPA - Bootstrap - Median",xlab="Median",xlim=c(2.8,3.45),ylim=c(0,7))
points(median(law[,2]),0)
hist(glawsamp.med,prob=T,main="GPA - Sample - Median",xlab="Median",xlim=c(2.8,3.45),ylim=c(0,7))
points(median(law82[,3]),0)

hist(glawboot.stddev,prob=T,main="GPA - Bootstrap - Std Deviation",xlab="Std Dev",xlim=c(0.08,0.32),ylim=c(0,14))
points(sqrt(var(law[,2])),0)
hist(glawsamp.stddev,prob=T,main="GPA - Sample - Std Deviation",xlab="Std Dev",xlim=c(0.08,0.32),ylim=c(0,14))
points(sqrt(var(law82[,3])*81/82),0)

lawbootr.index<-matrix(0,ncol=1000,nrow=15)
for (i in 1:1000)
  lawbootr.index[,i] <- sample(1:15,15,replace=T)

lawbootr.samp<-array(0,c(15,2,1000))
for (i in 1:1000) {
  lawbootr.samp[,,i] <- as.matrix(law[lawbootr.index[,i],]) }

lawbootr.corr <- rep(0,1000)
for (i in 1:1000) {
  lawbootr.corr[i] <- cor(lawbootr.samp[,1,i],lawbootr.samp[,2,i]) }

lbr.se <- rep(0,20)
for (i in 50:50:1000)
  lbr.se[i/50] <- sqrt(var(lawbootr.corr[1:i]))

lbr.bias <- rep(0,20)
for (i in 50:50:1000)
  lbr.bias[i/50] <- mean(lawbootr.corr[1:i]) - cor(law[,1],law[,2])

lawsampr.index<-matrix(0,ncol=1000,nrow=15)
for (i in 1:1000)
  lawsampr.index[,i] <- sample(1:82,15,replace=F)

lawsampr.samp<-array(0,c(15,2,1000))
for (i in 1:1000) {
  lawsampr.samp[,,i] <- as.matrix(law82[lawsampr.index[,i],2:3]) }

lawsampr.corr <- rep(0,1000)
for (i in 1:1000) {
  lawsampr.corr[i] <- cor(lawsampr.samp[,1,i],lawsampr.samp[,2,i]) }

lsr.se <- rep(0,20)
for (i in 50:50:1000)
  lsr.se[i/50] <- sqrt(var(lawsampr.corr[1:i]))

lsr.bias <- rep(0,20)
for (i in 50:50:1000)
  lsr.bias[i/50] <- mean(lawsampr.corr[1:i]) - cor(law82[,2],law82[,3])


par(mfrow=c(1,2))
hist(lawbootr.corr,prob=T,br=(-2:20)/20,main="Bootstrap - Correlation",xlab="Corr",xlim=c(-0.1,1),ylim=c(0,3.5))
hist(lawsampr.corr,prob=T,br=(-2:20)/20,main="Sample - Correlation",xlab="Corr",xlim=c(-0.1,1),ylim=c(0,3.5))

par(mfrow=c(1,2))
hist(lawbootr.corr,prob=T,br=(-2:20)/20,main="Bootstrap - Correlation",xlab="Corr",xlim=c(-0.1,1),ylim=c(0,3.5))
points(cor(law[,1],law[,2]),0)
hist(lawsampr.corr,prob=T,br=(-2:20)/20,main="Sample - Correlation",xlab="Corr",xlim=c(-0.1,1),ylim=c(0,3.5))
points(cor(law82[,2],law82[,3]),0)

par(mfrow=c(1,2))
hist(law82[,2],main="LSAT",xlab="LSAT")
points(law[,1],rep(0,15))
hist(law82[,3],main="GPA",xlab="GPA")
points(law[,2],rep(0,15))

par(mfrow=c(1,1))
plot(law82[,2],law82[,3],xlab="LSAT",ylab="GPA",col=gray(0.7))
points(law[,1],law[,2])
legend(650,2.7,c("Sampled","Unsampled"),col=gray(c(0,0.7)),pch=19)