Самодельная реализация метода наименьших квадратов в R, показывающая неожиданное поведение - PullRequest
Новая
       14

Самодельная реализация метода наименьших квадратов в R, показывающая неожиданное поведение

0 голосов
/ 24 мая 2018

Я строю пример, чтобы графически показать, как работает метод наименьших квадратов.Я применяю численный подход, где я задаю R ряд комбинаций возможных значений точки пересечения (a) и наклона (b), а затем вычисляю сумму квадратов (SSE) для всех возможных комбинаций.Комбинации a и b с ассоциированным самым низким SSE должны быть лучшими, но почему-то мои оценки a всегда не соответствуют действительности по сравнению с реальной стоимостьювычисляется с помощью lm ().Кроме того, моя оценка a чувствительна к диапазону возможных значений данных для R - чем шире диапазон, тем больше оценка a отключена.

Вот мой пример.Я использую набор данных "Longley", построенный в R: 

    data(longley)
    plot(GNP ~ Employed, data = longley,
        xlab="% employed adults",
        ylab="Gross National Product (million $?)",
        main="Money money money"
        )

scatterplot of the data from the longley dataset

    # 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 - кроме того, чтобы заставить мой цикл работать вечно, что я мог бы решить, используя 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
        )

Changes in SSE for each combination of a and b

Это говорит о том, что корреляция между суммой квадратов и возможными значениями a является примерно плоской, что объясняет, почему оценки a имеют тенденцию быть не на должном уровне.Это все еще озадачивает: если аналитический подход к методу наименьших квадратов прибивает оценки значений параметров, то следует применять и численный подход.

Разве это не так?

Заранее спасибо за вашеобратная связь.

РЕДАКТИРОВАТЬ

Было отмечено, что проблема является разрешающей.Я упустил из виду, что значение SSE, связанное с каждым значением a , не зависит от 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)

Correlations between a, b, and SSE

1 Ответ

0 голосов
/ 24 мая 2018

@ Бен-Болкер прав.Не совсем правильно говорить, что ваша «оценка 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
...