set.seed(1)
n1 <- 2
n2 <- 5
n3 <- 10
n4 <- 30 
reps <- 10000

# perform random sampling
samples_1 <- replicate(reps, rnorm(n1,2.5,1.118)) # 2 x 10000 sample matrix
samples_1_mean <- rep(NA, reps)
for (i in 1:reps) {
  samples_1_mean[i] <- mean(samples_1[,i])
}
mean(samples_1_mean)
var(samples_1_mean)

samples_2 <- replicate(reps, rnorm(n2,2.5,1.118)) # 5 x 10000 sample matrix
samples_2_mean <- rep(NA, reps)
for (i in 1:reps) {
  samples_2_mean[i] <- mean(samples_2[,i])
}
mean(samples_2_mean)
var(samples_2_mean)

samples_3 <- replicate(reps, rnorm(n3,2.5,1.118)) # 10 x 10000 sample matrix
samples_3_mean <- rep(NA, reps)
for (i in 1:reps) {
  samples_3_mean[i] <- mean(samples_3[,i])
}
mean(samples_3_mean)
var(samples_3_mean)

samples_4 <- replicate(reps, rnorm(n4,2.5,1.118)) # 30 x 10000 sample matrix
samples_4_mean <- rep(NA, reps)
for (i in 1:reps) {
  samples_4_mean[i] <- mean(samples_4[,i])
}
mean(samples_4_mean)
var(samples_4_mean)

par(mfrow=c(3,1))

#hist(samples_1_mean,  freq=F, col="grey", xlab="", xlim=c(0, 5), breaks=100)
#par(new=T)
#plot(density(samples_1_mean), axes=F, main="", xlim=c(0, 5), lwd=2, col="blue")

hist(samples_2_mean,  freq=F, col="grey", xlab="", xlim=c(0, 5), breaks=100)
par(new=T)
plot(density(samples_2_mean), axes=F, main="", xlim=c(0, 5), lwd=2, col="blue")

hist(samples_3_mean,  freq=F, col="grey", xlab="", xlim=c(0, 5), breaks=100)
par(new=T)
plot(density(samples_3_mean), axes=F, main="", xlim=c(0, 5), lwd=2, col="blue")

hist(samples_4_mean,  freq=F, col="grey", xlab="", xlim=c(0, 5), breaks=100)
par(new=T)
plot(density(samples_4_mean), axes=F, main="", xlim=c(0, 5), lwd=2, col="blue")