Пошаговая линейная алгебраическая вычисление регрессии наименьших квадратов

В качестве приквела к вопросу о линейно-смешанных моделях в R и для использования в качестве справочного материала для любителей статистики начального и среднего уровня, я решил опубликовать в качестве независимого «Q & A-style» шаги, связанные с «ручным» вычислением Коэффициенты и прогнозные значения простой линейной регрессии.

В качестве примера используется встроенный набор данных R mtcars, и он будет установлен в виде миль на галлон, потребляемых транспортным средством, действующим в качестве независимой переменной, с регрессией на вес автомобиля (непрерывная переменная) и количество цилиндров в виде фактор с тремя уровнями (4, 6 или 8) без взаимодействия.

РЕДАКТИРОВАТЬ: Если вы заинтересованы в этом вопросе, вы обязательно найдете подробный и удовлетворительный ответ в этом посте Мэтью Друри за пределами резюме .

Примечание . Я разместил расширенную версию этого ответа на своем веб-сайте .

Не могли бы вы опубликовать аналогичный ответ с действительным R-движком?

Конечно! Вниз по кроличьей норе мы идем.

Первый уровень - lmинтерфейс, предоставляемый программисту R. Вы можете посмотреть на источник для этого, просто набрав lmна консоли R. Большая часть этого (как и большая часть кода производственного уровня) занята проверкой входных данных, установкой атрибутов объекта и выдачей ошибок; но эта линия торчит

lm.fit(x, y, offset = offset, singular.ok = singular.ok, 
                ...)

lm.fitэто еще одна функция R, вы можете назвать ее самостоятельно. Хотя lmудобно работать с формулами и фреймом данных, lm.fitхочет матрицы, так что это один уровень абстракции удаляется. Проверка источника lm.fit, более занятая работа и следующая действительно интересная строка

z <- .Call(C_Cdqrls, x, y, tol, FALSE)

Теперь мы куда-то добираемся. .Callспособ вызова R в C-код. Где-то есть функция C, C_Cdqrls в источнике R, и нам нужно ее найти. Здесь .

Глядя на функцию C, мы снова обнаруживаем в основном проверку границ, очистку от ошибок и занятую работу. Но эта линия отличается

F77_CALL(dqrls)(REAL(qr), &n, &p, REAL(y), &ny, &rtol,
        REAL(coefficients), REAL(residuals), REAL(effects),
        &rank, INTEGER(pivot), REAL(qraux), work);

Итак, теперь мы находимся на нашем третьем языке, R назвал C, который вызывает в fortran. Вот код Фортрана .

Первый комментарий говорит обо всем

c     dqrfit is a subroutine to compute least squares solutions
c     to the system
c
c     (1)               x * b = y

(интересно, похоже, что название этой подпрограммы было изменено в какой-то момент, но кто-то забыл обновить комментарий). Таким образом, мы, наконец, находимся в точке, где мы можем сделать некоторую линейную алгебру и действительно решить систему уравнений. Это то, в чем действительно хорошо разбирается фортран, и это объясняет, почему мы прошли через так много слоев, чтобы попасть сюда.

Комментарий также объясняет, что будет делать код

c     on return
c
c        x      contains the output array from dqrdc2.
c               namely the qr decomposition of x stored in
c               compact form.

Итак, Фортран собирается решить систему, найдя разложение $QR$

Первое, что происходит, и, безусловно, самое важное, это

call dqrdc2(x,n,n,p,tol,k,qraux,jpvt,work)

Это вызывает функцию fortran dqrdc2в нашей входной матрице x. Что это?

 c     dqrfit uses the linpack routines dqrdc and dqrsl.

Итак, мы наконец-то добрались до linpack . Linpack - это библиотека линейной алгебры Фортрана, которая существует с 70-х годов. Наиболее серьезная линейная алгебра в конечном итоге находит свой путь в linpack. В нашем случае мы используем функцию dqrdc2

