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

set.seed(1234)

male = c( 107, 110, 116, 114, 113 )
female = c( 110, 111, 111, 108 )

X = mean( male ) - mean( female ) # the statistic observed

# Perform bootstrap replicas:
#
B = 1000
mean_diffs = rep( NA, B )
for( bi in 1:B ){
    b_m = sample( male, size=length(male), replace=TRUE )
    b_f = sample( female, size=length(female), replace=TRUE )
    mean_diffs[bi] = mean(b_m) - mean(b_f)
}

theta_star = X # the best estimate of the difference in the means
se_theta_star = sd( mean_diffs ) # the standard error of the mean difference
alpha = 0.1
z_crit = qnorm( 1 - alpha/2 )
ci = theta_star + c(-1,+1) * z_crit * se_theta_star
print( sprintf("%5.1f%% CI of the difference in mean tooth length is (%4.3f, %4.3f)", (1-alpha)*100, ci[1], ci[2] ) )

# What does the t-test say about this data:
print( t.test( male, female ) )