#
# Duplicates the analysis on the intervals between failures on Page 190 of the book.
#

set.seed(12345)

# Here "midway" is the treatement sample (we expect it to have a larger mean):
#
C = 10000 # at this point we censor
N_U = 2 # there are two censored observations in the midway group
S_U = 3820 + 7170 + 4800

X = S_U + C * N_U # this seems to just be equivalent to the sum of the failure time in the treated group

data = c( 880, 3430, 2860, C, 4750, 3820, C, 7170, C, 4800 )

# Generate the needed combinations of the data:
n_initial = 5 # the sumber of samples in the "initial" class
Perm_Test_Initial_Sample_Index = combn( length(data), n_initial )

# Compute the statistic on each pemutation:
#
n_perms = choose(length(data), n_initial)
perm_stats = rep(NA, n_perms)
for( si in 1:n_perms ){
    indx = Perm_Test_Initial_Sample_Index[, si]
    initial = data[indx]
    treated = data[-indx]

    perm_stats[si] = sum(treated)
}

#postscript("../../WriteUp/Graphics/Chapter11/dup_software_failure_times.eps", onefile=FALSE, horizontal=FALSE)
hist( perm_stats, breaks=20,
      xlab="S_U + C N_C = sum(treated failure times)",
      main="Possible Values for the Test Statistic under the Null Hypothesis" )
abline( v=X, col='red' )
#dev.off()

# What is the probability of getting the value obtained (or a larger one):
#
alpha = sum( perm_stats >= X ) / length( perm_stats )
print( sprintf("alpha (prob. type I error): %10.6f", alpha) )