source('https://myweb.uiowa.edu/pbreheny/7110/f20/notes/fun.R')

# Student's curves vs normal curves
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)
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=20), col=col[2], lwd=3)
toplegend(legend=c("Normal", "t"), col=col, lwd=3)

# Simulation
sim <- function(df, a, N=1000, n=100, d=5) {
  X <- matrix(runif(n*d), 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))
}

# df=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)

# df=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)

# df=3, a=5
set.seed(2)
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)

# box
ggbox(res$g, expression(X^T * X), '')

# n=1000, df=3
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)