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

# An illustration of MLE problems due to 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')

# Variance explodes as p > n/2
Fig1.1()

# Histogram of beta estimates
out <- cache(Ex1.1(bar=FALSE), 'e1-1.pqt')
Fig1.2(out)

# Estimation error
b <- out$estimate
mean(b^2)
(mse_marg <- mean(1/replicate(10000, var(runif(25))*24)))
mean(b^2)/(mse_marg)

# Prediction error
n <- 25
mean(out$pe)/n
mean(out$pex)/n

# Inference problems
median(out$p)
range(out$p)
mean(out$cover)