Calcul de l'algèbre linéaire pas à pas par régression des moindres carrés

22

En guise de préquelle à une question sur les modèles mixtes linéaires dans R, et à partager comme référence pour les aficionados de statistiques débutants / intermédiaires, j'ai décidé de publier en tant que "style Q&A" indépendant les étapes impliquées dans le calcul "manuel" du coefficients et valeurs prédites d'une régression linéaire simple.

L'exemple est avec l'ensemble de données intégré R mtcars, et serait configuré en miles par gallon consommé par un véhicule agissant comme variable indépendante, régressé sur le poids de la voiture (variable continue), et le nombre de cylindres en tant que facteur à trois niveaux (4, 6 ou 8) sans interactions.

EDIT: Si vous êtes intéressé par cette question, vous trouverez certainement une réponse détaillée et satisfaisante dans ce post de Matthew Drury hors CV .

Antoni Parellada
la source
Quand vous dites "calcul manuel", que recherchez-vous? Il est relativement simple de montrer une série d'étapes relativement simples pour obtenir des estimations de paramètres et ainsi de suite (via l'orthogonalisation de Gram-Schmidt, par exemple, ou par des opérateurs SWEEP), mais ce n'est pas ainsi que R effectue les calculs en interne; il (et la plupart des autres packages de statistiques) utilise la décomposition QR (discuté dans un certain nombre de publications sur le site - une recherche sur la décomposition QR révèle un certain nombre de publications, dont certaines pourraient vous procurer une valeur directe)
Glen_b -Reinstate Monica
Oui. Je crois que cela a été très bien traité dans la réponse de MD Je devrais probablement éditer mon message, en soulignant peut-être l'approche géométrique derrière ma réponse - espace de colonne, matrice de projection ...
Antoni Parellada
Oui! @Matthew Drury Voulez-vous que je supprime cette ligne dans l'OP ou que je mette à jour le lien?
Antoni Parellada
1
Je ne sais pas si vous avez ce lien, mais cela est étroitement lié, et j'aime vraiment la réponse de JM. stats.stackexchange.com/questions/1829/…
Haitao Du

Réponses:

51

Remarque : J'ai publié une version étendue de cette réponse sur mon site Web .

Envisageriez-vous de poster une réponse similaire avec le moteur R réel exposé?

Sûr! Dans le terrier du lapin, nous allons.

La première couche est lm, l'interface exposée au programmeur R. Vous pouvez regarder la source pour cela en tapant simplement lmsur la console R. La majorité (comme la majorité de la plupart des codes de niveau de production) est occupée à vérifier les entrées, à définir les attributs des objets et à lancer des erreurs; mais cette ligne ressort

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

lm.fitest une autre fonction R, vous pouvez l'appeler vous-même. Bien qu'il lmfonctionne de manière pratique avec des formules et un bloc de données, il lm.fitsouhaite des matrices, c'est donc un niveau d'abstraction supprimé. Vérification de la source lm.fit, travail plus chargé et ligne vraiment intéressante suivante

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

Nous arrivons maintenant quelque part. .Callest la façon dont R appelle le code C. Il y a une fonction C, C_Cdqrls quelque part dans la source R, et nous devons la trouver. Ça y est .

En regardant la fonction C, encore une fois, nous trouvons principalement la vérification des limites, le nettoyage des erreurs et le travail occupé. Mais cette ligne est différente

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

Alors maintenant, nous en sommes à notre troisième langage, R a appelé C qui appelle fortran. Voici le code fortran .

Le premier commentaire dit tout

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

