regresboot<-function(x,y,conflevel=0.95,nboots=120){
#This function estimates the CI for the slope in #
#linear regression using resampling from the residuals.#
alpha<-1-conflevel
b1hatstar<-rep(0,nboots)
model<-lm(y~x)
b0hat<-model$coefficients[1]
b1hat<-model$coefficients[2]
resids<-model$residuals
for (b in 1:nboots){
  ystar<-b0hat+b1hat*x+sample(resids,replace=T)
  b1hatstar[b]<-lm(ystar~x)$coefficients[2]
}
tstar<-sort(b1hatstar)
lowtstar<-0.5*(tstar[floor(nboots*alpha/2)]+tstar[ceiling(nboots*alpha/2)])
hitstar<-0.5*(tstar[floor(nboots*(1-alpha/2))]+tstar[ceiling(nboots*(1-alpha/2))])
c(lowtstar,hitstar)
}


regptsboot<-function(x,y,conflevel=0.95,nboots=120){
#This function estimates the CI for the slope in #
#linear regression using resampling of the observations.#
alpha<-1-conflevel
b1hatstar<-rep(0,nboots)
for (b in 1:nboots){
  pointstar<-sample(1:length(x),replace=T)
  b1hatstar[b]<-lm(y[pointstar]~x[pointstar])$coefficients[2]
}
tstar<-sort(b1hatstar)
lowtstar<-0.5*(tstar[floor(nboots*alpha/2)]+tstar[ceiling(nboots*alpha/2)])
hitstar<-0.5*(tstar[floor(nboots*(1-alpha/2))]+tstar[ceiling(nboots*(1-alpha/2))])
c(lowtstar,hitstar)
}


success<-matrix(0,100,2)
ci<-array(0,dim=c(100,2,2))
for (i in 1:100){
  x<-1:20
  y<-5+3*x+rnorm(20,0,0.5)
 ci[i,,1]<-regresboot(x,y,nboots=120)
 ci[i,,2]<-regptsboot(x,y,nboots=120)
 success[i,1]<-((ci[i,1,1]<3)&(ci[i,2,1]>3))
 success[i,2]<-((ci[i,1,2]<3)&(ci[i,2,2]>3))
}
apply(success,2,mean)
apply(ci[,2,]-ci[,1,],2,mean)