# Setup
source('https://myweb.uiowa.edu/pbreheny/7110/f24/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")
par(mar=c(4,4,2,0))
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)