source('http://myweb.uiowa.edu/pbreheny/7600/s16/notes/fun.R')

# Slide 14
set.seed(473)
x1 <- rnorm(20)
x2 <- rnorm(20, mean=x1, sd=.01)
y <- rnorm(20, mean=3+x1+x2)
lm(y~x1+x2)
require(MASS)
lm.ridge(y~x1+x2, lambda=1)

# Slide 20
f <- function(rho, lam, n) {
  G <- matrix(c(1, rho, rho,1), 2, 2)

  ## OLS
  bias.ols <- rep(0, 2)
  var.ols <- solve(G)/n
  mse.ols <- sum(diag(var.ols)) + crossprod(bias.ols)

  ## ridge
  W <- solve(G + lam/n*diag(2))
  var.ridge <- W %*% G %*% W / n
  bias.ridge <- lam/n * W %*% matrix(c(1,1), 2, 1)
  mse.ridge <- sum(diag(var.ridge)) + crossprod(bias.ridge)
  val <- c(mse.ols, mse.ridge, crossprod(bias.ridge), sum(diag(var.ridge)))
  names(val) <- c("mse.ols", "mse.ridge", "bias.ridge", "var.ridge")
  val
}
lam <- c(0, exp(seq(log(0.001), log(10), length=99)))
Y <- matrix(NA, 100, 4)
for (i in 1:100) {
  Y[i,] <- f(0.5, lam[i], 20)
}
matplot(lam, Y, type="l", col=c("gray50", pal(3)), lwd=3, lty=1, xlab=expression(lambda), las=1, ylab="", bty="n")
text(5.5, 0.10, "MSE", xpd=T)
text(9, 0.03, "Var", xpd=T)
text(9, 0.085, expression(Bias^2), xpd=T)

# Slide 25
b <- seq(-4, 4, len=99)
Coverage <- numeric(length(b))
for (i in 1:length(b)) {
  Coverage[i] <- 1 - pnorm(2*(b[i] - 2/sqrt(2)), b[i]) - pnorm(2*(b[i] + 2/sqrt(2)), b[i], lower.tail=FALSE)
}
plot(b, Coverage, type='l', las=1, bty='n', xlab=expression(beta), lwd=3, col=pal(2)[2], ylim=c(0,1))
abline(h=0.95, col="gray", lwd=3, lty=2)