#======================================================================
# Figure 12.6.R.
# R script to calculate intercept b[1] and slope b[2] for the
# least squares prediction line.  Script produces a scatterplot
# of the data along with overlayed prediction line.
#
# Data is apparent magnitudes and periods of nine Cepheid variable
# stars in the Small Magellanic Cloud.
#
# Data file should have two columns:  response variable to be
# predicted labeled "y" and predictor variable labeled "x".
# Change the R working directory to the location of the data file. 
# Re-label axes in the plot( ) function.
#====================================================================== 

smc=read.table("cepheids_smc.txt",header=TRUE)  # Change file name if
                                                   #  necessary.
attach(smc)
y=apmag
x=log(period)
n=length(y)               # Number of observations is n.
X=matrix(1,n,2)           # Form the
X[,2]=x;                  #  X matrix.
b=solve((t(X)%*%X))%*%t(X)%*%y   # Least squares estimates in b;
                                 #  t() is transpose.
plot(x,y,type="p",xlab="log(period)",
  ylab="apparent magnitude")  #Scatterplot of data.
yhat1=b[1]+b[2]*min(x)           # Calculate predicted y values at 
yhat2=b[1]+b[2]*max(x)           #  smallest and largest values of x.
yhat=rbind(yhat1,yhat2)
xvals=rbind(min(x),max(x))
points(xvals,yhat,type="l")    # Connect two predicted values with line.
"least squares intercept and slope:"  # Print the intercept and slope
b                                     #  to the console.
detach(smc)