n <- 40 iterate <- 10000 b0 <- 1 b1 <- 1 set.seed(10101) x <- rnorm(n) tstats <- rep(0,iterate) for (iter in 1:iterate) { u <- rchisq(n,1)-1 y <- b0+b1*x+u a <- coef(summary(lm(y~x))) tstats[iter] <- (a[2,1]-b1)/a[2,2] # (b1hat-b1)/se1 } hist(tstats,prob=T,ylim=c(0,1.1*dt(0,n-2))) lines(v <- -40:40/10,dt(v,n-2),lty="dotted")