source('../../Data/data_loaders.R')

DF = load_exercise_3_8_data()

# Part (a):
#
m = lm( Y ~ X, data=DF )

# Lets plot the data and models:
#
#postscript("../../WriteUp/Graphics/Chapter3/ex_3_8_part_a.eps", onefile=FALSE, horizontal=FALSE)
par(mfrow=c(1,4))
plot( DF$X, DF$Y, type='p', pch=19, cex=1.5, xlab='X', ylab='Y', main='data and linear model' )
abline(m, col='red')
grid()
plot( m$fitted.values, m$residuals, type='p', pch=19, cex=1.5, xlab='y_hat', ylab='residual', main='residual vs. fitted' )
grid()
plot( m, which=2 )
grid()
plot( m$fitted.values, DF$Y, xlab='y_hat', ylab='y', main='y vs. y_hat' )
grid()
par(mfrow=c(1,1))
#dev.off()


# Part (b):
#
library(MASS)
bc = boxcox(m)

# What lambda value gives the largest log-likelihood:
#
print(sprintf("Boxcox lambda estimate: %10.6f", bc$x[ which.max(bc$y) ]))

# Lets try Atkinsons method:
#
source('utils.R')

atk_res = Atkinsons( DF, 'Y', c( 'X' ) )

print(sprintf("Atkinsons lambda estimate: %10.6f (t-value %10.6f)", atk_res$lambda_hat, atk_res$gamma_t))


# Part (c):
#
DF$one_over_Y = 1/DF$Y
m_one_over = lm( one_over_Y ~ X, data=DF )

# Make some predictions with different models:
p1 = predict( m, newdata=data.frame( X=3 ), level=0.95, interval="prediction" )
p = predict( m_one_over, newdata=data.frame( X=3 ), level=0.95, interval="prediction" )
p2 = 1/p
t = p2[2]; p2[2] = p2[3]; p2[3] = t; # swap upper and lower bounds


# Put everything together:
#
R = rbind( p1, p2 )
rownames(R) = 1:dim(R)[1]
R = data.frame( R )
R$pi_range = R$upr - R$lwr
print( R )

# Plot the data and the three models:
#
#postscript("../../WriteUp/Graphics/Chapter3/ex_3_8_part_c.eps", onefile=FALSE, horizontal=FALSE)
plot( DF$X, DF$Y, type='p', pch=19, cex=1.5, xlab='X', ylab='Y' )
all_Xs = unique( DF$X )
points( all_Xs, predict( m, newdata=data.frame( X=all_Xs ) ), type='l', col='blue' )
p = predict( m_one_over, newdata=data.frame( X=all_Xs ) )
points( all_Xs, 1/p, type='l', col='green' )
legend( 'topleft',
       c('Y = beta_0 + beta_1 X', 'Y=1/(beta_0 + beta_1 X)'),
       col=c('blue', 'green'),
       lty=c(1, 1) )
grid()
#dev.off()