X <- read.csv("U:/public_html/588/Class1/CASE0101.csv", header = TRUE)     # replace by your own path here

V <- X[,1]    # all the scores in one "pile"
n1 <- 23
n2 <- 24

choose(47,23)      # how many combinations are there

M <- 10000            # size of the Monte Carlo study
Dstat <- rep(0, M)    # vector of simulation results
Tstat <- rep(0, M)    

for (mc in 1:M){
  Y <- sample(V, length(V), replace = FALSE)     # this is the permutation of the original data
  Y1 <- Y[1:n1]
  Y2 <-  Y[(n1+1):(n1+n2)]            # two parts of the sample
  Dstat[mc] <- mean(Y2) - mean(Y1)    # calculating the test statistic         
  # Tstat[mc] <- Dstat[mc]/sqrt( (var(Y1)*(n1-1) + var(Y2)*(n2-1))/(n1+n2-2))/sqrt(1/n1 + 1/n2)                            
}

hist(Dstat)
count <- sum( abs(Dstat) > 4.14)    # this provides the count of "exceptional" cases  
pvalue <- count/M
print( paste( "P-value for the difference statistic", pvalue))


hist(Tstat)
count2 <- sum( abs(Tstat) > 2.92)    
pvalue2 <- count2/M
print( paste( "P-value for the T statistic", pvalue2))