Репликация регрессии наименьших квадратов для проверки согласованности оценки и предсказания с истиной

Question

В этом проблема: Пять наблюдений на Y принимаются, если X = 4, 8, 12, 16, 20, соответственно. Истинная функция регрессии - E(y) = 20 + 4X, а ei являются независимыми N(O, 25).

1. Сформировать пять нормальных случайных чисел, со средним значением 0 и дисперсией 25. Рассмотрим эти случайные числа в качестве условий ошибок для пяти Y наблюдений на X = 4,8, 12, 16, 20 и рассчитать Y1, Y2, Y3, Y4 и Y5. Получите оценки наименьших квадратов bo и b1 при установке прямой линии в пять корпусов. Также рассчитать Yh, когда Xh = 10 и получить 95-процентный доверительный интервал для E(Yh), когда Xh = 10. Я сделал часть 1, но мне нужна помощь, чтобы повторить ее 200 раз.

2. Повторите часть (1) 200 раз, генерируя новые случайные числа каждый раз.

3. Сделать частотное распределение 200 оценок b1. Рассчитайте среднее и стандартное отклонение 200 оценок b1. Согласуются ли результаты с теоретическими ожиданиями?

4. Какая доля 200 доверительных интервалов для E(Yh), когда Xh = 10 включает E(Yh)? Согласуется ли этот результат с теоретическими ожиданиями?

Вот мой код до сих пор, я озадачен о том, как повторить часть 1 200 раз:

X <- matrix(c(4, 8, 12, 16, 20), nrow = 5, ncol = 1) 
e <- matrix(c(rnorm(5,0,sqrt(5))), nrow = 5, ncol = 1) 
Y <- 20 + 4 * X + e 
mydata <- data.frame(cbind(Y=Y, X=X, e=e)) 
names(mydata) <- c("Y","X","e") 
reg<-lm(Y ~ X, data = mydata) 
predict(reg, newdata = data.frame(X=10), interval="confidence") 

Answer 1

Существует ошибка в коде. Вам нужны независимые ошибки N(O, 25), но вы приняли sqrt(5) в качестве стандартной ошибки до rnorm(). Это должно быть 5.

Сначала мы завершим ваш код в функцию. Эта функция не принимает никаких данных, но запускает эксперимент один раз и возвращает коэффициенты регрессии b0, b1 и прогноз fit, lwr, upr в формате с именем.

sim <- function() { 
    x <- c(4, 8, 12, 16, 20) 
    y <- 20 + 4 * x + rnorm(5,0,5) 
    fit <- lm(y ~ x) 
    pred <- predict(fit, data.frame(x = 10), interval = "confidence") 
    pred <- setNames(c(pred), dimnames(pred)[[2]]) 
    ## return simulation result 
    c(coef(fit), pred) 
    } 

Например, давайте попробуем

set.seed(2016) 
sim() 
#(Intercept)   x   fit   lwr   upr 
# 24.222348 3.442742 58.649773 47.522309 69.777236 

Теперь мы используем replicate повторить такой эксперимент в 200 раз.

set.seed(0) 
z <- t(replicate(200, sim())) 

head(z) 
#  (Intercept)  x  fit  lwr  upr 
#[1,] 24.100535 3.987755 63.97808 57.61262 70.34354 
#[2,] 6.417639 5.101501 57.43265 52.44263 62.42267 
#[3,] 20.652355 3.797991 58.63227 52.74861 64.51593 
#[4,] 20.349829 3.816426 58.51409 52.59115 64.43702 
#[5,] 19.891873 4.095140 60.84327 57.49911 64.18742 
#[6,] 24.586749 3.589483 60.48158 53.64574 67.31743 

Будет создано 200 строк, для результатов 200 симуляций. Второй столбец содержит оценки для b1 при 200 симуляций, вычислим их среднее значение и стандартное отклонение:

mean(z[,2]) 
# [1] 3.976249 

sd(z[,2]) 
# [1] 0.4263377 

Мы знаем, что истинное значение 4, и очевидно, что наша оценка согласуется с истинными значениями.

Наконец, давайте проверим с доверительным интервалом 95% для прогноза на X = 10. Истинное значение 20 + 4 * 10 = 60, поэтому доля доверительного интервала, который охватывает этот истинный Вейл является

mean(z[, "lwr"] < 60 & z[, "upr"] > 60) 
## 0.95 

который точно 0,95.