#
# EPage 453
#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
#
#-----

# Ex 12.12 (Part d): 
#
S_x = 517 - 103 - 142
S_y = 346 - 75 - 90
S_x2 = 34095 - 103^2 - 142^2 
S_y2 = 17454 - 75^2 - 90^2 
S_xy = 25825 - 103*75 - 142*90
n = 12

beta_1 = ( S_xy - S_x * S_y / n ) / ( S_x2 - S_x^2 / n ) 
beta_0 = S_y / n - beta_1 * ( S_x / n ) 
SST = S_y2 - S_y^2/n 
SSE = 311.79
print( sprintf("beta_1= %10.6f; beta_0= %10.6f; r2= %10.6f", beta_1, beta_0, 1 - SSE / SST ) )


# Ex 12.13: 
#
xs = c( 20, 40, 60, 80 )
ys = c( .24, 1.2, 1.71, 2.22 )
plot( xs, ys ) 
print( summary( lm( ys ~ xs ) ) ) 


# Ex 12.14: 
#
DF = read.csv( "../../Data/CH12/ex12-01.txt", header=TRUE, quote="'" ) 
plot( DF$Temp, DF$Ratio )

m = lm( Ratio ~ Temp, data=DF )
abline(m) 


# Ex 12.15: 
#
DF = read.csv( "../../Data/CH12/ex12-15.txt", header=TRUE, quote="'" ) 

stem( DF$MOE ) 

m = lm( Strength ~ MOE, data=DF )

plot( DF$MOE, DF$Strength )
abline(m)


# Ex 12.16: 
#
DF = read.csv( "../../Data/CH12/ex12-16.txt", header=TRUE, quote="'" ) 

m = lm( y ~ x, data=DF )
plot( DF$x, DF$y )
abline(m)

print( summary(m) )


# Ex 12.19: 
#
DF = read.csv( "../../Data/CH12/ex12-19.txt", header=TRUE, quote="'" ) 
plot( DF$x, DF$y )

S_x = sum( DF$x ) 
S_x2 = sum( DF$x^2 )
S_y = sum( DF$y ) 
S_y2 = sum( DF$y^2 ) 
S_xy = sum( DF$x * DF$y ) 
n = length(DF$x)

beta_1 = ( S_xy - S_x * S_y / n ) / ( S_x2 - S_x^2 / n ) 
beta_0 = S_y / n - beta_1 * ( S_x / n ) 
print( sprintf("beta_1= %10.6f; beta_0= %10.6f", beta_1, beta_0) )

print( beta_0 + beta_1 * 225 )
print( beta_1 * (-50) ) 


# Ex 12.22
#
n = 15
S_x = 1425
S_x2 = 139037.25
S_y = 10.68
S_y2 = 7.8518
S_xy = 987.645

beta_1 = ( S_xy - S_x * S_y / n ) / ( S_x2 - S_x^2 / n ) 
beta_0 = S_y / n - beta_1 * ( S_x / n ) 

SSE = S_y2 - beta_0 * S_y - beta_1 * S_xy 

eps = sqrt( SSE / (n-2) ) 

SST = S_y2 - S_y^2 / n 

r2 = 1 - SSE / SST 

print( sprintf("beta_1= %10.6f; beta_0= %10.6f; SSE= %10.6f; eps= %10.6f; SST= %10.6f; r2= %10.6f", beta_1, beta_0, SSE, eps, SST, r2) )


# Ex 12.23
#
DF = read.csv( "../../Data/CH12/ex12-19.txt", header=TRUE, quote="'" ) 
m = lm( y ~ x, data=DF ) 

SSE_1 = sum( m$residuals^2 ) # one method for computing the SSE
SSE_2 = sum( DF$y^2 ) - m$coefficients[1] * sum( DF$y ) - m$coefficients[2] * sum( DF$x * DF$y ) # another method for computing the SSE
print( sprintf("SSE_1= %10.6f; SSE_2= %10.6f", SSE_1, SSE_2) )

SST = sum( DF$y^2 ) - sum( DF$y )^2 / n
r2 = 1 - SSE_1 / SST 
print( sprintf("SST= %10.6f r2= %10.6f", SST, r2) )


# Ex 12.24
#
DF = read.csv( "../../Data/CH12/ex12-24.txt", header=TRUE, quote="'" ) 

#postscript("../../WriteUp/Graphics/Chapter12/chap_12_ex_24_scatter_plot.eps", onefile=FALSE, horizontal=FALSE)
plot( DF$x, DF$y )

m_all = lm( y ~ x, data=DF ) 

# Lets drop the sample with the largest x value and refit
m_dropped = print( lm( y ~ x, data=DF[ -dim(DF)[1], ] ) )

abline(m_all, col='red')
abline(m_dropped, col='green')
#dev.off()


# Ex 12.28
#
DF = read.csv( "../../Data/CH12/ex12-20.txt", header=TRUE, quote="'" ) 

rng = range( c( DF$x - mean(DF$x), DF$x ) )

#postscript("../../WriteUp/Graphics/Chapter12/chap_12_ex_28_scatter_plot.eps", onefile=FALSE, horizontal=FALSE)
plot( DF$x, DF$y, col='green', xlim=rng, pch=19 )
lines( DF$x - mean(DF$x), DF$y, col='red', type='p', xlim=rng, pch=19 )
#dev.off()


# Ex 12.29
#
DF = read.csv( "../../Data/CH12/ex12-29.txt", header=TRUE, quote="'" ) 

#postscript("../../WriteUp/Graphics/Chapter12/chap_12_ex_29_scatter_plot.eps", onefile=FALSE, horizontal=FALSE)
par(mfrow=c(1,3))
for( ii in 1:3 ){
  df = DF[ DF$'Data.Set' == ii, ] 
  m = lm( y ~ x, data=df )
  sm = summary(m)
  plot( df$x, df$y, pch=19, main=sprintf("%3d: r2= %3.2f, sigma= %3.2f",ii,sm$r.squared,sm$sigma) )
  abline(m)
}
#dev.off()