#
#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

# fix the random.seed, so I'll always get the same answer:
set.seed(10131985)

# The t-based confidence intervals:
# 
x = seq(-3,3,le=5)
y = 2 + 4 * x + rnorm(5)
fit = lm(y~x)
fitsum = summary(fit)

fitsum$coefficients

# extract the simple linear regression beta_i estimates and their standard error:
#
beta0Hat = fitsum$coefficients[1,1]
beta0SE  = fitsum$coefficients[1,2]
beta1Hat = fitsum$coefficients[2,1]
beta1SE  = fitsum$coefficients[2,2]

alpha = 0.1
n = length(x)
c = qt( 1-0.5*alpha, n-2 ) # the critical value for these tests

# Compute 1-alpha = 0.95 confidence intervals on alpha_0, alpha_1
print( sprintf("the t-based %.2f%% confidence interval beta_0 is (%10.6f,%10.6f)",1-alpha,beta0Hat - c * beta0SE, beta0Hat + c * beta0SE) )
print( sprintf("the t-based %.2f%% confidence interval beta_1 is (%10.6f,%10.6f)",1-alpha,beta1Hat - c * beta1SE, beta1Hat + c * beta1SE) )



# The residual bootstrapped confidence intervals on beta_0 and beta_1
# 
Rdata = fit$residuals
nBoot = 2000
B=array(0,dim=c(nBoot,2))
for(i in 1:nBoot){
  ystar = y + sample(Rdata,replace=T) # add errors to response 
  Bfit = lm(ystar~x)
  B[i,] = Bfit$coefficients
}

# compute the 90% confidence interval from the bootstrap:
#
quantile( B[,1], c(0.05,0.95) ) # for beta_0 
quantile( B[,2], c(0.05,0.95) ) # for beta_1