(Fait intéressant, le nom de cette routine a changé à un moment donné, mais quelqu'un a oublié de mettre à jour le commentaire). Donc, nous sommes enfin au point où nous pouvons faire de l'algèbre linéaire et résoudre le système d'équations. C'est le genre de choses pour lesquelles fortran est vraiment bon, ce qui explique pourquoi nous avons traversé tant de couches pour arriver ici.

Le commentaire explique également ce que le code va faire

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 va donc résoudre le système en trouvant la décomposition QR

La première chose qui se produit, et de loin la plus importante, est

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

Cela appelle la fonction fortran dqrdc2sur notre matrice d'entrée x. Qu'est-ce que c'est ça?

 c     dqrfit uses the linpack routines dqrdc and dqrsl.

Nous l'avons donc finalement fait pour linpack . Linpack est une bibliothèque d'algèbre linéaire fortran qui existe depuis les années 70. La plus sérieuse des algèbres linéaires trouve finalement son chemin vers linpack. Dans notre cas, nous utilisons la fonction dqrdc2

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

C'est là que le travail réel est effectué. Il me faudrait une bonne journée complète pour comprendre ce que fait ce code, il est aussi bas que possible. Mais génériquement, nous avons une matrice et nous voulons la factoriser en un produit X = Q RQ est une matrice orthogonale et R est une matrice triangulaire supérieure. C'est une chose intelligente à faire, car une fois que vous avez Q et R, vous pouvez résoudre les équations linéaires de régressionXX=QRQRQR

XtXβ=XtY

très facilement. Effectivement

XtX=RtQtQR=RtR

de sorte que l'ensemble du système devient

RtRβ=RtQty

mais est triangulaire supérieur et a le même rang que X t X , donc tant que notre problème est bien posé, il est de rang complet, et nous pouvons tout aussi bien résoudre le système réduitRXtX

Rβ=Qty

Mais voici la chose impressionnante. est triangulaire supérieur, donc la dernière équation linéaire ici est juste , donc la résolution de β n est triviale. Vous pouvez ensuite remonter les lignes, une par une, et les remplacer par les β que vous connaissez déjà, en obtenant à chaque fois une équation linéaire simple à résoudre. Donc, une fois que vous avez Q et R , le tout s'effondre vers ce qu'on appelle la substitution à l'envers , ce qui est facile. Vous pouvez en savoir plus à ce sujet ici , où un petit exemple explicite est entièrement élaboré.Rconstant * beta_n = constantβnβQR

Matthew Drury
la source
4
C'était le court essai mathématique / de codage le plus amusant qu'on puisse imaginer. Je ne sais presque rien sur le codage, mais votre "visite" à travers les tripes d'une fonction R apparemment inoffensive était vraiment révélatrice. Excellente écriture! Puisque "gentiment" a fait l'affaire ... Pourriez-vous bien vouloir considérer celui-ci comme un défi connexe? :-)
Antoni Parellada
6
+1 Je n'avais jamais vu ça auparavant, joli résumé. Juste pour ajouter un peu d'informations au cas où @Antoni ne serait pas familier avec les transformations de Householder; c'est essentiellement une transformation linéaire qui vous permet de mettre à zéro une partie de la matrice R que vous essayez de réaliser sans nettoyer les parties que vous avez déjà traitées (tant que vous les parcourez dans le bon ordre), ce qui la rend idéale pour transformer des matrices en forme triangulaire supérieure (les rotations Givens font un travail similaire et sont peut-être plus faciles à visualiser, mais sont un peu plus lentes). Lorsque vous construisez R, vous devez en même temps construire Q
Glen_b -Reinstate Monica
2
Matthew (+1), je vous suggère de commencer ou de terminer votre message avec un lien vers votre rédaction beaucoup plus détaillée madrury.github.io/jekyll/update/2016/07/20/lm-in-R.html .
amibe dit Réintégrer Monica
3
-1 pour dégonfler et ne pas descendre au code machine.
S.Kolassa - Rétablir Monica
3
(Désolé, je plaisante ;-)
S. Kolassa - Réintègre Monica le
8

Les calculs réels étape par étape dans R sont magnifiquement décrits dans la réponse de Matthew Drury dans ce même fil. Dans cette réponse, je veux parcourir le processus de prouver à soi-même que les résultats dans R avec un exemple simple peuvent être atteints en suivant l'algèbre linéaire des projections sur l'espace de la colonne et le concept d'erreurs perpendiculaires (produit scalaire), illustré dans différents articles , et joliment expliqué par le Dr Strang dans Algèbre linéaire et ses applications , et facilement accessible ici .

β

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

À IDENTIQUE 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
la source
1
En général, en calcul numérique, je pense qu'il est préférable de résoudre l'équation linéaire au lieu de calculer la matrice inverse. Donc, je pense que beta = solve(t(X) %*% X, t(X) %*% y)c'est en pratique plus précis que solve(t(X) %*% X) %*% t(X) %*% y.
Matthew Drury
R ne le fait pas de cette façon - il utilise une décomposition QR. Si vous allez décrire l' algorithme utilisé, sur un ordinateur, je doute que quiconque utilise celui que vous montrez.
Rétablir Monica - G. Simpson
Pas après l'algorithme, juste pour essayer de comprendre les fondements de l'algèbre linéaire.
Antoni Parellada
@AntoniParellada Même dans ce cas, je trouve toujours la réflexion en termes d'équations linéaires plus éclairante dans de nombreuses situations.
Matthew Drury
1
Étant donné la relation périphérique de ce fil avec les objectifs de notre site, tout en jugeant utile d'illustrer l'utilisation de Rpour des calculs importants, je voudrais suggérer que vous envisagiez de le transformer en une contribution à notre blog.
whuber