Regresja najmniejszych kwadratów Obliczenia algebry liniowej krok po kroku

Jako prequel pytania o modele mieszane liniowo w R i jako odniesienie dla początkujących / średniozaawansowanych miłośników statystyki, postanowiłem opublikować jako niezależny styl „pytania i odpowiedzi” kroki związane z „ręcznym” obliczeniem współczynniki i przewidywane wartości prostej regresji liniowej.

Przykładem jest wbudowany zestaw danych R mtcars, który zostałby skonfigurowany jako mile na galon zużywane przez pojazd działający jako niezależna zmienna, regresowany względem masy samochodu (zmienna ciągła) i liczby cylindrów jako współczynnik z trzema poziomami (4, 6 lub 8) bez interakcji.

EDYCJA: Jeśli jesteś zainteresowany tym pytaniem, na pewno znajdziesz szczegółową i satysfakcjonującą odpowiedź w tym poście Matthew Drury poza CV .

Uwaga : opublikowałem rozszerzoną wersję tej odpowiedzi na mojej stronie internetowej .

Czy uprzejmie rozważysz opublikowanie podobnej odpowiedzi przy odsłoniętym silniku R?

Pewnie! W dół króliczej nory idziemy.

Pierwsza warstwa to lminterfejs udostępniony programatorowi R. Możesz zajrzeć do źródła tego po prostu pisząc lmna konsoli R. Większość z nich (podobnie jak większość kodu na poziomie produkcyjnym) zajmuje się sprawdzaniem danych wejściowych, ustawianiem atrybutów obiektu i zgłaszaniem błędów; ale ta linia wystaje

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

lm.fitto kolejna funkcja R, którą możesz nazwać samodzielnie. Chociaż lmwygodnie współpracuje z formułami i ramką danych, lm.fitchce macierzy, więc usunięto jeden poziom abstrakcji. Sprawdzanie źródła lm.fit, więcej pracy i następujący naprawdę interesujący wiersz

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

Teraz gdzieś idziemy. .Calljest sposobem wywoływania przez R. kodu C. Istnieje gdzieś funkcja C, C_Cdqrls w źródle R. I musimy ją znaleźć. Oto jest .

Patrząc na funkcję C, znowu znajdujemy głównie sprawdzanie granic, usuwanie błędów i zajęty pracę. Ale ta linia jest inna

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

Tak więc teraz jesteśmy w naszym trzecim języku, R nazywa C, który wzywa do fortranu. Oto kod fortran .

Pierwszy komentarz mówi wszystko

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

(co ciekawe, wygląda na to, że nazwa tej procedury została w pewnym momencie zmieniona, ale ktoś zapomniał zaktualizować komentarz). W końcu jesteśmy w punkcie, w którym możemy wykonać algebrę liniową i rozwiązać układ równań. To jest coś, w czym fortran jest naprawdę dobry, co wyjaśnia, dlaczego przeszliśmy przez tak wiele warstw, aby się tu dostać.

Komentarz wyjaśnia również, co zrobi kod

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

Fortran rozwiąże system, znajdując rozkład .$QR$

Pierwszą rzeczą, która się zdarza, i zdecydowanie najważniejszą, jest

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

To wywołuje funkcję fortran dqrdc2na naszej macierzy wejściowej x. Co to jest?

 c     dqrfit uses the linpack routines dqrdc and dqrsl.

W końcu dotarliśmy do Linpacka . Linpack jest fortranową biblioteką algebry liniowej, która istnieje od lat 70. Najpoważniejsza algebra liniowa ostatecznie trafia do linpack. W naszym przypadku korzystamy z funkcji dqrdc2

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

Właśnie tam wykonywana jest praca. Zajmie mi to cały dzień, żeby dowiedzieć się, co robi ten kod. Jest tak niski, jak to tylko możliwe. Ale generalnie mamy macierz i chcemy ją uwzględnić w produkcie gdzie jest macierzą ortogonalną, a jest górną trójkątną macierzą. Jest to mądra rzecz do zrobienia, ponieważ gdy masz i , możesz rozwiązać równania liniowe dla regresji$X$$X = QR$$Q$$R$$Q$$R$

bardzo łatwo. W rzeczy samej

cały system staje się

ale jest górnym trójkątem i ma taką samą rangę jak$R$$X^t X$ , więc dopóki nasz problem jest dobrze postawiony, ma pełną rangę i równie dobrze możemy po prostu rozwiązać zredukowany system

Ale oto niesamowita rzecz. jest górnym trójkątem, więc ostatnie równanie liniowe jest tutaj , więc rozwiązywanie β n jest banalne. Następnie możesz wchodzić po rzędach, jeden po drugim, i zamieniać β, które już znasz, za każdym razem uzyskując proste równanie liniowe o zmiennej zmiennej do rozwiązania. Tak więc, kiedy masz Q i R , cała sprawa zapada się do tak zwanej substytucji wstecznej , co jest łatwe. Możesz przeczytać o tym bardziej szczegółowo tutaj , gdzie w pełni opracowano wyraźny mały przykład.$R$constant * beta_n = constant$\beta_n$$\beta$$Q$$R$

Rzeczywiste obliczenia krok po kroku w R są pięknie opisane w odpowiedzi Matthew Drury'ego w tym samym wątku. W tej odpowiedzi chcę przejść przez proces udowodnienia sobie, że wyniki w R za pomocą prostego przykładu można uzyskać po algebrze liniowej rzutów na przestrzeń kolumny i koncepcji błędów prostopadłych (iloczynu iloczynowego), zilustrowanych w różnych postach i dobrze wyjaśnione przez dr Strang w Algebrze liniowej i jej zastosowaniach , i łatwo dostępne tutaj .

$\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

Identyczna: 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