library(hdrm)

# Simple example (Slide 14-15)
Ex1.2()

# MSE (Slide 20)
Fig1.3()

# Coverage (Slide 25)
# Based on assuming X'X/n=I, lam=1, sigma^2=n, in which case
# the posterior for beta is N(X'y/(2n), I/2)
b <- seq(-4, 4, len=99)
Coverage <- numeric(length(b))
for (i in 1:length(b)) {
  Coverage[i] <- 1 - pnorm(b[i]-4/sqrt(2)) - pnorm(b[i]+4/sqrt(2), 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)

# Expected coverage
f <- function(b) {
  dnorm(b) * (1 - pnorm(b-4/sqrt(2)) - pnorm(b+4/sqrt(2), lower.tail=FALSE))
}
integrate(f, -4, 4)