#
# 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.
#
# Computations used for the problems in this section of the text.
#
#-----

source('utils.R')

print( 'Ex 10.22:' )
DF = read.csv( "../../Data/CH10/ex10-22.txt", header=TRUE, quote="'" )
res = ANOVA_USS( DF$yield, DF$EC )


print( 'Ex 10.23:' )
xbar_i_dot = res$x_i_dot / res$Js # the group means
mean_differences = outer( xbar_i_dot, xbar_i_dot, "-" )
print( mean_differences )

q_crit = qtukey( 1-0.05, res$I, res$n-res$I )
vs = ( res$MSE / 2 ) * outer( 1/res$Js, 1/res$Js, "+" )

ws = q_crit * sqrt( vs ) # take only 1/2 of this matrix (the top or the bottom)
print( ws )


print( 'Ex 10.24:' )
I = 3
n = 74
MSTr = 76.09
SST = 1123.14

SSTr = (I-1)*MSTr
SSE = SST - SSTr
MSE = SSE / ( n-I )

f = MSTr / MSE

p_value = 1 - pf( f, I-1, n-I )



print( 'Ex 10.25:' )
Js = c( 8, 13, 17, 14 )
I = length(Js)
n = sum(Js)
xbars = c( 43.0, 42.4, 43.1, 43.5 )
ss = c( 1.5, 1.3, 1.2, 1.2 )

x_i_dot = xbars * Js
x_dot_dot = sum(x_i_dot)

SSTr = sum( x_i_dot^2 / Js ) - x_dot_dot^2 / n
MSTr = SSTr / (I-1)
SSE = sum( (Js-1) * ss^2 )
MSE = SSE / (n-I)

f = MSTr / MSE

p_value = 1 - pf( f, I-1, n-I )


print( 'Ex 10.26:' )
DF = read.csv( "../../Data/CH10/ex10-26.txt", header=TRUE, quote="'" )
res = ANOVA_USS( DF$PAPFUA., DF$Type )

# Get the confidence intervals for the differences in means (following the code for Ex. 10.23)
#
u_ci = mean_diff_CI_utils( res$x_i_dot, res$Js, res$MSE )
print( u_ci )

# Get the CI for the linear combination: 
cs = c( rep( 0.25, 4 ), rep( -0.5, 2 ) )
xbars = res$x_i_dot / res$Js # the means in the original sorted order

m = sum( cs * xbars )
v = res$MSE * sum( cs^2 / res$J )
t_crit = qt( 1 - 0.01/2, res$n - res$I )
print( "CI" )
print( m + sqrt( v ) * t_crit * c(-1,+1) )


print( 'Ex 10.27:' )

DF = read.csv( "../../Data/CH10/ex10-27.txt", header=TRUE, quote="'" ) 
res = ANOVA_USS( DF$folacin, DF$Brand.Tea )

# Get CI on the differences in the means
#
u_ci = mean_diff_CI_utils( res$x_i_dot, res$Js, res$MSE )
print( u_ci )


print( 'Ex 10.30:' )
phi2 = 2 * ( 4*(1/5)^2 + (4/5)^2 )
phi = sqrt(phi2)


print( 'Ex 10.32:' )
DF = read.csv( "../../Data/CH10/ex10-32.txt", header=TRUE, quote="'" ) 
DF$Flaws = sqrt( DF$Flaws ) # perform the needed transformation 
res = ANOVA( DF$Flaws, DF$Brand.Tape )