x <- c(22.0,23.9,20.9,23.8,25.0,24.0,21.7,23.8,22.8,23.1) y <- c(23.2,22.0,22.2,21.2,21.6,21.9,22.0,22.9,22.8) round(c(x.mean=mean(x),x.var=var(x),y.mean=mean(y),y.var=var(y)),1) n.sim <- 10000; m <- length(x); n <- length(y) phi <- (m-1)*var(x)/rchisq(n.sim,df=m-1) psi <- (n-1)*var(y)/rchisq(n.sim,df=n-1) delta <- rnorm(n.sim,(mean(x)-mean(y)),phi/m+psi/n) mean(delta>0)