## The number of simulation replications.
nrep = 1e4

## Sample sizes.
SS = seq(5,30,5)

## Storage for the results.
R = array(0, c(length(SS),2))

## Consider sample means for different sample sizes.
for (k in 1:length(SS))
{
  p = SS[k] ## The sample size for this iteration.

  ## Generate p different variances, store as an array.
  V = rexp(p)
  W = matrix(V, nrow=nrep, ncol=p, byrow=TRUE)

  ## Generate a 1000 x p array of normal values, with row i having 
  ## variance V[i].
  X = array(rnorm(nrep*p), c(nrep,p))
  X = sqrt(W) * X

  ## Get the sample mean of each row.
  M = apply(X, 1, mean)

  R[k,1] = var(M)     ## The simulation estimate
  R[k,2] = mean(V)/p  ## The theoretical value
}