r - R 中最小二乘法的自制实现显示意外行为

Question

我正在构建一个示例，以图形方式展示最小二乘法的工作原理。我正在应用一种数值方法，我向 R 提供截距 (a) 和斜率 (b) 的多种可能值组合，然后计算所有可能组合的平方和 (SSE)。 a 和 b 与最低 SSE 的组合应该是最好的组合，但不知何故，我对 a 的估计与由 lm() 计算的实际值。最重要的是，我对 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 估计值，给我一个估计值a 大约在 -1428 - 除了让我的循环永远工作之外，如果我不是一个懒惰的 SCSS ，我可以通过使用 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 的估计值往往偏离目标。这仍然让我感到困惑:如果最小二乘法的分析方法确定了参数值的估计值，那么数值方法也应该如此。

不应该吗？

提前感谢您的反馈。

编辑

有人指出该问题是一个解决方案。我忽略了与 a 的每个值关联的 SSE 值并不独立于 b；最重要的是，SSE 的变化受 b 变化的影响比受 a 变化的影响更大(或者至少这是我的理解)。结果是斜率 b 的估计值的近似值可以覆盖截距 a 的估计值。

以下三个图表显示了 a、b 和 SSE 之间更广泛和更稀疏值范围的相关性:

possible.a.vals <- seq(-3000,1000, by=10)
possible.b.vals <- seq(-30,60, by=2)

Answer 1

@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