# Montgomery Exercise 5.10
library(MPV)

# Bootstrap residuals, see Montgomery section 15.4
lm.fit <- lm(y ~ ., data = softdrink)
sum <- summary(lm.fit)
ms <- c(100, 500)
par(mfrow = c(2, 3))
for (m in ms) {
    coefs <- c()
    for (i in seq(m)) {
        n <- nrow(softdrink)
        indices <- sample(n, n, replace = TRUE)
        residual.boot <- sum$residuals[indices]
        y.boot <- lm.fit$fitted.values + residual.boot # New bootstrap samples, see (eq. 5.19)
        lm.fit.boot <- lm(y.boot ~ softdrink$x1 + softdrink$x2)
        coefs <- rbind(coefs, coef(lm.fit.boot))
    }

    param.sd.boot <- apply(coefs, 2, sd)
    print(param.sd.boot)
    # Also create confidence intervals accoring to the percentile method
    # presented in Section 15.4.2, Montgomery
    conf.ints <- c()
    for (k in seq( length(lm.fit$coefficients))) {
        hist(coefs[, k])
        quants <- quantile(coefs[, k], probs = c(0.025, 0.975))
        beta.est <- coef(lm.fit)[k]
        D1 <- beta.est - quants[1]
        D2 <- quants[2] - beta.est
        conf.ints <- rbind(conf.ints, c(beta.est - D2, beta.est + D1, beta.est))
    }

    colnames(conf.ints) <- c( names(quants), "beta est")
    rownames(conf.ints) <- names( coef(lm.fit))
    conf.ints
    confint(lm.fit)
}