### Check this code; no guarantee it is correct.
### To run this code, first save it in a file.
### Let XXXX denote the name of the file you saved it in.
### To run the code either
### (i) start R then type source("XXXX") or
### (ii) start R then cut and paste this code or
### (iii) do not start R and from a command line type:
### R CMD BATCH XXXX YYYY
### where YYYY is the name of the file in which you want to put the output.
library(MASS) #### loads the MASS library
simulate = function(n,p,mu,N){
### this function simulates data N times
### and computes the estimators
bayes = rep(0,N)
freq = rep(0,N)
Sigma = diag(rep(1,p))
for(i in 1:N){
x = mvrnorm(n,mu,Sigma)
xbar = apply(x,2,mean)
freq[i] = sum(xbar^2) - (p/n)
#### I am leaving this line for you to fill in bayes[i] =
}
return(list(freq=freq,bayes=bayes))
}
n = 100
p = 100
mu = rep(0,p)
N = 1000
out = simulate(n,p,mu,N)
pdf("plot.pdf")
par(mfrow=c(2,1))
hist(out$freq,nclass=25)
#hist(out$bayes,nclass=25)
abline(v=0,col="red",lwd=3)
dev.off()
cat("Plots are in plot.pdf \n")