gistfile1.r

# http://mikelove.wordpress.com/2013/11/07/empirical-bayes/

# Stein's estimation rule and its competitors - an empirical Bayes approach
# B Efron, C Morris, Journal of the American Statistical, 1973
n <- 1000
sigma.means <- 5
means <- rnorm(n, 0, sigma.means)
# sigma.y <- 5
library(manipulate)
manipulate({
  y <- rnorm(n,means,sigma.y)
  A <- sum((y/sigma.y)^2)/(n - 2) - 1
  B <- 1/(1 + A)
  
  eb <- (1 - B) * y
  par(mfrow=c(1,2),mar=c(5,5,3,1))
  plot(means, y, xlim=c(-20,20),ylim=c(-20,20),
       main="y ~ N(mean, sigma.y)\nmean ~ N(0, sigma.mean=5)")
  abline(a=0,b=1,col='gray')
  points(means, eb, col="red")
  legend("topleft","James-Stein estimators",col="red",pch=1)
  s <- seq(from=0,to=1,length=100)
  par(mar=c(5,5,1,1))
  plot(s, sapply(s, function(b) sum((means - (1 - b)*y)^2)),
       type="l", xlab="possible values for B", 
       ylab="sum squared error")
  points(B, sum((means - eb)^2),col="red",pch=16)
  legend("top","James-Stein B",col="red",pch=16)
  
  print(B)
  cat(sprintf("naive mse %g\n", mean((means-y)^2)))
  cat(sprintf("stein mse %g\n", mean((means-eb)^2)))
}, sigma.y = slider(1,10))