#
# Duplicates the analysis on the one-way table 8.1:
# The Effects of Fertilizer on Crop Yield on Page 137 of the book.
#
DF = data.frame(
    fertilizer=as.factor( c( rep( 'FastGro', 5 ), rep( 'NewGro', 3 ), rep( 'WunderGro', 3 ) ) ),
    yield=c( 27, 30, 55, 72, 18, 38, 12, 72, 75, 76, 54 )
)

# To start lets apply the standard ANOVA tests (from R) to see if there is a difference in group means:
#
res = anova( lm( DF$yield ~ DF$fertilizer ) )
print( sprintf("Classical ANOVA: prob. type I error: %10.6f", res$"Pr(>F)"[1]) )

# Apply a permutation test:
#
source('utils.R') # to get the function compute_F

F = compute_F( DF$fertilizer, DF$yield ) # the value of F under the original labeling

set.seed(1234) 
B = 500 # the sumber of bootstraps to run
all_Fs = rep( NA, B ) # the bootstrap values of the F_1
for( bi in 1:B ){
    all_Fs[bi] = compute_F( sample(DF$fertilizer), DF$yield ) # reorder the labels
}

plot( density( all_Fs ), xlab="statistic", main="" )
grid()
abline(v=F, col='red')

alpha = sum( all_Fs >= F ) / B
print( sprintf("Permutation Bootstrap (B: %d): prob. type I error: %10.6f", B, alpha) )