#=====================================================
# Figure 13.6.R.
# R script to draw contour plot of 3d surface.
# Example is sum of squares surface for slope and intercept
# of prediction line for Old Faithful data.
#=====================================================

#-----------------------------------------------------
# Input the data.
#-----------------------------------------------------
Geyser=read.table("old_faithful.txt",header=TRUE) # Change file name
                                                  #  if necessary.
attach(Geyser)

#-----------------------------------------------------
# Set up X matrix.
#-----------------------------------------------------
n=length(y)               # Number of observations is n.
X=matrix(1,n,2)           # Form the X matrix:  col 1 has 1’s,
X[,2]=x                   #  col 2 has predictor variable.

#-----------------------------------------------------
# Calculate range of intercept (b.int) and slope (b.slope)
# values.  Use 201X201 grid.
#-----------------------------------------------------
b.int=(0:200)*10/200+30   # Range from 30 to 40.
b.slope=(0:200)*4/200+10  # Range is from 10 to 14.

#-----------------------------------------------------
# Calculate matrix (ss) of sum of squares values.  Rows
# correspond to intercept values, columns correspond to
# slope values.
#-----------------------------------------------------
ss=matrix(0,length(b.int),length(b.slope))
for (i in 1:length(b.int)) {
  for (j in 1:length(b.slope)) {
    b=rbind(b.int[i],b.slope[j]) # Col vector of intercept, slope.
    ss[i,j]=sum((y-X%*%b)^2)     # Sum of squares;  X%*%b is col vec
                                 #  of predicted values.
  }
}

#-------------------------------------------------------
# Draw the sum of squares contour plot.  Calculate least
# squares solution for intercept and slope and add solution
# point to plot. 
#-------------------------------------------------------
contour(b.int,b.slope,ss,nlevels=30,xlab="intercept",ylab="slope")
bhat=solve(t(X)%*%X,t(X)%*%y)
points(bhat[1],bhat[2])
detach(Geyser)