set.seed(12345)

C = 24

drug=c( 23, C, C, C )
placebo = c( 15, 18, 19, 19, 20 )

X = mean(drug) - mean(placebo)

data = c( drug, placebo )

# Generate the needed combinations of the data:
n_drug = 5 # the sumber of samples in the "drug" class
Perm_Test_Drug_Sample_Index = combn( length(data), n_drug )

# Compute the statistic on each pemutation:
#
n_perms = choose(length(data), n_drug)
perm_stats = rep(NA, n_perms)
for( si in 1:n_perms ){
    indx = Perm_Test_Drug_Sample_Index[, si]
    boot_drug = data[indx]
    boot_placebo = data[-indx]

    perm_stats[si] = mean(boot_drug) - mean(boot_placebo)
}

hist( perm_stats, breaks=20,
     xlab="mean(drug) - mean(placebo)",
     main="Possible Values for the Test Statistic under the Null Hypothesis" )
abline( v=X, col='red' )

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