Chapter 15. Shrinkage Estimation

Figure 15.1. MSE of James-Stein Estimator

x <- as.matrix(seq(0,5,by=0.01))
G <- length(x)

K <- c(4,6,12,48)
n <- length(K)
a <- as.matrix(1:500)

mse <- matrix(0,G,n)
for (ki in 1:n){    
  k <- K[ki]
  xk <- x*k
  mse[,ki] <- 1 - (((k-2)^2)/k)
  for (gi in 1:G){
    lambda <- xk[gi]
    lambda2a <- cumprod((lambda/2)/a)
    EQ <- (1/(k-2) + sum(lambda2a/(k-2+a*2)))*exp(-lambda/2)
    mse[gi,ki] <- 1 - EQ*(((k-2)^2)/k)
  }
}

plot(x,mse[,1],type="l",lty=1,xaxs="i",yaxs="i",ylim=c(0,1),ylab="MSE/K",xlab=expression(lambda/K),cex.lab=.8,cex.axis=.75)
lines(x,mse[,2],lty=2)
lines(x,mse[,3],lty=6)
lines(x,mse[,4],lty=5)
legend("bottomright",legend=c(expression(K==4),expression(K==6),expression(K==12),expression(K==48)),lty=c(1,2,6,5),cex=.75,bty="n") 

Figure 15.2. MSE of Positive-Part James-Stein Estimator

x <- as.matrix(seq(0,5,by=0.01))
G <- length(x)

K <- c(4,12)
n <- length(K)
a <- as.matrix(1:200)

mse <- matrix(0,G,n)
mse2 <- matrix(0,G,n)
for (ki in 1:n){    
  k <- K[ki]
  F <- pchisq(k-2,k-2+a*2)
  xk <- x*k
  for (gi in 1:G){
    lambda <- xk[gi]
    lambda2a <- cumprod((lambda/2)/a)
    EQ <- (1/(k-2) + sum(lambda2a/(k-2+a*2)))*exp(-lambda/2)
    ms <- 1 - EQ*(((k-2)^2)/k)
    mse[gi,ki] <- ms
    JK <- (pchisq(k-2,k-2)/(k-2) + sum(F*lambda2a/(k-2+a*2)))*exp(-lambda/2)
    mse2[gi,ki] <- ms - 2*pchisq(k-2,k,lambda) + pchisq(k-2,k+2,lambda) + lambda/k*pchisq(k-2,k+4,lambda) + ((k-2)^2)/k*JK
  }
}
 
plot(x,mse[,1],type="l",lty=1,xaxs="i",yaxs="i",ylim=c(0,1),ylab="MSE/K",xlab=expression(lambda/K),cex.lab=.8,cex.axis=.75)
lines(x,mse2[,1],lty=2)
lines(x,mse[,2],lty=6)
lines(x,mse2[,2],lty=5)
legend("bottomright",legend=c(expression(paste("JS,  ",K==4)),expression(paste("JS+, ",K==4)),
                              expression(paste("JS,  ",K==12)),expression(paste("JS+, ",K==12)))
       ,lty=c(1,2,6,5),cex=.75,bty="n")