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 .
źródło
Odpowiedzi:
Uwaga : opublikowałem rozszerzoną wersję tej odpowiedzi na mojej stronie internetowej .
Pewnie! W dół króliczej nory idziemy.
Pierwsza warstwa to
lm
interfejs udostępniony programatorowi R. Możesz zajrzeć do źródła tego po prostu pisząclm
na 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 wystajelm.fit
to kolejna funkcja R, którą możesz nazwać samodzielnie. Chociażlm
wygodnie współpracuje z formułami i ramką danych,lm.fit
chce macierzy, więc usunięto jeden poziom abstrakcji. Sprawdzanie źródłalm.fit
, więcej pracy i następujący naprawdę interesujący wierszTeraz gdzieś idziemy.
.Call
jest 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
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
(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
Fortran rozwiąże system, znajdując rozkład .QR
Pierwszą rzeczą, która się zdarza, i zdecydowanie najważniejszą, jest
To wywołuje funkcję fortran
dqrdc2
na naszej macierzy wejściowejx
. Co to jest?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
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 regresjiX 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ę jakR XtX , 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 βn β Q R
constant * beta_n = constant
źródło
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 .
lm
Identyczna:
coef(lm(mpg ~ wt + as.factor(cyl)-1))
.y_hat <- HAT %*% mpg
cyl <- as.factor(cyl); OLS <- lm(mpg ~ wt + cyl); predict(OLS)
:źródło
beta = solve(t(X) %*% X, t(X) %*% y)
jest w praktyce dokładniejszy niżsolve(t(X) %*% X) %*% t(X) %*% y
.R
ważnych obliczeń, chciałbym zasugerować, aby zastanowić się nad przekształceniem go w nasz blog.