gen.y<-function(theta){rnorm(length(theta),theta,1)}
gen.theta<-function(c){c*c(-1,-.5,0,.5,1)}
theta.hat.MLE<-function(y){y}
theta.hat.EB<-function(y){(1-(length(y)-1)/sum(y^2))*y}
loss<-function(theta,theta.hat){mean((theta-theta.hat)^2)}

n.c<-25
min.c<-0
max.c<-10
M<-10000

c.grid<-seq(min.c,max.c,length=n.c)
risk.F<-risk.EB<-rep(0,n.c)

for(j in 1:n.c){
   theta<-gen.theta(c.grid[j])
   loss.F<-loss.EB<-rep(0,M)
   for(s in 1:M){
      y<-gen.y(theta)
      loss.F[s]<-loss(theta,theta.hat.MLE(y))
      loss.EB[s]<-loss(theta,theta.hat.EB(y))
    }
   risk.F[j]<-mean(loss.F)
   risk.EB[j]<-mean(loss.EB)
}

plot(c.grid,risk.F,type="l",ylim=c(0.4,1.1),xlab="c",ylab="Frequentist Risk")
lines(c.grid,risk.EB,lty=2)
legend("bottomright",c("MLE","EB"),lty=1:2,inset=0.05)