## Population Parameters:
mu=4
sigma2=3

reps=1e4
n=30
X = array(rnorm(n*reps, mean=mu, sd=sqrt(sigma2)),c(n,reps))
Xbar = apply(X,2,mean)
EXbar = mean(Xbar)
VXbar = var(Xbar)


MSE.Xbar = mean((Xbar - mu)^2)

## Population Parameters
mu=4
sigma2=3
reps=1e4

MSE=NULL
BIAS2=NULL
VAR=NULL

for ( n in c(5,10,30,100,200,500) )
{
  ## Generate data
  X = array(rnorm(n*reps, mean=mu, sd=sqrt(sigma2)),c(n,reps))
  
  ## Compute the two estimators
  sigmaHat2 = apply(X,2,var)
  sigmaTilde2 = sigmaHat2*(n-1)/n
  
  ## Compute MSE, Bias, and Variance
  MSE.sH2 = mean((sigmaHat2-sigma2)^2)
  BIAS.sH2 = mean(sigmaHat2) - sigma2
  VAR.sH2 = MSE.sH2 - BIAS.sH2^2
  
  MSE.sT2 = mean((sigmaTilde2 - sigma2)^2)
  BIAS.sT2 = mean(sigmaTilde2) - sigma2
  VAR.sT2 = MSE.sT2 - BIAS.sT2^2

  ## store the results:  
  MSE = rbind(MSE , c(n, MSE.sH2, MSE.sT2))
  BIAS2 = rbind(BIAS2, c(n, BIAS.sH2^2,  BIAS.sT2^2))
  VAR = rbind(VAR, c(n, VAR.sH2, VAR.sT2)) 
}