## Column 1 is the plug-in Z-statistic p-value, column 2 is the
## t-statistic p-value, column 3 is the proportion of 1000 null
## plug-in Z-statistics exceeding the actual plug-in Z-statistic, and
## column 4 is the proportion of simulated null t-statistics exceeding
## the actual t-statistic.
nrep = 1e3
Q = array(0, c(nrep,4))

## Sample sizes for the two groups.
m = 10
n = 5

## Simulate nrep null data sets.
X = array(rnorm(n*nrep), c(nrep,n))
Y = array(rnorm(m*nrep), c(nrep,m))

## Calculate plug-in Z-statistics for the simulated data sets.
MD = apply(X, 1, mean) - apply(Y, 1, mean)
VX = apply(X, 1, var)
VY = apply(Y, 1, var)
TZS = MD / sqrt(VX/n + VY/m)

## Calculate t-statistics for the simulated data sets.
Sp2 = ((n-1)*VX + (m-1)*VY) / (n+m-2)
TTS = MD / sqrt((m+n)*Sp2/(m*n))

## Do nrep replications.
for (r in 1:nrep)
{
  ## Population means and variances for populations A and B.
  ## The minimum separation between the means is 1.
  M = c(2, runif(1))
  V = 2*runif(2)

  ## Generate the sample data.
  X = M[1] + sqrt(V[1])*rnorm(n)
  Y = M[2] + sqrt(V[2])*rnorm(m)

  ## Get the plug-in Z-statistic p-value.
  MD = mean(X) - mean(Y)
  VX = var(X)
  VY = var(Y)
  TZ = MD / sqrt(VX/n + VY/m)
  Q[r,1] = 1-pnorm(TZ)

  ## Get the T-statistic p-value.
  Sp2 = ((n-1)*VX + (m-1)*VY) / (n+m-2)
  TT = MD / sqrt(Sp2*(m+n)/(m*n))
  Q[r,2] = 1-pt(TT, n+m-2)

  ## Get the proportion of null test statistics exceeding the actual
  ## test statistic, for both the plug-in Z statistic (column 3) and
  ## the t-statistic (column 4).
  Q[r,3] = mean(TZS > TZ)
  Q[r,4] = mean(TTS > TT)
}