#
# 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 13.15:
#
DF = read.csv( "../../Data/CH13/ex13-15.txt", header=TRUE, quote="'" ) 
plot( DF$x, DF$y )

plot( log(DF$x), log(DF$y) )

DF$lx = log(DF$x)
DF$ly = log(DF$y)
m = lm( ly ~ lx, data=DF )

preds = predict( m, newdata=data.frame(lx=log(20)), interval='prediction', level=0.95 )
exp(preds)

# Look at the residual plots in the transformed space:
#
plot( m$fitted.values, rstandard(m) )



# Ex 13.16:
#
DF = read.csv( "../../Data/CH13/ex13-16.txt", header=TRUE, quote="'" ) 
plot( DF$Load, DF$Time ) # this has y values that are spread all over (one large many small)
plot( DF$Load, log(DF$Time) ) # taking the log compresses the scale of the y axis

DF$lt = log(DF$Time)

m = lm( lt ~ Load, data = DF )
abline(m)

sm = summary(m)

# Make a prediction of y' when load percentage is given
pred = predict( m, newdata=data.frame(Load=80), interval='prediction', level=0.95 )

# Undo the nonlinear transformation to get y:
exp(pred)


# Ex 13.17:
#
# EPage 517
# 
DF = read.csv( "../../Data/CH13/ex13-17.txt", header=TRUE, quote="'" ) 
plot( DF$x, DF$y )


plot( log(DF$x), log(DF$y) )

DF$lx = log(DF$x)
DF$ly = log(DF$y)

m = lm( ly ~ lx, data=DF )
sm = summary(m) 

plot( m$fitted.values, rstandard(m) ) # look at a residual plot 

beta_1_0 = 4/3 # the hypothetical value for beta_1 i.e. H_0: \beta_1 = \beta_1_0

t = ( sm$coefficients[2,1] - beta_1_0 ) / sm$coefficients[2,2]

n = length(DF$x)
alpha = 0.05
t_crit = qt( 1 - alpha, n-2 )


beta_1_0 = 1 # the hypothetical value for beta_1 i.e. H_0: \beta_1 = \beta_1_0

t = ( sm$coefficients[2,1] - beta_1_0 ) / sm$coefficients[2,2]

n = length(DF$x)
alpha = 0.05
t_crit = qt( 1 - alpha/2, n-2 )



# Ex 13.18:
#
# EPage 518
# 
DF = read.csv( "../../Data/CH13/ex13-18.txt", header=TRUE, quote="'" ) 
plot( DF$Cycfail, DF$Strampl ) # note the large dynamic range of Cycfail 

# these two models look very similar:
# 
plot( log(DF$Cycfail), DF$Strampl )
plot( log(DF$Cycfail), log(DF$Strampl) )

# compute the R^2 of each model and determine which one is larger:
#
DF$lCycfail = log(DF$Cycfail)
m = lm( Strampl ~ lCycfail, data=DF )
summary( m )$r.squared 
summary( lm( log(Strampl) ~ log(Cycfail), data=DF ) )$r.squared

# plot a residual plot:
#
plot( m$fitted.values, rstandard(m) ) 

# make predictions:
#
pred = predict( m, newdata=data.frame(lCycfail = log(5000.)), interval='prediction', level=0.95 )
print( pred )


# Ex 13.19:
#
# EPage 518
# 
DF = read.csv( "../../Data/CH13/ex13-19.txt", header=TRUE, quote="'" ) 
plot( DF$Temp, DF$Lifetime )

plot( 1/DF$Temp, log( DF$Lifetime ) )

DF$iTemp = 1/DF$Temp
DF$lLifetime = log( DF$Lifetime )

m = lm( lLifetime ~ iTemp, data=DF )

# Make the required prediction:
# 
exp( predict( m, newdata=data.frame(iTemp = 1/220) ) )

# Look at whether we can reject H_0 (a linear model is adequate):
#
source('H0_linear_model.R')
H0_linear_model( DF$iTemp, DF$lLifetime )



# Ex 13.20:
#
# EPage 518
# 
DF = read.csv( "../../Data/CH13/ex13-14.txt", header=TRUE, quote="'" ) 

# Plot each of the suggested transformations
x = DF$x
y = DF$y

plot( x, y ) # seems to have power low or exponential decay 
plot( log(x), y ) # just changes x scale
plot( x, log(y) ) # better at linearizing
plot( log(x), log(y) )
plot( 1/x, y )
plot( 1/x, log(y) )



# Ex 13.21:
#
# EPage 518
# 
DF = read.csv( "../../Data/CH13/ex13-21.txt", header=TRUE, quote="'" ) 

plot( DF$x, DF$y )
plot( 10^4 * (1/DF$x), DF$y )

DF$tx = 10^4 * (1/DF$x) # our transformed x value

m = lm( y ~ tx, data=DF )
summary(m)

predict( m, newdata=data.frame(tx=(10^4)/500) )