x<-RocketProp$Age
y<-RocketProp$Strength
n<-length(x)

plot(x,y,pch=19,xlab='Age',ylab='Shear Strength')
fit.RP<-lm(y~x)
summary(fit.RP)

coef(fit.RP)
abline(coef(fit.RP),col='red')
title('Line of best fit for Rocket Propulsion Data')

####################################################
#How does the sum-of-squares vary for different beta values ?

be0<-seq(2600,2650,by=2)
be1<-seq(-40,-35,by=0.05)

Resid.func<-function(b0,b1,xv,yv){

	return((yv-(b0+b1*xv))^2)


}

Smat<-outer(be0,be1,'*')*0
for(i in 1:n){

	Smat<-Smat+outer(be0,be1,Resid.func,xv=x[i],yv=y[i])

}

zval<-quantile(Smat,prob=c(1:19)/20)
image(be0,be1,Smat,col=rev(heat.colors(200)))
contour(be0,be1,Smat,add=T,levels=zval)
abline(v=coef(fit.RP)[1],h=coef(fit.RP)[2])

be.hat<-coef(fit.RP)
be0.hat<-be.hat[1]
be1.hat<-be.hat[2]
persp(be0,be1,Smat,theta = 60, phi = 20,zlab='S',ticktype='detailed')->res
points(trans3d(be0.hat, be1.hat, min(Smat), pmat = res), col = 2, pch = 16,cex=2)