c     dqrdc2 uses householder transformations to compute the qr
c     factorization of an n by p matrix x.

Это где фактическая работа сделана. Мне потребовался бы хороший полный день, чтобы понять, что делает этот код, он так же низок, как и они. Но, как правило, у нас есть матрица и мы хотим разложить ее на произведение X = Q R, где Q - ортогональная матрица, а R - верхняя треугольная матрица. Это умная вещь, потому что, получив Q и R, вы можете решить линейные уравнения для регрессии.$X$$X = QR$$Q$$R$$Q$$R$

очень легко. Верно

так что вся система становится

но имеет верхнюю треугольную форму и имеет тот же ранг, что и X t X , поэтому, пока наша задача хорошо поставлена, она является полным рангом, и мы можем также просто решить приведенную систему$R$$X^t X$

Но вот потрясающая вещь. верхняя треугольная, поэтому последнее линейное уравнение здесь просто , поэтому решение для β n тривиально. Затем вы можете пройти вверх по строкам, одну за другой, и подставить в β, которые вы уже знаете, каждый раз, получая простое линейное уравнение с одной переменной для решения. Итак, когда у вас есть Q и R , все это сводится к тому, что называется обратной заменой , что легко. Вы можете прочитать об этом более подробно здесь , где явный маленький пример полностью проработан.$R$constant * beta_n = constant$\beta_n$$\beta$$Q$$R$

Фактические пошаговые вычисления в R прекрасно описаны в ответе Мэтью Друри в этой же теме. В этом ответе я хочу пройти процесс доказательства себя, что результаты в R на простом примере могут быть достигнуты, следуя линейной алгебре проекций на пространство столбцов и концепции ошибок перпендикулярного (точечного произведения), проиллюстрированной в разных постах. и хорошо объяснено доктором Странгом в « Линейной алгебре и ее приложениях» и легко доступно здесь .

$\small \beta$

$\small D1$$\small D2$$\small X$

attach(mtcars)    
x1 <- wt

    x2 <- cyl; x2[x2==4] <- 1; x2[!x2==1] <-0

    x3 <- cyl; x3[x3==6] <- 1; x3[!x3==1] <-0

    x4 <- cyl; x4[x4==8] <- 1; x4[!x4==1] <-0

    X <- cbind(x1, x2, x3, x4)
    colnames(X) <-c('wt','4cyl', '6cyl', '8cyl')

head(X)
        wt 4cyl 6cyl 8cyl
[1,] 2.620    0    1    0
[2,] 2.875    0    1    0
[3,] 2.320    1    0    0
[4,] 3.215    0    1    0
[5,] 3.440    0    0    1
[6,] 3.460    0    1    0

$\small [*]$lm

$\small\beta$$\small ProjMatrix = \small (X^{T}X)^{-1}X^{T}$$\small [ProjMatrix] \,[y]\, =\, [RegrCoef's]$$\small (X^{T}X)^{-1}X^{T}\,y = \beta$

X_tr_X_inv <- solve(t(X) %*% X)    
Proj_M <- X_tr_X_inv %*% t(X)
Proj_M %*% mpg

          [,1]
wt   -3.205613
4cyl 33.990794
6cyl 29.735212
8cyl 27.919934

Идентичный: coef(lm(mpg ~ wt + as.factor(cyl)-1)).

$\small Hat Matrix = \small X(X^{T}X)^{-1}X^{T}$

HAT <- X %*% X_tr_X_inv %*% t(X)

$\hat{y}$$\small X(X^{T}X)^{-1}X^{T}\,y$y_hat <- HAT %*% mpg

cyl <- as.factor(cyl); OLS <- lm(mpg ~ wt + cyl); predict(OLS):

y_hat <- as.numeric(y_hat)
predicted <- as.numeric(predict(OLS))
all.equal(y_hat,predicted)
[1] TRUE