# Doing it the hard way  -- only because I can't make the script using apply
#   come out with believable results.

# The data are set up with the ID is col 1, the observations in col 2,
# and the group membership in col 3.
# For Example:
#ID Arousal Group
#01 39.1 1
#02 38.0 1
#03 14.9 1

# There are column headings, but these are then ignored after they are read.
# The p value that is given is the two-tailed p.
# For an example with 35 cases in one group and 29 in the other, with
# nreps = 10,000, this took 21 seconds.

data <- read.table(file.choose(), header = TRUE)
#Sys.time()  Uncomment this line (and the last line) if you want to time the run
x <- data[,2]
group <- data[,3]
nrep <- 100000
t <- numeric(nrep)
N <- length(x)
n1 <- length(x[group == 1])
n2 <- length(x[group == 2])

 # Now set up function to compute t values
 ttest <- function(x,group, n1, n2)  {

   mean1 <- mean(x[group == 1])
   mean2 <- mean(x[group == 2])
   var1 <- var(x[group == 1])
   var2 <- var(x[group == 2])
   s2pooled <- ((n1-1)*var1 + (n2-1)*var2)/(n1 + n2 - 2)
   t <- (mean1 - mean2)/sqrt(s2pooled/n1 + s2pooled/n2)
   return(t)
   }
 tobt <- ttest(x,group,n1,n2)


start <- proc.time()
  xx <- numeric(N)
  t <- numeric(nrep)
 for (i in 1:nrep) {
    xx <- sample(x)
    t[i] <- ttest(xx,group, n1, n2)
    }
finish <- proc.time()
time <- finish-start
cat("Timing:  ",time, "\n")
    

tneg <- -1*abs(tobt)
tpos <- abs(tobt)
lower <- length(t[t <= tneg])
upper <- length(t[t >= tpos])
lowerp <- lower/nrep
upperp <- upper/nrep
twotailed <- lowerp + upperp
cat("Lower one-tailed probability is ", lowerp ,"\n")
cat("Upper one-tailed probability is ", upperp, "\n" )
cat("The two-tailed probability is ", twotailed," \n")

# Set up histogram
dist <- seq(-6, 6,.1)
hist(t, breaks = 50, freq = F, xlim = c(-5,5),main = "Distribution of t for randomization",
  col = "gray")
  # Superimpose the t distribution just for fun
quant <- seq(min(t), max(t), length = 100)
y <- dt(quant,  df = 41, ncp = 0)
lines(quant,y)
  
legend(2, .3, paste("t = ", round(tobt, digits = 4)), col=1,  cex = 0.8)
legend(2, .25, paste("p = ", round(twotailed, digits = 4)), col=1, cex = 0.8)

# Sys.time()