Regression der kleinsten Quadrate Schrittweise lineare Algebra-Berechnung

22

Als Vorbemerkung zu einer Frage zu linear gemischten Modellen in R und als Referenz für Kenner der Statistik für Anfänger und Fortgeschrittene entschied ich mich, die Schritte, die bei der "manuellen" Berechnung der R-Werte erforderlich sind, als unabhängigen "Q & A-Stil" zu veröffentlichen Koeffizienten und vorhergesagte Werte einer einfachen linearen Regression.

Das Beispiel bezieht sich auf den in R eingebauten Datensatz mtcarsund wird als Meilen pro Gallone festgelegt, die von einem Fahrzeug verbraucht werden, das als unabhängige Variable fungiert, die über das Gewicht des Fahrzeugs (kontinuierliche Variable) und die Anzahl der Zylinder als a zurückgegangen ist Faktor mit drei Ebenen (4, 6 oder 8) ohne Wechselwirkungen.

BEARBEITEN: Wenn Sie an dieser Frage interessiert sind, werden Sie in diesem Beitrag von Matthew Drury außerhalb des Lebenslaufs definitiv eine ausführliche und zufriedenstellende Antwort finden .

Antoni Parellada
quelle
Wenn Sie "manuelle Berechnung" sagen, wonach suchen Sie dann? Es ist relativ einfach, eine Reihe relativ einfacher Schritte zu zeigen, um Parameterschätzungen usw. zu erhalten (z. B. durch Gram-Schmidt-Orthogonalisierung oder durch SWEEP-Operatoren), aber so führt R die Berechnungen nicht intern durch. Es (und die meisten anderen Statistikpakete) verwendet eine QR-Zerlegung (die in einer Reihe von Posts auf der Website beschrieben wird - eine Suche nach einer QR-Zerlegung
bringt
Ja. Ich glaube, das war eine sehr schöne Antwort in der Antwort von MD. Ich sollte wahrscheinlich meinen Beitrag überarbeiten und vielleicht den geometrischen Ansatz hinter meiner Antwort betonen - Spaltenraum, Projektionsmatrix ...
Antoni Parellada
Ja! @ Matthew Drury Soll ich diese Zeile im OP löschen oder den Link aktualisieren?
Antoni Parellada
1
Ich bin mir nicht sicher, ob Sie diesen Link haben, aber dieser ist eng verwandt, und ich liebe die Antwort von JM wirklich. stats.stackexchange.com/questions/1829/…
Haitao Du

Antworten:

51

Hinweis : Ich habe eine erweiterte Version dieser Antwort auf meiner Website veröffentlicht .

Würden Sie freundlicherweise in Betracht ziehen, eine ähnliche Antwort mit dem tatsächlichen R-Motor zu veröffentlichen?

Sicher! Das Kaninchenloch runter gehen wir.

Die erste Schicht ist lmdie Schnittstelle, die dem R-Programmierer ausgesetzt ist. Sie können die Quelle dafür anzeigen, indem Sie einfach lman der R-Konsole tippen. Der größte Teil davon (wie der größte Teil des Codes auf Produktionsebene) ist damit beschäftigt, Eingaben zu überprüfen, Objektattribute festzulegen und Fehler zu werfen. aber diese Linie ragt heraus

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

lm.fitist eine andere R-Funktion, die Sie selbst aufrufen können. Funktioniert zwar lmbequem mit Formeln und Datenrahmen, lm.fitmöchte aber Matrizen, so dass eine Abstraktionsebene entfernt wird. lm.fitSuchen Sie nach der Quelle , nach mehr Arbeit und nach der folgenden wirklich interessanten Zeile

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

Jetzt kommen wir voran. .Callist Rs Art, C-Code aufzurufen. Es gibt eine C-Funktion, C_Cdqrls, irgendwo in der R-Quelle, und wir müssen sie finden. Hier ist es .

Wenn wir uns die C-Funktion noch einmal ansehen, finden wir meistens eine Überprüfung der Grenzen, eine Fehlerbereinigung und viel Arbeit. Aber diese Linie ist anders

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

Jetzt sind wir in unserer dritten Sprache, R hat C gerufen, was fortran anruft. Hier ist der Fortran-Code .

Der erste Kommentar sagt alles

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

(Interessanterweise wurde der Name dieser Routine irgendwann geändert, aber jemand hat vergessen, den Kommentar zu aktualisieren). Damit sind wir endlich an dem Punkt angelangt, an dem wir eine lineare Algebra erstellen und das Gleichungssystem tatsächlich lösen können. Dies ist die Art von Dingen, in denen Fortran wirklich gut ist, was erklärt, warum wir so viele Schichten durchlaufen haben, um hierher zu gelangen.

Der Kommentar erklärt auch, was der Code tun wird

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

Also wird fortran das System lösen, indem es die Zerlegung findet.Q.R

Das Erste, was passiert, und bei weitem das Wichtigste, ist

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

Dies ruft die fortran -Funktion dqrdc2in unserer Eingabematrix auf x. Was ist das?

 c     dqrfit uses the linpack routines dqrdc and dqrsl.

Also haben wir es endlich geschafft, Linpack zu machen . Linpack ist eine Fortran-Bibliothek für lineare Algebra, die es seit den 70er Jahren gibt. Die schwerste lineare Algebra findet schließlich ihren Weg zum Linpack. In unserem Fall verwenden wir die Funktion dqrdc2

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

Hier wird die eigentliche Arbeit erledigt. Es würde einen guten ganzen Tag dauern, bis ich herausgefunden hätte, was dieser Code tut. Er ist so niedrig wie er kommt. Aber im Allgemeinen haben wir eine Matrix und wollen sie in ein Produkt zerlegen, wobei eine orthogonale Matrix und eine obere Dreiecksmatrix ist. Dies ist eine kluge Sache, denn sobald Sie und , können Sie die linearen Gleichungen für die Regression lösenX = Q R Q R Q RXX=Q.RQ.RQ.R

XtXβ=XtY.

sehr leicht. Tatsächlich

XtX=RtQ.tQ.R=RtR

so wird das ganze system

RtRβ=RtQ.ty

aber ist das obere Dreieck und hat den gleichen Rang wie , solange unser Problem gut gestellt ist, ist es der volle Rang, und wir können das reduzierte System genauso gut lösenX t XRXtX

Rβ=Q.ty

Aber hier ist die tolle Sache. ist das obere Dreieck, also ist die letzte lineare Gleichung gerade , also ist das Auflösen nach trivial. Sie können dann die Zeilen nacheinander durchlaufen und die bereits bekannten Werte ersetzen , wobei Sie jedes Mal eine einfache lineare Gleichung mit einer Variablen zur Lösung erhalten. Sobald Sie also und , kollabiert das Ganze zu einer so genannten Rückwärtssubstitution , was einfach ist. Sie können hier ausführlicher darüber lesen , wo ein explizites kleines Beispiel vollständig ausgearbeitet ist.β n β Q RRconstant * beta_n = constantβnβQ.R

Matthew Drury
quelle
4
Dies war der unterhaltsamste mathematisch / codierende Kurzaufsatz, den man sich vorstellen kann. Ich weiß so gut wie nichts über Codierung, aber Ihre "Tour" durch die Eingeweide einer scheinbar harmlosen R-Funktion war wirklich augenöffnend. Exzellentes Schreiben! Da „freundlich“ hat der Trick ... Könnten Sie freundlich betrachten diese ein als nahe stehende Herausforderung? :-)
Antoni Parellada
6
+1 Ich hatte das noch nie gesehen, nette Zusammenfassung. Nur um ein paar Informationen hinzuzufügen, falls @Antoni nicht mit den Transformationen von Hausbesitzern vertraut ist. Es handelt sich im Wesentlichen um eine lineare Transformation, mit der Sie einen Teil der R-Matrix, die Sie erreichen möchten, auf Null setzen können, ohne die bereits behandelten Teile zu zerstören (sofern Sie diese in der richtigen Reihenfolge durchlaufen). Dies macht sie ideal zur Umwandlung von Matrizen in obere Dreiecksform (Givens-Rotationen verrichten einen ähnlichen Job und sind möglicherweise einfacher zu visualisieren, jedoch etwas langsamer). Wie Sie R bauen, müssen Sie Q zugleich Konstrukt
Glen_b -Reinstate Monica
2
Matthew (+1), ich schlage vor, dass Sie Ihren Beitrag mit einem Link zu Ihrer ausführlicheren Beschreibung beginnen oder beenden . Madrury.github.io/jekyll/update/2016/07/20/lm-in-R.html .
Amöbe sagt Reinstate Monica
3
-1 für Hühnerei und nicht auf Maschinencode.
S. Kolassa - Reinstate Monica
3
(Sorry, nur ein Scherz ;-)
S. Kolassa - Reinstate Monica
8

