### 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")