nsamp = 20     ## The number of data values in each data set
nrep = 1000    ## The number of data sets to generate
pop_mean = 0   ## The population mean of each data value
pop_var = 1    ## The population variance of each data value

## Generate a nrep x nsamp array of standard normal draws.
D = rnorm(nrep*nsamp, mean=pop_mean, sd=sqrt(pop_var)) ## **
X = array(D, c(nrep,nsamp))

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

## Calculate the sample variance of Y.
V = var(Y)

## Generate a histogram of the raw data.
h1 = hist(D[1:nrep])

## Generate a histogram of the sample means.
h2 = hist(Y)

## Plot the raw data histogram in blue.
plot(h1, col='blue', ylim=c(0,max(h1$counts, h2$counts)), 
     xlab='', main='')

## Overplot the sample means histogram in red.
lines(h2, col='red')

## Add a legend to the plot.
legend(x='topright', legend=c('Raw data', 'Averages'), 
       col=c('blue', 'red'), lty=c(1,1))