library(hdrm)

# Simple example of ridge and collinearity
# This is the same as
Ex1.2()
set.seed(11)

x1 <- rnorm(20)
x2 <- rnorm(20, mean = x1, sd = 0.01)
y <- rnorm(20, mean = 3 + x1 + x2)
lm(y ~ x1 + x2) |> coef()

ridge(y ~ x1 + x2, lambda = 0.1) |> coef()

# Sketch of proof
Fig1.3()

# Coverage of Bayes intervals
# 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=breheny::pal(2)[2], ylim=c(0,1))
abline(h=0.95, col="gray", lwd=3, lty=2)

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