#=====================================================
# Figure 13.5.R.
# R script to draw 3d surface (wire mesh plot).
# 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 21X21 grid.
#-----------------------------------------------------
b.int=(0:20)*10/20+30     # Range from 30 to 40.
b.slope=(0:20)*4/20+10    # Range 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 surface.
#-------------------------------------------------------
persp(b.int,b.slope,ss)  # R function for drawing wire mesh surface.
detach(Geyser)