## ——————————————
##
## ST4060 / ST6015 / ST6040
## R script – Wedn 29 Sept 2021 lecture
## Eric Wolsztynski
##
## ——————————————
## ——————————————
# (1) Example of a Monte Carlo experiment:
# evaluation of the distribution of the
# sample mean for a sample of size N=100
M = 10000 # number of MC repetitions
N = 100 # sample size
xbar = numeric(M) # just initialising a vector
# xbar will be used to store sample mean values
# as we calculate them
# repeat the same experiment M times…
for(m in 1:M){
x = rnorm(N, mean=10, sd=2)
xbar[m] = mean(x)
}
hist(xbar, breaks=41)
mean(xbar) # this should approximate E(X)=10
sd(xbar) # this should approximate SD(X)/sqrt(N)=2/10
2/sqrt(N)
## ——————————————
# (2) Example of a Monte Carlo experiment:
# evaluation of the distribution of the
# slope estimator in a linear regression model
# Monte carlo sample size:
M = 10000
# sample size:
N = 100
# true model parameter value :
thbar = 10
# (we pretend we don’t know it!)
# storing vectors:
est1 = numeric(M)
est2 = numeric(M)
for(m in 1:M){
# generate a random bunch of values for X:
x = runif(N, 2, 10)
# generate values for Y according to the
# theoretic model:
y = thbar*x + rnorm(N, mean=0, sd=2)
# fit the linear regresion model:
lmo = lm(y~x+0)
# store the model coefficient estimate:
est1[m] = coef(lmo)
}
# view the estimator’s sample distribution:
hist(est1, breaks=21)
# estimate the estimator’s distribution moments:
mean(xbar)
sd(xbar)
# Complete the above code to include another estimator,
# then compare their performance!