# Montgomery 2.23
m <- 500
n <- 20
coeffs <- c()
mean.ests <- c()
conf.ints <- c()
mean.cints <- c()
for(i in seq(m)){
    x <- seq(1, 10, 0.5)
    eps <- rnorm(length(x), mean = 0, sd = 1)
    y <- 50 + 10 * x + eps
    fit <- lm(y ~ x)
    s <- summary(fit)
    # For a)
    coeffs <- rbind(coeffs, s$coefficients[,1])
    # For b)
    mean.est <- predict(fit, data.frame("x" = c(5)))
    mean.ests <- rbind(mean.ests, mean.est)
    # For c)
    conf.ints <- rbind(conf.ints, confint(fit)[2,])
    # For d)
    mean.cint <- predict(fit, data.frame("x" = c(5)),
                         interval="confidence",
                         level=0.95)
    mean.cints <- rbind(mean.cints, mean.cint)
}
# a) Plot parameter histograms
par(mfrow=c(2,2))
hist(coeffs[,1],20)
hist(coeffs[,2],20)

# b) Mean histogram
hist(mean.ests,20)

fit <- lm(y ~ x)
summary(fit)
tmp <- predict(fit, data.frame("x" = c(5)), interval="confidence", level=0.95)
tmp

# c) beta1 confidence interval. Cont the number of intervals that covers
# the true value beta1=10.
sum((conf.ints[,1] < 10) & (conf.ints[,2] > 10)) / m

# d) conficence interval for E[Y|x=5]. Cont the number of intervals that covers
# the true value 100.
sum(mean.cints[,2] < 100 & mean.cints[,3] > 100) / m