m <- 50
n <- 50
mu <- 0
sigma1 <- 1
sigma2 <- 1

vals <- seq(-3,3,0.01)

numIter <- 1000
testStat <- rep(0,numIter);


for (i in 1:numIter) {

	x <- rnorm(m,mu,sigma1);
	y <- rnorm(m,mu,sigma2);

	testStat[i] <- (mean(x)-mean(y))/((var(x)/m) + (var(y)/n))^0.5


	df <- ((var(x)/m) + (var(y)/n))^2 / 
	 	(   ((m-1)^(-3) * (var(x)*(m-1)/m)^2) +
		    ((n-1)^(-3) * (var(y)*(n-1)/n)^2)   )
	
	cat(df,"\n");	
	plot(vals,dt(vals,df),type="l",ylim=c(0,0.5),col=1);
	par(new=TRUE)
}


plot(density(testStat,bw=0.2,from=-3,to=3),ylim=c(0,0.5),col=2)