set.seed(12)

n = 20

p=5
popVar=1

numReps = 1000

v=0
vstein=0
vME=0

sumMSE=0
sumStein=0
sumME=0

for(m in c(1:numReps)){  # varying PopulationMeans weights
PopulationMeans =rnorm(p, mean=0, sd=1)

SampleMeans<-c()
RepSampleMeans<-c()

	for(i in c(1:p)) {
	x<- rnorm(n, mean = PopulationMeans[i], sd = sqrt(popVar))
	SampleMeans<-c(SampleMeans, mean(x))
	Repx<- rnorm(n, mean = PopulationMeans[i], sd = sqrt(popVar))
	RepSampleMeans<-c(RepSampleMeans, mean(Repx))
	  }
 
  # Stein estimates
  adjust = (1- ((p-2)*popVar/n)/sum(SampleMeans*SampleMeans)  ) 
#	if(adjust <0){adjust =0}
  Steinestimates = adjust * SampleMeans
  
	# Morris-Efron estimates
	avg = mean(SampleMeans)
	MEadjust = (1- ((p-2)*popVar/n)/sum((SampleMeans-avg)*(SampleMeans-avg)) )  
	if(MEadjust<0){MEadjust=0}
	if(sum((SampleMeans-avg)*(SampleMeans-avg))   ==0){ MEadjust=0}
	MEestimates = MEadjust * (SampleMeans-avg) + avg
  
  # Compare estimates
  MSESamples = 0.0
  MSEStein = 0
  MSEME=0  # Morris -Efron estimator
  
  for(i in c(1:p)){
    MSESamples = MSESamples + (SampleMeans[i] - RepSampleMeans[i])^2	
    MSEStein = MSEStein + (Steinestimates[i] - RepSampleMeans[i])^2	
    MSEME = MSEME + (MEestimates[i] - RepSampleMeans[i])^2	
    
  #  cat(i, abs(SampleMeans[i] - PopulationMeans[i]), abs(Steinestimates[i] - PopulationMeans[i]), abs(MEestimates[i] - PopulationMeans[i]),"\n")
  }
  sumMSE = sumMSE+MSESamples
  sumStein = sumStein+MSEStein
  sumME = sumME + MSEME
  
  if(m<=-100){
  title<-sprintf("MSE=%f, MSE_JS=%f", MSESamples, MSEME)
#  title<-sprintf("MSE=%f", MSESamples)
  
  par(bg="lightblue")
  range = c(min(c(PopulationMeans, SampleMeans, MEestimates)), max(c(PopulationMeans, SampleMeans, MEestimates)))
 # range = c(min(c(PopulationMeans, SampleMeans)), max(c(PopulationMeans, SampleMeans)))
  plot(PopulationMeans, SampleMeans, main= title, xlim=range, ylim=range)
  points(PopulationMeans, MEestimates, pch=19)
  filename<-sprintf("MEplot%d.png",m)
dev.copy(png,filename)
dev.off()
}
#  cat(SampleMeans[1], SampleMeans[2], "\n")
#  cat(Steinestimates[1], Steinestimates[2], "\n")
  
  
  #cat(m, MSESamples, MSEStein, MSEME, "\n")
  
  if (MSEStein<MSESamples){ vstein= vstein+1}
  if (MSEME<MSESamples){ vME= vME+1}
    

}

cat("Average MSE\n")
cat("Samples = ", sumMSE/numReps, "\n")
cat("Stein = ", sumStein/numReps, "\n")
cat("ME = ", sumME/numReps, "\n")

cat("\nStein beats OLS ", (vstein/numReps), "of the time.\n")
cat("Morris-Efron beats OLS ", (vME/numReps), "of the time.\n")