#
# 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.30: 
#
xs = c( 1000, 1500, 2000, 2500, 3000, 3500, 4000 )
print( mean(xs) )
S_xx = sum( ( xs - mean(xs) )^2 ) 
print( S_xx )
sigma = 350
sigma_beta_1 = sigma/sqrt(S_xx) 
print( sigma_beta_1 )

n = 7
print( pt( ( 1.5 - 1.25 )/sigma_beta_1, n-2 ) - pt( ( 1.0 - 1.25 )/sigma_beta_1, n-2 ) )

xs = seq( 2000, 3000, by=100 ) 
S_xx = sum( ( xs - mean(xs) )^2 ) 


# Ex 12.31: 
#
n = 15 # numbers taken from 12.18
S_x = 1425
S_y = 10.68
S_x2 = 139037.25
S_xy = 987.645
S_y2 = 7.8518

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 

s = sqrt( SSE / (n-2) ) # the estimate of sigma
 
s_beta = s / sqrt( S_x2 - S_x^2/n ) # Part (a): the estimate of sigma_{\hat{\beta}}

t_value = qt( 1-0.025, 13 )
beta_1 + c(-1,+1) * t_value * s_beta # Part (b): CI on beta_1 


# Ex 12.32: 
#
t_value = qt( 1-0.025, 16-2 )

ci = 0.82697 + c(-1,+1) * t_value * (0.03652)


# Ex 12.33: 
#
t_value = qt( 1-0.025, 27-2 )

ci = 0.10748 + c(-1,+1) * t_value * 0.0128


# Ex 12.34: 
#
t_value = qt( 1-0.005, 13-2 )

ci = 0.9668 + c(-1,+1) * t_value * 0.1829


# Ex 12.35: 
#
n = 17
S_x = 222.1
S_y = 193
S_x2 = 3056.69
S_xy = 2759.6
S_y2 = 2975

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 

s = sqrt( SSE / (n-2) ) # the estimate of sigma
 
s_beta = s / sqrt( S_x2 - S_x^2/n ) # the estimate of sigma_{\hat{\beta}}

t_value = qt( 1-0.025, n-2 )
CI = beta_1 + c(-1,+1) * t_value * s_beta # CI on beta_1 

t_value = beta_1 / s_beta

1 - pt( t_value, n-2 ) 

# maybe the definition includes an absolute value?
# print( 1 - ( pt( t_value, n-2 ) - pt( -t_value, n-2 ) ) )
#

# Drop the sample (6.0,2.5)
n = n-1
S_x = S_x - 6.0
S_y = S_y - 2.5
S_x2 = S_x2 - 6.0^2 
S_y2 = S_y2 - 2.5^2
S_xy = S_xy - 6.0*2.5

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 

s = sqrt( SSE / (n-2) ) # the estimate of sigma
 
s_beta = s / sqrt( S_x2 - S_x^2/n ) # the estimate of sigma_{\hat{\beta}}

t_value = qt( 1-0.025, n-2 )
CI = beta_1 + c(-1,+1) * t_value * s_beta # CI on beta_1 


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

m = lm( y ~ x, data=DF ) 
sm = summary(m)
print(sm)

beta_1_times_900 = m$coefficients[2] * 900 
se_beta_1_times_900 = sm$coefficients[2,2] * 900 
t_value = qt( 1-0.0025, length(DF$x)-2 ) 
CI = beta_1_times_900 + c(-1,+1) * t_value * se_beta_1_times_900 



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

# Part (a): Lets look if the difference between the two levels is significant: 
t_stat = mean( DF$Level.2 - DF$Level.1 ) / ( sd( DF$Level.2 - DF$Level.1 ) / sqrt( length(DF$Level.2) ) )
print( t_stat ) 

# Part (b): Regresss one level on the other level: 
plot( DF$Level.1, DF$Level.2 )

m = lm( Level.2 ~ Level.1, data=DF ) 
sm = summary(m)
print(sm)

# Lets compute the 95% confidence interval between the two Levels: 
#
t_value = qt( 1-0.0025, length(DF$Level.1)-2 ) 
CI = m$coefficients[2] + c(-1,+1) * t_value * sm$coefficients[2,2]
print( CI ) 



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

m = lm( DF$y ~ DF$x, data = DF ) 
sm = summary(m)

t_value = qt( 1-0.025, length(DF$x)-2 ) 
10 * sm$coefficients[2,1] + c(-1,+1) * t_value * 10 * sm$coefficients[2,2] 


# Ex 12.39: 
#
SSE = 7.968
SST = 124039.58 - 1574.8^2 / 20 
SSR = SST - SSE 
n = 20

f_value = SSR / ( SSE/(n-2) ) 

F_alpha_1_n_minus_2 = qf( 1-0.05, 1, n-2 )

qf( 1-0.05, 1, 20-2 ) 
1 - pf( 71.97, 1, 20-2 ) 

beta_1 = 0.04103377
S_x = 2817.9
S_x2 = 415949.85
n = 20 
S_xx = S_x2 - S_x^2 / n 
SSE = 7.968 
s = sqrt( SSE / (n-2) )
s_beta = s / sqrt( S_xx ) 
t_stat = beta_1 / s_beta 
t_stat^2 

# Lets compute the significance probablity using the t-distribution: 
1 - ( pt( t_stat, n-2 ) - pt( -t_stat, n-2 ) ) 


# Ex 12.43: 
#
numer = abs( 1 - 2 ) 
denom = 4 * sqrt( (15 - 1)/( 11098 - 402^2/15 ) )
d = numer / denom
print(d)