#
# 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.
#
#-----

save_plots = F

# Ex 9.23 (EPage 341):
#
H = c( 1.2, 0.9, 0.7, 1.0, 1.7, 1.7, 1.1, 0.9, 1.7, 1.9, 1.3, 2.1,
       1.6, 1.8, 1.4, 1.3, 1.9, 1.6, 0.8, 2.0, 1.7, 1.6, 2.3, 2.0 )
P = c( 1.6, 1.5, 1.1, 2.1, 1.5, 1.3, 1.0, 2.6 )

# Make QQ/Normal Plots:
#
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter9/chap_9_ex_23_qqplots.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(1,2))
qqnorm( H ); qqline( H, col=2 ); grid()
qqnorm( P ); qqline( P, col=2 ); grid()
par(mfrow=c(1,1))
if( save_plots ){ dev.off() }

# Make boxplots:
#
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter9/chap_9_ex_23_boxplot.eps", onefile=FALSE, horizontal=FALSE) }
h_df = data.frame( extensibility=H, type=rep( 'high', length(H) ) )
p_df = data.frame( extensibility=P, type=rep( 'poor', length(P) ) )
df = rbind( h_df, p_df )
boxplot( extensibility ~ type, data=df )
if( save_plots ){ dev.off() }

d = mean( H ) - mean( P )
s = sqrt( var(H)/length(H) + var(P)/length(P) )
t = d / s
nu = two_sample_t_test_nu( sd(H), length(H), sd(P), length(P) )
p_value = 2*(1-pt( abs(t), nu )) # for a two-tailed test