## ——————————————
##
## ST4060 / ST6015 / ST6040
## R script – Wedn 6 Oct 2021 lecture
## Eric Wolsztynski
##
## ——————————————
## ——————————————
# (1) Monte Carlo simulation: illustrating the CLT for the sample mean
set.seed(4060)
M = 100
mu = 10
sigma = 3
# different sample sizes:
Ns = c(20, 30, 50, 100, 200, 500)
# to store all sample mean estimates:
means = matrix(NA, nrow=M, ncol=length(Ns))
for(i in 1:length(Ns)){ # i.e. for each sample size:
N = Ns[i]
X = matrix( rnorm(N*M, mean=mu, sd=sigma), nrow=M, ncol=N)
means[,i] = apply(X, 1, mean)
}
# now inspect the distributions with respect to sample size N:
boxplot(means, names=Ns, xlab=”sample size N”)
abline(h=mu, lwd=4, col=’red’)
# equivalent but different coding approach:
means = matrix(NA, nrow=M, ncol=length(Ns))
for(i in 1:length(Ns)){ # i.e. for each sample size:
N = Ns[i]
for(m in 1:M){
x = rnorm(N, mean=mu, sd=sigma)
means[m, i] = mean(x)
}
}
# now inspect the distributions with respect to sample size N:
boxplot(means, names=Ns, xlab=”sample size N”)
abline(h=mu, lwd=4, col=2)
# We can add the theoretic values of the confidence intervals around
# the sample mean estimates, for each of these distributions,
# since we know the theoretic distribution of the sample mean:
lower.bounds = mu – 1.96*sigma/sqrt(Ns)
upper.bounds = mu + 1.96*sigma/sqrt(Ns)
lower.bounds = apply(means,2,mean) – 1.96*sigma/sqrt(Ns)
upper.bounds = apply(means,2,mean) + 1.96*sigma/sqrt(Ns)
lines(lower.bounds, col=’navy’, lwd=2)
lines(upper.bounds, col=’navy’, lwd=2)