# Setup
source("https://pbreheny.org/7110/f26/notes/fun.R")
library(ggplot2)
theme_set(theme_minimal())
sim <- function(df, a, N = 1000, n = 100, d = 5) {
  X <- matrix(runif(n * d, -1, 1), n, d)
  X[1:2, 1] <- c(-a, a)
  s <- sqrt(solve(crossprod(X))[1, 1])
  b <- numeric(N)
  for (i in 1:N) {
    e <- rt(n, df) * sqrt((df - 2) / df)
    b[i] <- coef(lm(e ~ 0 + X))[1]
  }
  list(b = b, s = s, N = N, g = apply(X, 1, crossprod))
}

# Illustration of the two conditions (nu=3, a=5)
x <- seq(-4, 4, len = 101)
col <- c("#BEBEBE", "#6A5ACD")
plot(
  x,
  dnorm(x),
  type = "l",
  ylab = "",
  bty = "n",
  ylim = c(0, 0.5),
  col = col[1],
  xlab = "",
  lwd = 3,
  yaxt = "n"
)
lines(x, dt(x, df = 3), col = col[2], lwd = 3)
toplegend(legend = c("Normal", "t"), col = col, lwd = 3)

# Boxplot of x'x (rows of x)
set.seed(2)
res <- sim(3, 5)
ggbox(res$g, expression(X^T * X)) + ylab('')

# Results: nu=50, a=5
set.seed(1)
res <- sim(50, 5)
qqplot(
  qnorm(1:res$N / (res$N + 1)),
  res$b / res$s,
  pch = 19,
  col = 'gray50',
  bty = 'n',
  las = 1,
  xlab = 'Theoretical',
  ylab = 'Observed'
)
abline(0, 1, lwd = 2, lty = 2)

# Results: nu=3, a=1
set.seed(1)
res <- sim(3, 1)
qqplot(
  qnorm(1:res$N / (res$N + 1)),
  res$b / res$s,
  pch = 19,
  col = 'gray50',
  bty = 'n',
  las = 1,
  xlab = 'Theoretical',
  ylab = 'Observed'
)
abline(0, 1, lwd = 2, lty = 2)

# Results: nu=3, a=5
set.seed(1)
res <- sim(3, 5)
qqplot(
  qnorm(1:res$N / (res$N + 1)),
  res$b / res$s,
  pch = 19,
  col = 'gray50',
  bty = 'n',
  las = 1,
  xlab = 'Theoretical',
  ylab = 'Observed'
)
abline(0, 1, lwd = 2, lty = 2)

# Results: nu=3, a=5
set.seed(1)
res <- sim(3, 5, n = 1000)
qqplot(
  qnorm(1:res$N / (res$N + 1)),
  res$b / res$s,
  pch = 19,
  col = 'gray50',
  bty = 'n',
  las = 1,
  xlab = 'Theoretical',
  ylab = 'Observed'
)
abline(0, 1, lwd = 2, lty = 2)