Die tatsächlichen schrittweisen Berechnungen in R werden in der Antwort von Matthew Drury in diesem Thread sehr schön beschrieben. In dieser Antwort möchte ich durch den Prozess des Beweises gehen, dass die Ergebnisse in R mit einem einfachen Beispiel erreicht werden können, indem man der linearen Algebra von Projektionen auf den Spaltenraum und dem Konzept der senkrechten (Skalarprodukt-) Fehler folgt, die in verschiedenen Beiträgen dargestellt werden und erklärte schön Dr. Strang in Lineare Algebra und ihre Anwendungen und leicht zugänglich hier .

β

mpG=ichntercept(cyl=4)+β1weichGht+D1ichntercept(cyl=6)+D2ichntercept(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

βPrOjMeintrichx=(XTX)-1XT[PrOjMeintrichx][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

Identisch mit: coef(lm(mpg ~ wt + as.factor(cyl)-1)).

HeintMeintrichx=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
quelle
1
Im Allgemeinen glaube ich, dass es beim numerischen Rechnen am besten ist, die lineare Gleichung zu lösen, anstatt die inverse Matrix zu berechnen. Also, ich denke beta = solve(t(X) %*% X, t(X) %*% y)ist in der Praxis genauer als solve(t(X) %*% X) %*% t(X) %*% y.
Matthew Drury
R macht das nicht so - es verwendet eine QR-Zerlegung. Wenn Sie den verwendeten Algorithmus beschreiben wollen , bezweifle ich, dass auf einem Computer jemand den von Ihnen gezeigten verwendet.
Setzen Sie Monica - G. Simpson
Nicht nach dem Algorithmus, nur um die Grundlagen der linearen Algebra zu verstehen.
Antoni Parellada
@AntoniParellada Auch in diesem Fall finde ich das Denken in linearen Gleichungen in vielen Situationen aufschlussreicher.
Matthew Drury
1
In Anbetracht der peripheren Beziehung dieses Threads zu den Zielen unserer Website und der Tatsache, dass die Verwendung Rfür wichtige Berechnungen von Nutzen ist , möchte ich Sie bitten, ihn in einen Beitrag zu unserem Blog zu
whuber