我正在构建一个示例,以图形方式显示最小二乘法的工作原理。
我正在应用一种数值方法,在该方法中,我将截距(a)和斜率(b)的可能值的多个组合输入给R,然后为所有可能的组合计算平方和(SSE)。与最低SSE相关的a和b组合应该是最好的,但是与lm()计算的实际值相比,我对a的估计总是不合理。最重要的是,我对a的估计对给定R的可能值的范围很敏感-范围越宽,a的估计就越不正确。
这是我的例子。我正在使用内置于R中的数据集“longley”:
data(longley)
plot(GNP ~ Employed, data = longley,
xlab="% employed adults",
ylab="Gross National Product (million $?)",
main="Money money money"
)
# ranges of a and be where we think their true value lies:
possible.a.vals <- seq(-1431,-1430, by=0.01)
possible.b.vals <- seq(27,28.5, by=0.01)
# all possible combinations of a and b:
possible.ab <- expand.grid(possible.a.vals = possible.a.vals,
possible.b.vals = possible.b.vals
)
possible.ab.SSE <- as.data.frame(possible.ab)
head(possible.ab.SSE); tail(possible.ab.SSE)
possible.ab.SSE$SSE <- rep(NA, length.out = length(possible.ab.SSE[,1]))
for (i in 1:length(possible.ab.SSE[,1])){
predicted.GNP <- possible.ab.SSE$possible.a.vals[i] + possible.ab.SSE$possible.b.vals[i] * longley$Employed
possible.ab.SSE$SSE[i] <- sum((longley$GNP - predicted.GNP)^2)
}
possible.ab.SSE$possible.a.vals[which(possible.ab.SSE$SSE == min(possible.ab.SSE$SSE))]
possible.ab.SSE$possible.b.vals[which(possible.ab.SSE$SSE == min(possible.ab.SSE$SSE))]
# Estimate of a = -1430.73
# estimate of b = 27.84
# True values of a and b:
# a = -1430.48
# b = 27.84
我对b的估计是正确的,但a略有偏离。
此外,扩展a和b的可能值范围会产生对a的估计值,该估计值甚至比实际值还要远,这给了我大约-1428的估计值-除了使我的循环永远有效之外,我可以使用apply( )(如果我不是懒惰的驴子)。
# plot in 3d:
require(akima) # this helps interpolating the values of a,b, and SSE to create a surface
x= possible.ab.SSE$possible.a.vals
y= possible.ab.SSE$possible.b.vals
z=possible.ab.SSE$SSE
s=interp(x,y,z)
persp(x = s$x,
y = s$y,
z = s$z,
theta =50, phi = 10,
xlab="a", ylab="b", zlab="SSE",
box=T
)
这表明平方和与可能的a值之间的相关性大致是平坦的,这解释了为什么对趋势的估计往往没有意义。这仍然让我感到困惑:如果采用最小二乘法的分析方法确定参数值的估计值,那么数值方法也应该如此。
应该不是吗?
预先感谢您的反馈。
编辑
有人指出,这个问题是解决方案之一。我忽略了与a的每个值相关的SSE值并不独立于b;最重要的是,与b的变化相比,b的变化对SSE的变化影响更大(或者至少是我对此的理解)。结果是,斜率b的估计值的近似值可以抵消截距a的估计值。
以下是三个图表,显示了较大,较小范围的值在a,b和SSE之间的相关性:
possible.a.vals <- seq(-3000,1000, by=10)
possible.b.vals <- seq(-30,60, by=2)
最佳答案
@ ben-bolker是正确的。说您的“b的估计是正确的”并不是完全正确的。最小化示例中的SSE的值27.84
与OLS估计值27.83626
之间的差值会显着影响截距估计值。
data(longley)
# ranges of a and be where we think their true value lies:
possible.a.vals <- seq(-1431,-1430, by = 0.005)
possible.b.vals <- seq(27.5,28, by = 0.00001)
# all possible combinations of a and b:
possible.ab.SSE <- expand.grid(possible.a.vals = possible.a.vals,
possible.b.vals = possible.b.vals)
possible.ab.SSE <- as.matrix(possible.ab.SSE)
out <- tcrossprod(cbind(1, longley$Employed), possible.ab.SSE)
possible.ab.SSE <- as.data.frame(possible.ab.SSE)
possible.ab.SSE$SSE <- colSums((out - longley$GNP)^2)
possible.ab.SSE[order(possible.ab.SSE$SSE), ][1, ]
# possible.a.vals possible.b.vals SSE
# 6758127 -1430.48 27.83622 4834.891
coef(lm(GNP ~ Employed, data = longley))
# (Intercept) Employed
# -1430.48231 27.83626
关于r - R中的最小二乘方法的自制酿造实现,显示了意外的行为,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/50508425/