复制最小二乘回归以检查估计和预测与真实情况的一致性
Replicate least squares regression to check consistency of estimation and prediction with truth
问题来了:在X = 4, 8, 12, 16, 20
时分别对Y
进行五次观察。真正的回归函数是E(y) = 20 + 4X
,ei
是独立的N(O, 25)
.
生成五个正态随机数,均值为 0,方差为 25。将这些随机数视为 X = 4,8, 12, 16, 20
处五个 Y
观测值的误差项,并计算 Y1
、Y2
、Y3
、Y4
和 Y5
。将直线拟合到五个案例时,获得最小二乘估计 bo
和 b1
。当 Xh = 10
时还计算 Yh
并在 Xh = 10
时获得 E(Yh)
的 95% 置信区间。我做了第 1 部分,但我需要帮助才能重复 200 次。
重复第(1)部分200次,每次生成新的随机数。
对 200 个估计值进行频率分布 b1
。计算 200 个估计值的均值和标准差 b1
。结果是否符合理论预期?
当 Xh = 10
时 E(Yh)
的 200 个置信区间中包含 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")
您的代码有误。您想要独立的 N(O, 25)
错误,但您将 sqrt(5)
作为标准错误传递给了 rnorm()
。应该是 5
.
我们首先将您的代码封装到一个函数中。此函数不接受输入,但 运行 实验一次,returns 回归系数 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 次模拟的结果。
第二列包含 200 次模拟下 b1
的估计值,我们计算它们的均值和标准误差:
mean(z[,2])
# [1] 3.976249
sd(z[,2])
# [1] 0.4263377
我们知道真实值为4
,可见我们的估计与真实值是一致的
最后,让我们检查 X = 10
处预测的 95% 置信区间。真值是20 + 4 * 10 = 60
,所以覆盖这个真值的置信区间比例是
mean(z[, "lwr"] < 60 & z[, "upr"] > 60)
## 0.95
正好0.95。
问题来了:在X = 4, 8, 12, 16, 20
时分别对Y
进行五次观察。真正的回归函数是E(y) = 20 + 4X
,ei
是独立的N(O, 25)
.
生成五个正态随机数,均值为 0,方差为 25。将这些随机数视为
X = 4,8, 12, 16, 20
处五个Y
观测值的误差项,并计算Y1
、Y2
、Y3
、Y4
和Y5
。将直线拟合到五个案例时,获得最小二乘估计bo
和b1
。当Xh = 10
时还计算Yh
并在Xh = 10
时获得E(Yh)
的 95% 置信区间。我做了第 1 部分,但我需要帮助才能重复 200 次。重复第(1)部分200次,每次生成新的随机数。
对 200 个估计值进行频率分布
b1
。计算 200 个估计值的均值和标准差b1
。结果是否符合理论预期?当
Xh = 10
时E(Yh)
的 200 个置信区间中包含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")
您的代码有误。您想要独立的 N(O, 25)
错误,但您将 sqrt(5)
作为标准错误传递给了 rnorm()
。应该是 5
.
我们首先将您的代码封装到一个函数中。此函数不接受输入,但 运行 实验一次,returns 回归系数 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 次模拟的结果。
第二列包含 200 次模拟下 b1
的估计值,我们计算它们的均值和标准误差:
mean(z[,2])
# [1] 3.976249
sd(z[,2])
# [1] 0.4263377
我们知道真实值为4
,可见我们的估计与真实值是一致的
最后,让我们检查 X = 10
处预测的 95% 置信区间。真值是20 + 4 * 10 = 60
,所以覆盖这个真值的置信区间比例是
mean(z[, "lwr"] < 60 & z[, "upr"] > 60)
## 0.95
正好0.95。