Regresja najmniejszych kwadratów Obliczenia algebry liniowej krok po kroku

22

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 .

Antoni Parellada
źródło
Kiedy mówisz „obliczenia ręczne”, czego szukasz? Przedstawienie serii względnie prostych kroków w celu uzyskania oszacowania parametrów i tak dalej (na przykład przez ortogonalizację Gram-Schmidta lub przez operatorów SWEEP) jest stosunkowo proste, ale R nie wykonuje wewnętrznych obliczeń; to (i większość innych pakietów statystyk) wykorzystuje rozkład QR (omówiony w wielu postach na stronie - wyszukiwanie w rozkładzie QR
pokazuje
Tak. Uważam, że były to bardzo miłe adresy w odpowiedzi MD. Prawdopodobnie powinienem edytować swój post, być może podkreślając geometryczne podejście do mojej odpowiedzi - przestrzeń kolumny, matryca projekcji ...
Antoni Parellada
Tak! @Matthew Drury Czy chcesz, żebym usunął tę linię w OP lub zaktualizował link?
Antoni Parellada,
1
Nie jestem pewien, czy masz ten link, ale jest to ściśle powiązane i naprawdę podoba mi się odpowiedź JM. stats.stackexchange.com/questions/1829/…
Haitao Du

Odpowiedzi:

51

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 regresjiXX=QRQRQR

XtXβ=XtY

bardzo łatwo. W rzeczy samej

XtX=RtQtQR=RtR

cały system staje się

RtRβ=RtQty

ale jest górnym trójkątem i ma taką samą rangę jakRXtX , 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

Rβ=Qty

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.Rconstant * beta_n = constantβnβQR

Matthew Drury
źródło
4
To był najbardziej zabawny krótki esej matematyczno-kodowy, jaki można sobie wyobrazić. Prawie nic nie wiem o kodowaniu, ale twoja „wycieczka” po wnętrznościach pozornie nieszkodliwej funkcji R naprawdę otworzyła oczy. Doskonałe pisanie! Ponieważ „uprzejmie” załatwiło sprawę ... Czy mógłbyś łaskawie pod uwagę ten jeden jako powiązany wyzwanie? :-)
Antoni Parellada
6
+1 Nie widziałem tego wcześniej, fajne podsumowanie. Wystarczy dodać trochę informacji na wypadek, gdyby @Antoni nie znała transformacji Householder; jest to w zasadzie transformacja liniowa, która pozwala wyzerować jedną część macierzy R, którą próbujesz osiągnąć, bez pomijania części, z którymi już sobie poradziłeś (pod warunkiem, że przejdziesz ją we właściwej kolejności), dzięki czemu jest idealna do przekształcania matryc w górną trójkątną formę (rotacje Givens wykonują podobną pracę i być może łatwiej je zwizualizować, ale są nieco wolniejsze). Kiedy budujesz R, musisz jednocześnie zbudować Q
Glen_b
2
Matthew (+1), sugeruję, aby rozpocząć lub zakończyć swój post linkiem do znacznie bardziej szczegółowego zapisu madrury.github.io/jekyll/update/2016/07/20/lm-in-R.html .
ameba mówi Przywróć Monikę
3
-1 za strzepywanie i nie schodzenie do kodu maszynowego.
S. Kolassa - Przywróć Monikę
3
(Przepraszam, tylko żartowałem ;-)
S. Kolassa - Przywróć Monikę
8

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 .

β

mpg=intercept(cyl=4)+β1weight+D1intercept(cyl=6)+D2intercept(cyl=8)[]

D1D2X

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

[]lm

βProjMatrix=(XTX)1XT[ProjMatrix][y]=[RegrCoefs](XTX)1XTy=β

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)).

HatMatrix=X(XTX)1XT

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

y^X(XTX)1XTyy_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
Antoni Parellada
źródło
1
Ogólnie rzecz biorąc, uważam, że w obliczeniach numerycznych najlepiej jest rozwiązać równanie liniowe zamiast obliczyć macierz odwrotną. Myślę więc, że beta = solve(t(X) %*% X, t(X) %*% y)jest w praktyce dokładniejszy niż solve(t(X) %*% X) %*% t(X) %*% y.
Matthew Drury
R nie robi tego w ten sposób - wykorzystuje rozkład QR. Jeśli masz zamiar opisać zastosowany algorytm , wątpię, czy ktoś używa tego, który pokazujesz.
Przywróć Monikę - G. Simpson
Nie po algorytmie, po prostu próbując zrozumieć podstawy algebry liniowej.
Antoni Parellada
@AntoniParellada Nawet w takim przypadku myślenie w kategoriach równań liniowych jest bardziej pouczające w wielu sytuacjach.
Matthew Drury
1
Biorąc pod uwagę peryferyjny związek tego wątku z celami naszej witryny, a jednocześnie widząc wartość ilustrującą użycie Rważnych obliczeń, chciałbym zasugerować, aby zastanowić się nad przekształceniem go w nasz blog.
whuber