set.seed(123456)
par(mfrow=c(2,2))

n1 <- 10
n2 <- 20
n3 <- 30
n4 <- 100

reps <- 10000

vbar10<-numeric(10000)
samples_1 <- replicate(reps, rnorm(n1,10,2)) # 10 x 10000 sample matrix
vbar10 <- rep(NA, reps)
for (i in 1:reps) {
  vbar10[i] <- var(samples_1[,i])
}
mean(vbar10)
var(vbar10)
hist(vbar10, freq=F, xlab="", breaks=100,xlim=c(0,10), ylim=c(0,0.8))
par(new=T)
plot(density(vbar10), axes=F, xlim=c(0,10), ylim=c(0,0.8), main="", col="blue")

vbar20<-numeric(10000)
samples_2 <- replicate(reps, rnorm(n2,10,2)) # 20 x 10000 sample matrix
vbar20 <- rep(NA, reps)
for (i in 1:reps) {
  vbar20[i] <- var(samples_2[,i])
}
mean(vbar20)
var(vbar20)
hist(vbar20, freq=F, xlab="", breaks=100,xlim=c(0,10),ylim=c(0,0.8))
par(new=T)
plot(density(vbar20), axes=F, xlim=c(0,10), ylim=c(0,0.8),main="", col="blue")

vbar30<-numeric(10000)
samples_3 <- replicate(reps, rnorm(n3,10,2)) # 30 x 10000 sample matrix
vbar30 <- rep(NA, reps)
for (i in 1:reps) {
  vbar30[i] <- var(samples_3[,i])
}
mean(vbar30)
var(vbar30)
hist(vbar30, freq=F, xlab="", breaks=100,xlim=c(0,10), ylim=c(0,0.8))
par(new=T)
plot(density(vbar30), axes=F, xlim=c(0,10), ylim=c(0,0.8),main="", col="blue")

vbar100<-numeric(10000)
samples_4 <- replicate(reps, rnorm(n4,10,2)) # 100 x 10000 sample matrix
vbar100 <- rep(NA, reps)
for (i in 1:reps) {
  vbar100[i] <- var(samples_4[,i])
}
mean(vbar100)
var(vbar100)
hist(vbar100, freq=F, xlab="", breaks=100,xlim=c(0,10), ylim=c(0,0.8))
par(new=T)
plot(density(vbar100), axes=F,xlim=c(0,10), ylim=c(0,0.8), main="", col="blue")