# template.txt    Regression Artifact example
# Testing improvement for top 50 and bottom 50
# So df=49, and critical value is qt(0.95,49) = 1.676551
rm(list=ls())
source("https://www.utstat.toronto.edu/brunner/Rfunctions/rmvn.txt")
outfile = "simoutX.txt" # Upper case x will be replaced by 1, 2, 3, 4
nsim=5000 # Monte Carlo sample size
N = c(100, 200, 500, 1000); nn = length(N)
RHO = seq(from=0.0, to=0.9, by = 0.1); nrho = length(RHO)
output = matrix(NA,nsim*nn*nrho,4) 
colnames(output) = c("n","Correlation","High_t","Low_t")

rownumber = 0
for(i in 1:nn)
    {
    for(j in 1:nrho)
        {
        for(sim in 1:nsim)
            {
            rownumber = rownumber+1
            n = N[i]; rho = RHO[j]
            kor = rbind(c(1,rho),  # Correlation matrix of the 
                        c(rho,1))  # "Before" and "After" measurements.
            x = rmvn(n,c(0,0),kor) # Col 1 is "Before," col 2 is "After."
            # Sort on before measurements, smallest to largest
            orderbefore = order(x[,1])
            x = x[orderbefore,] 
            x = round(15*x + 100) # mu=100, sigma=15, like IQ scores
            improvement = x[,2]-x[,1] # After minus before
            elite = improvement[(n-49):n]   # Improvement of the top 50
            disad = improvement[1:50]       # Improvement of the bottom 50
            topt = t.test(elite); 
            bottomt = t.test(disad)
            t1 = topt$statistic; t2 = bottomt$statistic
            thisrow = c(n,rho,t1,t2)
            output[rownumber,] = thisrow
            } # Next simulation 
        cat("n = " ,n, "  Correlation = ", rho, "\n")
        } # Next rho
    } # Next n
write.table(output,outfile,quote=F,col.names=F)