## To study confidence and prediction intervals for cars data

> plot(cars,xlim=c(0,25))

> X <- cars[,1]
> Y <- cars[,2]
> n <- length(X)
> Design <- cbind(X,X^2)

> Mod <- lm(Y~Design-1)                # Fits Y_i = beta x_i + gamma x_i^2 + epsilon_i
> beta.hat <- Mod$coef
> beta.hat

> x <- seq(0,30,length=151)
> lines(x,beta.hat[[1]]*x + beta.hat[[2]]*x^2)

> p <- 2
> sigma.tilde.sqd <- sum((Y - Design%*%beta.hat)^2)/(n-p)

## First let’s compute standard errors for the coefficient estimates

> Std.err1 <- sqrt(sigma.tilde.sqd*(solve(t(Design)%*%Design))[1,1])
> Std.err2 <- sqrt(sigma.tilde.sqd*(solve(t(Design)%*%Design))[2,2])
> Std.err1
> Std.err2

> summary(Mod)    # Compare the estimates and standard errors

## Now let's compute a confidence interval for the expected stopping distance at initial speed 21mph

> xstar <- c(21,21^2)
> tau.sqd <- t(xstar)%*%solve(t(Design)%*%Design)%*%xstar
> tau.sqd

> Est <- t(xstar)%*%beta.hat
> Est 

> Lower <- Est - sqrt(sigma.tilde.sqd)*sqrt(tau.sqd)*qt(0.975,n-p)    # qt(0.975,n-p) is the upper 0.025 
                                                                      # point of the t_{n-p} distribution

> Lower
> Upper <- Est + sqrt(sigma.tilde.sqd)*sqrt(tau.sqd)*qt(0.975,n-p)
> Upper

> segments(xstar[1],Lower,xstar[1],Upper)             # Draws a straight line from 
                                                      # (xstar[1],Lower) to 
                                                      # (xstar[1],Upper)


## Finally let's compute a prediction interval for the stopping distance of another car at initial speed 21mph

> Lower.Pred <- Est - sqrt(sigma.tilde.sqd)*sqrt(tau.sqd+1)*qt(0.975,n-p)
> Lower.Pred
> Upper.Pred <- Est + sqrt(sigma.tilde.sqd)*sqrt(tau.sqd+1)*qt(0.975,n-p)
> Upper.Pred

> segments(xstar[1]+0.1,Lower.Pred,xstar[1]+0.1,Upper.Pred)  # x coordinate perturbed by 0.1 for visual clarity
							     # The prediction interval is # wider, because it takes into account the uncertainty in the new observation, as well as 
# our uncertainty in the value of x^{*T} \beta, which is the mean response at x^*.