Generate normally distributed data with rnorm.

a <- rnorm(10000,mean=65,sd=sqrt(159))
hist(a, col="#0000ff40", border="#0000ff", xlim=c(20,100), prob=T)
lines(density(a), col = "red", lwd = 2)

samplemean<-c()
size<-10
for(i in c(1:10000)){
  samplemean<-c(samplemean,mean(sample(a,size)))
}

hist(samplemean, col="#0000ff40", border="#0000ff", xlim=c(20,100), prob=T)
lines(density(samplemean), col = "red", lwd = 2)

samplemean<-c()
size<-10
for(i in c(1:10000)){
  samplemean<-c(samplemean,mean(sample(a,size)))
}

# Mean of the sample means
mean(samplemean)
## [1] 64.82349
# Variance of the sample means
sum((samplemean-mean(samplemean))^2)/length(samplemean)
## [1] 15.36882
# Population variance divided by the size of a sample
159/size
## [1] 15.9