#ST810A                                        Spring '09
############################################################:
#  Simulation to compare sampling properties of three      #:
#  different estimators for the mean of a distribution     #:
#  based on an iid sample of size n:                       #:
#  -  sample mean                                          #:
#  -  20% trimmed mean                                     #:
#  -  median                                               #:
#  S = total number of simulated data sets                 #:
############################################################:






############################################################:
###    Define functions to be called in the simulation   ###:
############################################################:

#functions to generate S data sets of size n:

generate.gamma <- function(S,n,mu,sigma){
   a <- 1/(sigma^2)
   s <- mu/a
   return(matrix(rgamma(n*S,shape=a,scale=s),ncol=n,byrow=T)) 
}

generate.normal <- function(S,n,mu,sigma){
  return(matrix(rnorm(n*S,mu,sigma),ncol=n,byrow=T)) 
}

generate.t <- function(S,n,mu,df){
 return(matrix(mu + rt(n*S,df),ncol=n,byrow=T)) 
}

#  calculates summary statistics across the 1000 data sets 
simsum <- function(dat,trueval){
   MCmean <- mean(dat)
   MCbias <- MCmean-trueval
   MCrelbias <- MCbias/trueval
   MCstddev <- sd(dat)
   MCMSE <- mean((dat-trueval)^2)  
   
   return(list(MCmean=MCmean,MCbias=MCbias,MCrelbias=MCrelbias,
          MCstddev=MCstddev,MCMSE=MCMSE))
}
  
#  function to compute the 20% trimmed mean
trimmean <- function(Y){mean(Y,0.2)}







############################################################:
###                 Start the simulation!                ###:
############################################################:

set.seed(3)#  set the seed for the simulation

#  Generate Data
S <- 1000 
n <- 15
mu <- 1
sigma <- sqrt(5/3)

#out <- generate.gamma(S,n,mu,sigma) 
#out <- generate.t(S,n,mu,3) 
out <- generate.normal(S,n,mu,sigma) 

#get the sample means, trimmed means, and medians from each data set

outsampmean <- apply(out,1,mean)
outtrimmean <- apply(out,1,trimmean)
outmedian <- apply(out,1,median)

#  now we can summarize -- remember, mu is the true value of the mean
#  for the true distribution of the data

sum.mean<-simsum(outsampmean,mu)
sum.trimmean<-simsum(outtrimmean,mu)
sum.median<-simsum(outmedian,mu)



############################################################:
###                    Display results                   ###:
############################################################:

print(sum.mean$MCMSE)
print(sum.trimmean$MCMSE)
print(sum.median$MCMSE)

#Compute mse and its standard error the sample mean
mse<-(outsampmean-mu)^2
boxplot(mse)
print("Estimated MSE of the sample mean")
round(mean(mse),4)
print("Standard error of the estimated mse")
round(sd(mse)/sqrt(S),4)

#Compute mse and its standard error for the sample median
mse<-(outmedian-mu)^2
boxplot(mse)
print("Estimated MSE of the sample median")
round(mean(mse),4)
print("Standard error of the estimated mse")
round(sd(mse)/sqrt(S),4)