# remotes::install_github("pbreheny/hdrm")
library(hdrm)

# Slide 11: Flat likelihood
library(rgl)
a <- rnorm(100)
X <- cbind(a+rnorm(100, sd=0.3), a+rnorm(100, sd=0.3))
cor(X)
b <- seq(-2, 2, len=99)
Z <- matrix(NA, 99, 99)
for (i in 1:99) {
  for (j in 1:99) {
    Z[i,j] <- -crossprod(X %*% c(b[i], b[j]))
  }
}
persp3d(b, b, Z, col='gray')

# Slide 12
Fig1.1()

# Slide 19
out <- Ex1.1()
Fig1.2(out)

# Slide 21
n <- 25
b <- out$results[,"Estimate"]
mean(b^2)
12/n
mean(b^2)/(12/n)

# Slide 23
mean(out$pred^2)

# Slide 24
median(out$results[,"Pr(>|t|)"])
range(out$results[,"Pr(>|t|)"])
mean(out$cover)