Régression linéaire

Objectifs du chapitre

À la fin de ce chapitre, vous serez capable de:

poser la régression linéaire comme une minimisation du risque empirique pour la perte carrée, et dériver la droite des moindres carrés en annulant les dérivées partielles;
écrire le problème sous forme matricielle, établir les équations normales $\tilde{X}^{⊤} \tilde{X} w = \tilde{X}^{⊤} y$ et les lire comme une projection orthogonale de $y$ sur l'espace engendré par les colonnes de $\tilde{X}$ ;
calculer et interpréter le coefficient de détermination $R^{2}$ , et lire les coefficients d'un modèle avec leurs unités, sans les prendre pour des effets causaux;
étendre le modèle à plusieurs caractéristiques (régression multiple) et à des caractéristiques polynomiales, en comprenant pourquoi le modèle reste linéaire en ses paramètres;
démontrer que, sous un bruit gaussien, les moindres carrés coïncident avec le maximum de vraisemblance;
reconnaître une colinéarité qui rend $\tilde{X}^{⊤} \tilde{X}$ singulière ou presque, et en mesurer les conséquences sur les coefficients.

Le problème: la meilleure droite

Au chapitre 1, la règle de l'agence, $f (x) = 500 + 20 x$ , prédisait les loyers du jeu A avec une erreur quadratique moyenne de $20550$ CHF². Elle avait été fixée d'avance, par l'expérience d'un agent immobilier. La question de ce chapitre est simple à poser: parmi toutes les droites $f (x) = w_{0} + w_{1} x$ , laquelle a la plus petite erreur quadratique moyenne sur les huit appartements? Et peut-on la calculer exactement, plutôt que de la chercher à tâtons?

La réponse est oui, et la méthode, vieille de plus de deux siècles, reste l'une des plus utilisées de toute la statistique appliquée. Elle est aussi la porte d'entrée de presque tout ce cours: la régression logistique (chapitre 4), la régularisation (chapitre 6), les méthodes à noyau (chapitre 9) et la dernière couche d'un réseau de neurones (chapitre 10) sont des variations sur le modèle linéaire. Il vaut la peine de le comprendre sous tous ses angles — analytique, matriciel, géométrique et probabiliste — car chacun de ces angles éclaire une généralisation différente.

Rappelons les données. Ce sont des données fictives, construites pour le cours: huit appartements lausannois, leur surface $x$ en m² et leur loyer mensuel $y$ en CHF, charges comprises.

Logement	1	2	3	4	5	6	7	8
Surface $x$ (m²)	30	45	50	62	70	85	95	110
Loyer $y$ (CHF)	1 200	1 550	1 380	1 840	1 750	2 350	2 250	2 900

En Python, comme au chapitre 1:

# Jeu de données A: (surface en m², loyer en CHF). Données fictives.
loyers = [(30, 1200), (45, 1550), (50, 1380), (62, 1840),
          (70, 1750), (85, 2350), (95, 2250), (110, 2900)]

C'est exactement le principe de minimisation du risque empirique du chapitre 1, appliqué au modèle affine et à la perte carrée: $J (w)$ est le risque empirique $\hat{R} (f_{w})$ . Diviser par $n$ ne change pas le minimiseur; on garde $SCR$ pour les calculs à la main, où il évite des fractions, et $J$ quand on compare des jeux de tailles différentes.

Pourquoi le carré des résidus? Le chapitre 1 a donné les raisons de fond: la perte carrée pénalise fortement les grandes erreurs, ce qui convient quand une grosse erreur coûte cher, et sa minimisation sur une constante donne la moyenne. Il y en a deux autres, plus techniques, que ce chapitre va rendre précises. La première est que $SCR$ est un polynôme du second degré en $(w_{0}, w_{1})$ : ses dérivées sont affines, et annuler des fonctions affines revient à résoudre un système linéaire, ce que l'on sait faire exactement. La seconde est probabiliste: si le bruit qui sépare le loyer de la droite est gaussien, les moindres carrés sont l'estimateur du maximum de vraisemblance (théorème 2.5). Avec la perte absolue, aucune de ces deux propriétés ne tient — il n'y a pas de formule close, et l'on doit recourir à la programmation linéaire.

Question 2.1

Pourquoi ne choisit-on pas plutôt la droite qui rend la somme des résidus $\sum_{i} e_{i}$ aussi proche que possible de zéro?

Les moindres carrés en dimension un

La fonction $SCR (w_{0}, w_{1})$ est un polynôme du second degré en deux variables. Si elle a un minimum, ses deux dérivées partielles s'y annulent (annexe A rappelle ce que sont les dérivées partielles et le gradient). Calculons-les.

Démonstration. Conditions du premier ordre. La dérivée d'un carré $(y_{i} - w_{0} - w_{1} x_{i})^{2}$ par rapport à vaut , et par rapport à elle vaut . En sommant,

\frac{\partial SCR}{\partial w _{0}} = - 2 i = 1 \sum n (y_{i} - w_{0} - w_{1} x_{i}), \frac{\partial SCR}{\partial w _{1}} = - 2 i = 1 \sum n x_{i} (y_{i} - w_{0} - w_{1} x_{i}) .

Annuler ces deux dérivées donne un système de deux équations linéaires en $(w_{0}, w_{1})$ , les équations normales:

n w_{0} + (i \sum x_{i}) w_{1} (i \sum x_{i}) w_{0} + (i \sum x_{i}^{2}) w_{1} = i \sum y_{i}, = i \sum x_{i} y_{i} . (2.4)

La première, divisée par $n$ , s'écrit $w_{0} + w_{1} \overset{x}{ˉ} = \overset{y}{ˉ}$ , d'où : la droite passe par le point moyen. Substituons dans la seconde, après l'avoir réécrite :

i \sum x_{i} ((y_{i} - \overset{y}{ˉ}) - w_{1} (x_{i} - \overset{x}{ˉ})) = 0.

Comme $\sum_{i} (y_{i} - \overset{y}{ˉ}) = 0$ et $\sum_{i} (x_{i} - \overset{x}{ˉ}) = 0$ , on peut remplacer le facteur par sans changer la somme (on retranche fois une somme nulle). Il reste . Les n'étant pas tous égaux, , et . Le système (2.4) a donc une solution unique, (2.3).

C'est bien un minimum global. Annuler les dérivées ne suffit pas: il faut encore montrer que ce point critique est un minimum. Fixons $w_{1}$ . La fonction $w_{0} \mapsto \sum_{i} ((y_{i} - w_{1} x_{i}) - w_{0})^{2}$ est le risque quadratique d'une sur les nombres ; par le théorème 1.2, elle est minimale en leur moyenne , et seulement là. Pour tout , on a donc . Or

SCR (\overset{y}{ˉ} - w_{1} \overset{x}{ˉ}, w_{1}) = i \sum ((y_{i} - \overset{y}{ˉ}) - w_{1} (x_{i} - \overset{x}{ˉ}))^{2} = S_{y y} - 2 w_{1} S_{x y} + w_{1}^{2} S_{xx},

avec $S_{y y} = \sum_{i} (y_{i} - \overset{y}{ˉ})^{2}$ . C'est un trinôme en de coefficient dominant , minimal en son sommet et seulement là. Le minimum de est donc atteint en et en aucun autre point.

L'hypothèse «les $x_{i}$ ne sont pas tous égaux» n'est pas une précaution de style. Si tous les appartements avaient la même surface, il n'y aurait aucune information sur l'effet de la surface: toute pente conviendrait, pourvu que la droite passe par le point moyen. C'est le premier exemple de colinéarité, que la dernière section du chapitre étudie en général.

La figure 2.1 montre la droite et ses résidus. Trois logements (3, 5 et 7) sont nettement sous la droite, de 150 à 190 CHF; trois autres (2, 6 et 8) au-dessus, de 110 à 160 CHF. Ces écarts n'ont rien d'anormal: ils traduisent tout ce que la surface ne dit pas, l'étage, le quartier, l'état de l'immeuble. Nous verrons à la section sur la régression multiple qu'une deuxième caractéristique peut en expliquer une bonne partie.

Explorateur 2.1 · Choisir la droite à la main

Réglez l'ordonnée à l'origine $w_{0}$ et la pente $w_{1}$ . Chaque résidu est un segment vertical, et son carré est dessiné comme un vrai carré: l'erreur quadratique moyenne est l'aire moyenne de ces huit carrés. La droite en tirets est celle des moindres carrés. Les deux dernières valeurs affichées restent égales quoi que vous fassiez: c'est l'identité de Pythagore de la démonstration du théorème 2.2.

Ordonnée à l'origine w₀ (CHF)1 000

Pente w₁ (CHF/m²)12,0

MSE de votre droite (CHF²)

68 602

MSE des moindres carrés (CHF²)

19 309,17

Optimum (w₀; w₁)

524,2; 20,16

RMSE de votre droite (CHF)

261,9

Excès sur le minimum (CHF²)

49 292,8

Écart moyen des prédictions au carré (CHF²)

49 292,8

Dans l'explorateur, chaque carré a pour côté la longueur d'un résidu; l'erreur quadratique moyenne est l'aire moyenne de ces carrés, et chercher les moindres carrés revient littéralement à rendre ces carrés aussi petits que possible en moyenne. Essayez de rejoindre la droite en tirets à la main: vous verrez que l'on peut diminuer beaucoup de carrés à la fois en tournant la droite autour du point moyen, et que l'erreur est bien moins sensible à $w_{0}$ qu'à $w_{1}$ . Une variation de 1 CHF/m² sur la pente déplace la prédiction à 110 m² de 110 CHF, alors qu'une variation de 1 CHF sur l'ordonnée à l'origine ne la déplace que de 1 CHF. Ce déséquilibre d'échelle entre les deux paramètres n'a aucune importance pour la formule close (2.3), mais il en aura une grande pour la descente de gradient du chapitre 3.

Question 2.2

Complétez moindres_carres(xs, ys), qui renvoie le couple (w0, w1) de la droite des moindres carrés par les formules (2.3). Le squelette renvoie la droite nulle. Le programme affiche la droite du jeu A, la prédiction pour 80 m² et la somme des carrés des résidus.

moindres_carres.py

Lire les coefficients

Un modèle linéaire est apprécié parce qu'on peut le lire. Encore faut-il le lire correctement.

Les unités. La pente $\overset{w}{^}_{1} = 20, 1573$ s'exprime en CHF/m²: c'est l'augmentation du loyer prédit pour un m² supplémentaire. L'ordonnée à l'origine $\overset{w}{^}_{0} = 524, 245$ s'exprime en CHF: c'est le loyer prédit pour une surface nulle. Une phrase correcte sur le modèle est donc: «selon ce modèle, dans la gamme de 30 à 110 m² couverte par les données, chaque m² supplémentaire s'accompagne en moyenne de 20,16 CHF de loyer en plus».

L'ordonnée à l'origine n'est souvent qu'une extrapolation. Il n'existe pas d'appartement de 0 m², et le jeu A ne descend pas en dessous de 30 m². Les 524 CHF ne sont pas un «loyer de base» observé; c'est le point où la droite, prolongée bien au-delà des données, coupe l'axe. Le chapitre 1 parlait de «loyer de base» pour la règle de l'agence, ce qui est une façon commode de parler, pas une constatation. Si l'on veut un coefficient constant interprétable, il vaut mieux centrer la caractéristique: avec $x^{'} = x - 70$ , l'écart à 70 m², la pente ne change pas, et l'ordonnée à l'origine devient le loyer prédit pour 70 m², $1935, 26$ CHF, une valeur au milieu des données (exercice 2.2).

Changer d'unité change les coefficients, pas le modèle. Si l'on exprime les loyers en milliers de francs, tous les coefficients sont divisés par mille, $\overset{w}{^}_{1} = 0, 0201573$ kCHF/m², et les prédictions aussi: c'est le même modèle, écrit autrement. Si l'on standardise la surface, $z = (x - \overset{x}{ˉ}) / s_{x}$ avec l'écart type de population m², les moindres carrés en donnent et : un écart type de surface en plus, c'est 509,84 CHF de loyer en plus. Les prédictions sont, encore une fois, identiques. Le chapitre 3 montrera que cette standardisation, sans effet sur la solution exacte, change tout pour un algorithme itératif.

La pente et la corrélation. Divisons numérateur et dénominateur de (2.3) par $n$ : $\overset{w}{^}_{1}$ est le rapport de la covariance empirique de $x$ et $y$ à la variance empirique de $x$ (Probabilités et statistique, chapitre 7). En introduisant le coefficient de corrélation $r = S_{x y} / S_{xx} S_{y y}$ , on obtient

\overset{w}{^}_{1} = r \frac{s _{y}}{s _{x}}, (2.5)

où $s_{x}$ et $s_{y}$ sont les écarts types. Sur le jeu A, $r = 0, 9648$ . La pente est la corrélation, remise à l'échelle des unités: une corrélation proche de 1 dit que les points sont presque alignés, la pente dit de combien $y$ change avec $x$ . Les deux informations sont différentes; une pente faible peut aller avec une corrélation parfaite.

Un coefficient de régression n'est pas un effet causal

La pente de 20,16 CHF/m² décrit une association dans les données: les grands appartements de l'échantillon sont plus chers. Elle ne dit pas qu'agrandir un appartement de 10 m² augmenterait son loyer de 201,57 CHF. Les grands appartements diffèrent des petits par bien d'autres choses — le quartier, l'âge de l'immeuble, la vue — et la pente mélange l'effet de la surface avec celui de tout ce qui varie avec elle. La section sur la régression multiple montre un cas où ce mélange retourne carrément le signe d'un coefficient. Établir un effet causal demande une expérience contrôlée, ou des hypothèses explicites sur la façon dont les données ont été produites; un ajustement, aussi bon soit-il, ne le fait pas à lui seul.

Question 2.3

Une analyste réestime la droite du jeu A en mesurant la surface en dm² (1 m² = 100 dm²) au lieu de m², sans toucher aux loyers. Qu'obtient-elle?

La forme matricielle

Avec une seule caractéristique, les sommes de l'exemple 2.1 suffisent. Avec $d$ caractéristiques — la surface, la distance au centre, l'étage, l'âge de l'immeuble —, il faut une écriture qui ne dépende pas de $d$ . C'est le calcul matriciel.

Pour le jeu A, $d = 1$ et

\tilde{X} = 111 ⋮ 1 304550 ⋮ 110 \in R^{8 \times 2}, y = 120015501380 ⋮ 2900 .

Le produit $\tilde{X}^{⊤} \tilde{X}$ est une matrice $2 \times 2$ dont les coefficients sont les sommes de l'exemple 2.1: $n$ , $\sum_{i} x_{i}$ et . De même, a pour composantes et . Le système (2.4) n'était autre que . C'est le résultat général.

Démonstration. Conditions du premier ordre. Écrivons $SCR (w) = \sum_{i = 1}^{n} (y_{i} - \sum_{k = 0}^{d} \tilde{x}_{ik} w_{k})^{2}$ , où . La dérivée partielle par rapport à vaut

\frac{\partial SCR}{\partial w _{j}} = - 2 i = 1 \sum n \tilde{x}_{ij} (y_{i} - k = 0 \sum d \tilde{x}_{ik} w_{k}) = - 2 [\tilde{X}^{⊤} (y - \tilde{X} w)]_{j} .

Les $d + 1$ dérivées s'annulent simultanément si et seulement si $\tilde{X}^{⊤} (y - \tilde{X} w) = 0$ , c'est-à-dire si et seulement si $w$ vérifie (2.8). Tout minimiseur, qui est un point critique de la fonction dérivable $SCR$ , vérifie donc les équations normales.

Toute solution est un minimiseur global. Soit $\hat{w}$ une solution de (2.8) et $\hat{e} = y - \tilde{X} \hat{w}$ , de sorte que $\tilde{X}^{⊤} \hat{e} = 0$ . Pour tout ,

y - \tilde{X} w = \hat{e} + \tilde{X} (\hat{w} - w) .

Le produit scalaire des deux termes est nul: $(\tilde{X} (\hat{w} - w))^{⊤} \hat{e} = (\hat{w} - w)^{⊤} \tilde{X}^{⊤} \hat{e} = 0$ . Par le théorème de Pythagore, , ce qui est (2.10). Le second terme est positif ou nul: est un minimiseur global, et un autre fait aussi bien si et seulement si .

Existence et unicité. Montrons d'abord que $\tilde{X}^{⊤} \tilde{X}$ et $\tilde{X}$ ont le même noyau. Si $\tilde{X} v = 0$ , alors . Réciproquement, si , alors , donc . Par le théorème du rang (Algèbre linéaire, chapitre 5), les deux matrices, qui ont colonnes, ont donc le même rang, et ce rang est aussi celui de . Or l'image de est contenue dans celle de ; deux sous-espaces emboîtés de même dimension sont égaux. Le vecteur , qui est dans l'image de , est donc dans celle de : le système (2.8) a au moins une solution. Si les colonnes de sont indépendantes, le noyau commun est réduit à : la matrice carrée est inversible, la solution de (2.8) est unique et vaut (2.9).

La démonstration contient plus que l'énoncé. L'identité (2.10) dit de combien un $w$ quelconque est moins bon que l'optimum: exactement la norme au carré de l'écart entre ses prédictions et celles des moindres carrés. C'est la généralisation de l'identité (1.5) du chapitre 1, et c'est elle qu'affichent les deux dernières lectures de l'explorateur 2.1: l'excès d'erreur quadratique moyenne de votre droite sur l'optimum est égal, à l'arrondi près, à l'écart quadratique moyen entre vos prédictions et celles de la droite en tirets. Elle dit aussi que l'ensemble des minimiseurs est $\hat{w} + ker \tilde{X}$ : un seul point si les colonnes sont indépendantes, une droite ou un plan de solutions sinon, toutes avec les mêmes prédictions sur les données.

Démonstration. Les équations normales s'écrivent $\tilde{X}^{⊤} (y - \tilde{X} \hat{w}) = 0$ , soit $\tilde{X}^{⊤} \hat{e} = 0$ . La composante de ce vecteur est le produit scalaire de la colonne de avec . Pour , la colonne est le vecteur de uns, et la composante vaut ; pour , la colonne est celle de la caractéristique , et la composante vaut . Toutes sont nulles, ce qui est (2.11). Comme , la première égalité donne . Enfin, la prédiction au point moyen vaut .

Il faut bien voir le rôle de la colonne de uns dans ce résultat. C'est elle, et elle seule, qui force la somme des résidus à être nulle. Un modèle sans ordonnée à l'origine, $f (x) = w_{1} x$ , a des résidus orthogonaux à $x$ mais pas de somme nulle: ses erreurs sont, en moyenne, d'un côté. Le problème guidé 2.1 le calcule sur le jeu A.

La lecture géométrique: une projection

Les équations normales ont une interprétation géométrique qui vaut tous les calculs. Changeons de point de vue: au lieu de regarder les $n$ exemples comme $n$ points du plan, regardons les colonnes de $\tilde{X}$ et le vecteur $y$ comme des vecteurs de $R^{n}$ — un espace dont chaque coordonnée correspond à un exemple.

Quand $w$ parcourt $R^{d + 1}$ , le vecteur des prédictions $\tilde{X} w$ parcourt l'espace des colonnes $col (\tilde{X})$ , le sous-espace de $R^{n}$ engendré par les colonnes de . Pour le jeu A, c'est un plan de , engendré par le vecteur de uns et le vecteur des surfaces . Minimiser , c'est chercher . La réponse est connue (Algèbre linéaire, chapitre 7): c'est la de sur , caractérisée par le fait que est orthogonal à tout le sous-espace — donc à chacune de ses colonnes génératrices. C'est exactement . Les équations normales tirent d'ailleurs leur nom de là: elles disent que le résidu est , c'est-à-dire perpendiculaire, à l'espace des colonnes.

L'idempotence se vérifie en une ligne: $H^{2} = \tilde{X} (\tilde{X}^{⊤} \tilde{X})^{- 1} (\tilde{X}^{⊤} \tilde{X}) (\tilde{X}^{⊤} \tilde{X})^{- 1} \tilde{X}^{⊤} = H$ . Elle a une traduction concrète: si l'on remplace les loyers par les loyers ajustés et que l'on refait la régression, on retrouve exactement la même droite, puisque les nouveaux loyers sont déjà dans le plan.

Le cas $n = 3$ se dessine. Prenons trois observations, $x = (1, 2, 3)$ et $y = (1, 3, 2)$ . Les formules (2.3) donnent $\overset{x}{ˉ} = 2$ , $\overset{y}{ˉ} = 2$ , , , d'où et . Les valeurs ajustées sont , un vecteur du plan engendré par et , et le résidu vaut . On vérifie à la main que et .

Le théorème de Pythagore, dans ce triangle rectangle, s'écrit $∥ y ∥^{2} = ∥ \hat{y} ∥^{2} + ∥ \hat{e} ∥^{2}$ ; sur l'exemple, $14 = 12, 5 + 1, 5$ . La matrice de projection de cet exemple vaut

H = \frac{1}{6} 52 - 1 222 - 1 25,

et sa trace vaut $2$ , la dimension du plan sur lequel elle projette. C'est un fait général: la trace de $H$ est le nombre $d + 1$ de paramètres du modèle. Ce nombre de paramètres est une première mesure de la flexibilité d'un modèle, notion que le chapitre 6 développe.

La lecture géométrique éclaire trois questions d'un coup. Pourquoi la solution existe toujours: un sous-espace a toujours une projection orthogonale. Quand elle est unique: la projection $\hat{y}$ est toujours unique, mais les coefficients $\hat{w}$ qui l'expriment comme combinaison des colonnes ne le sont que si les colonnes forment une base, c'est-à-dire sont indépendantes. Ce que les moindres carrés ne font pas: ils ne «trouvent» aucune vérité sur $y$ ; ils l'approchent au mieux par ce que les colonnes permettent d'écrire, et tout ce qui est orthogonal à ces colonnes reste dans le résidu.

Le coefficient de détermination

L'erreur quadratique moyenne de 19 309,17 CHF² ne dit pas grand-chose seule. Le chapitre 1 a donné le bon réflexe: la comparer à celle du prédicteur constant, qui vaut $279243, 75$ CHF². C'est ce que fait le coefficient de détermination.

Comme $SCT / n$ est l'erreur quadratique moyenne du prédicteur constant $\overset{y}{ˉ}$ (théorème 1.2), $R^{2}$ mesure la fraction de l'erreur de ce prédicteur de référence que le modèle fait disparaître: $R^{2} = 1$ pour un modèle qui passe par tous les points, $R^{2} = 0$ pour un modèle qui ne fait pas mieux que la moyenne. Pour un modèle des moindres carrés avec ordonnée à l'origine, il a une seconde lecture, que voici.

Démonstration. Écrivons $y - \overset{y}{ˉ} 1 = (\hat{y} - \overset{y}{ˉ} 1) + \hat{e}$ . Le vecteur appartient à , puisque y appartient et que en est une colonne — c'est ici que l'ordonnée à l'origine sert. Par le théorème 2.3, est orthogonal à toutes les colonnes de , donc à toute combinaison de ces colonnes, en particulier à . Pythagore donne , qui est (2.14). Les trois termes étant positifs, et (on suppose , c'est-à-dire des non tous égaux).

En régression simple, $\overset{y}{^}_{i} - \overset{y}{ˉ} = \overset{w}{^}_{1} (x_{i} - \overset{x}{ˉ})$ par (2.3), donc , et .

Sur le jeu A, $SCT = 2233950$ CHF², et la décomposition s'écrit

2233950 = 2079476, 62 + 154473, 38, R^{2} = 1 - \frac{154 473 , 38}{2 233 950} = 0, 9309.

La surface «explique» 93 % de la variabilité des loyers de l'échantillon, au sens précis de (2.14). Et l'on vérifie $r^{2} = 0, 964 8^{2} = 0, 9309$ .

R^{2}

d'entraînement élevé ne promet rien sur des données nouvelles

Le $R^{2}$ calculé sur les données qui ont servi à ajuster le modèle est une erreur d'entraînement déguisée, et il hérite de tous ses défauts (chapitre 1). En particulier, ajouter une caractéristique ne peut jamais le faire baisser, même si cette caractéristique est du pur bruit: le modèle élargi peut toujours mettre son coefficient à zéro et retrouver l'ancien ajustement (exercice 2.4). Avec autant de paramètres que d'exemples, on obtient en général $R^{2} = 1$ sur l'entraînement, quelles que soient les données. Pour juger un modèle, on calcule son — ou son erreur quadratique — sur un ensemble de test (chapitre 5). Sur des données de test, la décomposition (2.14) ne tient plus, et peut être : c'est le cas d'un modèle qui prédit plus mal que la simple moyenne.

Question 2.4

La règle de l'agence $f (x) = 500 + 20 x$ a, sur le jeu A, une somme des carrés des résidus de $164400$ CHF² (exemple 1.1). Calculez son coefficient de détermination $R^{2} = 1 - SCR / SCT$ .

Régression multiple

Les résidus de la figure 2.1 ne sont pas du hasard pur: ils traduisent ce que la surface ne dit pas. Supposons que l'agence ait aussi relevé la distance de chaque appartement à la gare de Lausanne, en km. Ce relevé est propre à ce chapitre; ce sont, comme le reste, des données fictives, construites pour le cours.

Logement	1	2	3	4	5	6	7	8
Surface $x_{1}$ (m²)	30	45	50	62	70	85	95	110
Distance $x_{2}$ (km)	1,0	1,0	3,0	1,5	3,0	1,5	4,0	2,5
Loyer $y$ (CHF)	1 200	1 550	1 380	1 840	1 750	2 350	2 250	2 900

Les trois logements sous la droite de la figure 2.1 — les logements 3, 5 et 7 — sont aussi les plus éloignés du centre. Le modèle à deux caractéristiques s'écrit $f (x) = w_{0} + w_{1} x_{1} + w_{2} x_{2}$ , sa matrice de conception a trois colonnes (les uns, les surfaces, les distances), et les équations normales (2.8) forment un système .

C'est la lecture correcte d'un coefficient de régression multiple: $\overset{w}{^}_{j}$ est la variation de la prédiction quand $x_{j}$ augmente d'une unité, toutes les autres caractéristiques restant fixées (ceteris paribus). Elle explique aussi pourquoi la pente de la surface a changé, de $20, 16$ à $23, 59$ CHF/m²: la régression simple attribuait à la surface une partie de l'effet de la distance, qui jouait en sens contraire. L'exercice 2.5 établit la relation exacte entre les deux pentes; c'est le ().

Les caractéristiques catégorielles entrent dans le même cadre par les indicateurs du chapitre 1 (one-hot encoding), avec une précaution. Si un appartement est soit à Ouchy, soit au Flon, soit à la Sallaz, les trois indicateurs somment à 1 pour chaque ligne: leur somme est la colonne de uns. Avec l'ordonnée à l'origine, les colonnes de $\tilde{X}$ sont alors dépendantes. On omet donc une modalité, dite de référence, et le coefficient de chacune des autres se lit comme un écart à la référence (exercice 2.3).

Question 2.5

Complétez equations_normales(X, y), qui reçoit la matrice de données X (une liste de lignes, sans la colonne de uns) et les étiquettes y, et renvoie le couple (A, b) avec $A = \tilde{X}^{⊤} \tilde{X}$ et $b = \tilde{X}^{⊤} y$ , où $\tilde{X}$ est X précédée d'une colonne de uns. La fonction resoudre(A, b), une élimination de Gauss, est fournie. Le programme ajuste le modèle surface + distance de l'exemple 2.3. Attention: ne modifiez pas les lignes de X qu'on vous passe.

equations_normales.py

def resoudre(A, b):
  # Elimination de Gauss avec pivot partiel (fournie, ne pas modifier).
  n = len(A)
  M = [[float(v) for v in ligne] + [float(bi)] for ligne, bi in zip(A, b)]
  for c in range(n):
      p = max(range(c, n), key=lambda r: abs(M[r][c]))
      M[c], M[p] = M[p], M[c]
      for r in range(c + 1, n):
          f = M[r][c] / M[c][c]
          for k in range(c, n + 1):
              M[r][k] -= f * M[c][k]
  w = [0.0] * n
  for i in range(n - 1, -1, -1):
      s = sum(M[i][k] * w[k] for k in range(i + 1, n))
      w[i] = (M[i][n] - s) / M[i][i]
  return w

def equations_normales(X, y):
  # A completer: construire les lignes de X tilde ([1] + ligne de X),
  # puis A[j][k] = somme sur i de xt[i][j] * xt[i][k]
  # et b[j] = somme sur i de xt[i][j] * y[i].
  return [], []

# surface (m2), distance a la gare (km); loyer (CHF). Donnees fictives.
X = [[30, 1.0], [45, 1.0], [50, 3.0], [62, 1.5],
   [70, 3.0], [85, 1.5], [95, 4.0], [110, 2.5]]
y = [1200, 1550, 1380, 1840, 1750, 2350, 2250, 2900]
A, b = equations_normales(X, y)
w = resoudre(A, b)
print("w =", [round(v, 3) for v in w])
if len(w) == 3:
  print("prediction pour 80 m2 a 2 km: %.2f CHF" % (w[0] + w[1] * 80 + w[2] * 2))

Caractéristiques polynomiales

Le mot «linéaire» dans «régression linéaire» porte sur les paramètres, pas sur les caractéristiques. Rien n'empêche de calculer, à partir d'une caractéristique brute $t$ , de nouvelles caractéristiques $t^{2}$ , $t^{3}$ , $ln t$ ou $t$ , et de les mettre dans les colonnes de . Le modèle est une parabole en , mais une fonction de , et tout ce qui précède s'applique sans changer une ligne: équations normales, projection, .

Exemple 2.4 · Chauffer l'hiver, refroidir l'été

Un immeuble de bureaux est chauffé quand il fait froid et climatisé quand il fait chaud. Huit jours de relevés donnent sa consommation d'énergie en fonction de la température extérieure moyenne. Ce sont des données fictives, construites pour le cours, propres à ce chapitre.

Température $t$ (°C)	−5	0	5	10	15	20	25	30
Énergie $y$ (kWh/jour)	90	63	49	34	32	29	39	49

La droite des moindres carrés, avec la seule caractéristique $t$ , vaut $\overset{y}{^} = 62, 083 - 1, 1167 t$ . Elle a une RMSE de kWh/jour et un de : elle capte que l'on consomme plus par temps froid, et rien d'autre. Ses résidus le disent: ils valent à −5 °C, sont négatifs de 5 à 20 °C et redeviennent positifs, , à 30 °C. Un résidu qui change de signe de façon ordonnée le long de signale une forme que le modèle ne sait pas suivre.

La figure 2.3 montre l'idée en une image: le modèle linéaire n'est pas condamné aux droites. Avec $φ (t) = (t, t^{2})$ , il suit une parabole; avec des puissances plus élevées, des courbes plus tourmentées; avec des indicateurs d'intervalles, des fonctions en escalier; avec des produits $x_{1} x_{2}$ , des interactions — l'effet de la surface sur le loyer peut alors dépendre du quartier. Le chapitre 9 poussera l'idée à sa limite, avec des transformations à une infinité de composantes, calculées implicitement par un noyau.

Mais cette liberté a un prix, que la figure 2.3 laisse deviner. Huit points ont toujours un polynôme de degré 7 qui passe exactement par eux tous — c'est l'interpolation de Lagrange —, et son $R^{2}$ d'entraînement vaut 1. Ce polynôme oscille entre les points et part à l'infini hors de l'intervalle des données; il a appris le bruit avec le signal. Choisir le degré est un problème de surapprentissage, que l'on ne résout pas sur les données d'entraînement: c'est l'objet du chapitre 6. Même la parabole, honnête entre −5 et 30 °C, doit être prise avec prudence au-delà: elle prédit $68, 63$ kWh/jour à 35 °C et $118, 88$ à −10 °C, des extrapolations que rien dans les données ne vient confirmer.

Question 2.6

Complétez polynomiales(t, degre), qui renvoie la ligne $(1, t, t^{2}, \dots, t^{m})$ de la matrice de conception pour le degré $m$ , et r_carre(y, y_chapeau), qui renvoie $R^{2} = 1 - SCR / SCT$ . La fonction ajuster, qui résout les équations normales, est fournie. Le programme compare les degrés 1 et 2 sur les données de l'exemple 2.4.

polynomes.py

def resoudre(A, b):
  # Elimination de Gauss avec pivot partiel (fournie).
  n = len(A)
  M = [[float(v) for v in ligne] + [float(bi)] for ligne, bi in zip(A, b)]
  for c in range(n):
      p = max(range(c, n), key=lambda r: abs(M[r][c]))
      M[c], M[p] = M[p], M[c]
      for r in range(c + 1, n):
          f = M[r][c] / M[c][c]
          for k in range(c, n + 1):
              M[r][k] -= f * M[c][k]
  w = [0.0] * n
  for i in range(n - 1, -1, -1):
      s = sum(M[i][k] * w[k] for k in range(i + 1, n))
      w[i] = (M[i][n] - s) / M[i][i]
  return w

def ajuster(lignes, y):
  # Moindres carres sur des lignes de X tilde deja construites (fournie).
  p = len(lignes[0])
  A = [[sum(r[j] * r[k] for r in lignes) for k in range(p)] for j in range(p)]
  b = [sum(r[j] * yi for r, yi in zip(lignes, y)) for j in range(p)]
  return resoudre(A, b)

def polynomiales(t, degre):
  # A completer: renvoyer [1, t, t**2, ..., t**degre] (degre + 1 valeurs).
  return [1]

def r_carre(y, y_chapeau):
  # A completer: 1 - SCR / SCT, avec SCT calculee autour de la moyenne de y.
  return 0.0

temperatures = [-5, 0, 5, 10, 15, 20, 25, 30]
energie = [90, 63, 49, 34, 32, 29, 39, 49]
for degre in [1, 2]:
  lignes = [polynomiales(t, degre) for t in temperatures]
  w = ajuster(lignes, energie)
  y_chapeau = [sum(wj * fj for wj, fj in zip(w, f)) for f in lignes]
  print("degre %d: R2 = %.4f" % (degre, r_carre(energie, y_chapeau)))

La lecture probabiliste

Jusqu'ici, les moindres carrés étaient un choix: celui de la perte carrée. On peut aussi les déduire d'une hypothèse sur la façon dont les données sont produites. C'est l'objet de cette section, et c'est l'argument qui fait de la régression linéaire un modèle statistique et pas seulement un ajustement de courbe.

Le modèle dit que les loyers sont dispersés autour d'une droite, symétriquement, avec la même amplitude typique $σ$ pour tous les appartements, et que les grandes déviations sont rares — la loi normale donne une probabilité d'environ 0,3 % à un écart de plus de trois écarts types. Si l'on accepte ce modèle, le principe du maximum de vraisemblance (Probabilités et statistique, chapitre 10) dit comment estimer $w$ : choisir la valeur sous laquelle les données observées sont les plus probables.

Démonstration. Sous (2.16), $y_{i}$ a pour densité $\frac{1}{2 π σ ^{2}} exp (- \frac{( y _{i} - x ~ _{i}^{⊤} w ) ^{2}}{2 σ ^{2}})$ . Les étant indépendants, la densité conjointe des étiquettes — la vraisemblance — est le produit de ces densités:

L (w, σ^{2}) = i = 1 \prod n \frac{1}{2 π σ ^{2}} exp (- \frac{( y _{i} - x ~ _{i}^{⊤} w ) ^{2}}{2 σ ^{2}}) .

Le logarithme transforme le produit en somme: $ln L = - \frac{n}{2} ln (2 π σ^{2}) - \frac{1}{2 σ ^{2}} \sum_{i} (y_{i} - \tilde{x}_{i}^{⊤} w)^{2}$ , qui est (2.17). Pour fixé, le premier terme ne dépend pas de , et le second est , avec . Maximiser en revient donc à minimiser : les maximiseurs sont exactement les solutions des moindres carrés, quel que soit . Le logarithme étant croissant, maximiser ou est la même chose.

Remplaçons $w$ par $\hat{w}$ et notons $s = SCR (\hat{w}) > 0$ . La fonction $g (v) = - \frac{n}{2} ln (2 π v) - \frac{s}{2 v}$ , pour , a pour dérivée

g^{'} (v) = - \frac{n}{2 v} + \frac{s}{2 v ^{2}} = \frac{s - n v}{2 v ^{2}},

positive pour $v < s / n$ et négative pour $v > s / n$ . Elle est donc maximale en $v = s / n$ , ce qui est (2.18). Comme pour tout $(w, v)$ on a $ln L (w, v) \leq ln L (\hat{w}, v) \leq ln L (\hat{w}, s / n)$ , le couple maximise la vraisemblance.

Sur le jeu A, l'écart type du bruit estimé par maximum de vraisemblance est la RMSE, $\overset{σ}{^} = 138, 96$ CHF. On démontre en statistique que cet estimateur de $σ^{2}$ est biaisé vers le bas — les résidus sont plus petits que le bruit, puisque la droite a été choisie pour les rendre petits — et que $SCR / (n - d - 1)$ ne l'est pas; nous l'admettons ici. Avec $n - d - 1 = 6$ , cela donne CHF², soit un écart type de CHF. Le diviseur compte les des résidus: résidus, liés par les contraintes d'orthogonalité (2.11).

Le théorème 2.5 a une portée plus large que la régression linéaire. Il dit qu'une perte est une hypothèse sur le bruit déguisée. Si l'on suppose un bruit de Laplace, de densité $\frac{1}{2 b} e^{- ∣ ε ∣ / b}$ , le même calcul donne une log-vraisemblance $- n ln (2 b) - \frac{1}{b} \sum_{i} ∣ y_{i} - \tilde{x}_{i}^{⊤} w ∣$ , et le maximum de vraisemblance minimise la somme des des résidus: la perte absolue du chapitre 1, dont le minimiseur constant était la médiane. La loi de Laplace a des queues plus lourdes que la loi normale; l'hypothèse correspondante accepte mieux les valeurs extrêmes, et l'estimateur qui en résulte y est moins sensible. Le chapitre 4 refera ce raisonnement pour une étiquette binaire, et en tirera l'entropie croisée.

Ce que l'hypothèse gaussienne n'autorise pas

Les moindres carrés calculent toujours quelque chose, que le bruit soit gaussien ou non; l'hypothèse ne sert qu'à leur donner une justification probabiliste, et à en tirer des intervalles de confiance et des tests (Probabilités et statistique, chapitres 11 et 12). Elle est souvent fausse en détail. Sur des loyers, l'écart à la droite grandit probablement avec la surface — se tromper de 150 CHF sur un studio n'est pas la même chose que sur un cinq-pièces —, ce qui contredit la variance constante; on parle d'hétéroscédasticité. Une erreur de saisie produit un écart que la loi normale juge presque impossible, et les moindres carrés lui donnent un poids énorme (exemple 1.2). Avant de croire un intervalle de confiance tiré de (2.16), regardez les résidus.

Question 2.7

Remettez dans l'ordre les étapes de la démonstration du théorème 2.5 (moindres carrés et maximum de vraisemblance).

Glissez les éléments pour les mettre dans le bon ordre

Constater que seul le terme

- SCR (w) / (2 σ^{2})

dépend de

w

Prendre le logarithme pour transformer le produit en somme

Annuler la dérivée en

σ^{2}

pour obtenir

\overset{σ}{^}^{2} = SCR (\hat{w}) / n

Multiplier les

n

densités, ce que permet l'indépendance des bruits

Conclure que maximiser en

w

revient à minimiser la somme des carrés des résidus

Écrire la densité normale de chaque

y_{i}

, d'espérance

\tilde{x}_{i}^{⊤} w

et de variance

σ^{2}

Quand $\tilde{X}^{⊤} \tilde{X}$ n'est pas inversible

Le théorème 2.2 suppose les colonnes de $\tilde{X}$ indépendantes. Que se passe-t-il quand elles ne le sont pas?

Trois situations la produisent couramment. Une caractéristique redondante: la surface en m² et la même surface en pieds carrés, ou un prix TTC à côté du prix HT avec un taux unique. Des indicateurs complets avec ordonnée à l'origine: les indicateurs de toutes les modalités d'une variable catégorielle somment à la colonne de uns, comme on l'a vu. Plus de paramètres que d'exemples: si $d + 1 > n$ , les $d + 1$ colonnes vivent dans $R^{n}$ et ne peuvent pas être indépendantes. Ce dernier cas n'est pas exotique; il est la règle en génomique, où l'on mesure des milliers de gènes sur quelques centaines de patients, et en traitement du texte.

Colinéarité exacte. Supposons que chaque appartement du jeu A ait un balcon de 8 m², et que l'on ajoute au modèle la surface totale $x_{3} = x + 8$ . La troisième colonne de $\tilde{X}$ est la deuxième plus huit fois la première: avec $v = (8, 1, - 1)$ , on a $\tilde{X} v = 8 \cdot 1 + x - (x + 8 \cdot 1) = 0$ . La matrice

\tilde{X}^{⊤} \tilde{X} = 854761154742519468956114689551783

a un déterminant nul. Le théorème 2.2 dit ce qui se passe alors: les équations normales ont une infinité de solutions, $\hat{w} + t v$ pour tout réel $t$ , à partir par exemple de $\hat{w} = (524, 245; 20, 1573; 0)$ , et toutes donnent exactement les mêmes prédictions. Le modèle à trois coefficients ne sait pas répartir l'effet de la surface entre «surface habitable» et «surface totale», puisque ces deux colonnes portent la même information. Les prédictions restent bien définies; ce sont les coefficients qui ne le sont plus.

Quasi-colinéarité. Le cas exact est rare dans des données réelles; le cas presque exact est fréquent, et plus traître. Supposons maintenant que le logement 8 ait un balcon de 9 m², tous les autres gardant 8 m². Les colonnes ne sont plus dépendantes, $\tilde{X}^{⊤} \tilde{X}$ est inversible, et la solution est unique:

Balcon du logement 8	$\overset{w}{^}_{0}$ (CHF)	$\overset{w}{^}_{1}$ (surface habitable)	$\overset{w}{^}_{3}$ (surface totale)

Un mètre carré de balcon en plus ou en moins, sur un seul appartement, suffit à faire basculer tous les coefficients, qui prennent des valeurs absurdes: avec un balcon de 9 m², chaque m² habitable retirerait 277,62 CHF de loyer, tandis que chaque m² total en ajouterait 295,37. La matrice $\tilde{X}^{⊤} \tilde{X}$ a des valeurs propres de $94521$ et de $0, 00896$ : un rapport d'environ $1, 05 \cdot 1 0^{7}$ . Dans la direction propre de la plus petite valeur propre, la somme des carrés est presque plate, et un minuscule changement des données déplace le minimum très loin le long de cette direction. Les prédictions sur les huit appartements restent raisonnables — le est le même dans les deux lignes du tableau —, mais pour un nouvel appartement de 80 m² avec un balcon de 9 m², les deux modèles prédisent et CHF.

Même sans colinéarité, le conditionnement de $\tilde{X}^{⊤} \tilde{X}$ compte. Pour le jeu A seul, ses valeurs propres valent environ $42526$ et $0, 963$ , dans un rapport d'environ $44170$ , simplement parce que la colonne de uns et la colonne des surfaces n'ont pas la même échelle. La formule close s'en accommode; le chapitre 3 montre que la descente de gradient, elle, en souffre beaucoup, et comment la standardisation y remédie.

Les moindres carrés en pratique

Pour information, voici les équivalents usuels, que vous n'avez pas à exécuter ici. numpy.linalg.lstsq résout le problème des moindres carrés par une décomposition en valeurs singulières, sans former $\tilde{X}^{⊤} \tilde{X}$ , et renvoie la solution de plus petite norme en cas de colinéarité; scikit-learn ajoute lui-même l'ordonnée à l'origine.

import numpy as np
from sklearn.linear_model import LinearRegression

Synthèse

La régression linéaire minimise la somme des carrés des résidus sur le modèle $f (x) = \tilde{x}^{⊤} w = w_{0} + w_{1} x_{1} + \dots + w_{d} x_{d}$ . En dimension un, et ; sur le jeu A, , avec une RMSE de CHF et une prédiction de CHF 2 136,83 pour 80 m².

Série d'exercices du chapitre 2Exercice 1 sur 5

Question 2.8

Avec la droite des moindres carrés du jeu A, $\overset{y}{^} = 524, 245 + 20, 1573 x$ , quel loyer prédit-on pour un appartement de 100 m², en CHF?

Exercices

Vous pouvez afficher le corrigé directement sous chaque énoncé après avoir cherché la solution.

Exercice 2.1 · Une régression à la main

On considère quatre observations: $(1, 3)$ , $(2, 4)$ , $(4, 8)$ et $(5, 9)$ .

Calculez $\overset{x}{ˉ}$ , $\overset{y}{ˉ}$ , et , puis la droite des moindres carrés.

Exercice 2.2 · Changer d'unités et centrer

On reprend la droite des moindres carrés du jeu A, $\overset{w}{^}_{0} = 524, 245$ CHF et $\overset{w}{^}_{1} = 20, 1573$ CHF/m².

Exercice 2.3 · Une variable catégorielle

Six appartements (données fictives) sont situés dans trois quartiers: deux au Flon, loués 1 800 et 2 000 CHF; trois à Ouchy, loués 2 600, 2 400 et 2 500 CHF; un à la Sallaz, loué 1 500 CHF. On code le quartier par des indicateurs $u_{F}$ , $u_{O}$ , $u_{S}$ .

Exercice 2.4 · Ajouter

x^{2}

au jeu A

On ajuste sur le jeu A le modèle $f (x) = w_{0} + w_{1} x + w_{2} x^{2}$ . Les moindres carrés donnent , et , avec une somme des carrés des résidus de CHF².

Exercice 2.5 · Le biais de variable omise (démonstration)

On considère deux modèles sur les mêmes données: la régression simple de $y$ sur $x$ , de coefficients $(\tilde{w}_{0}, \tilde{w}_{1})$ , et la régression de sur et une seconde caractéristique , de coefficients . On note enfin les coefficients de la régression simple de sur . Toutes les matrices de conception sont supposées de rang plein.

Références

James, G., Witten, D., Hastie, T. et Tibshirani, R., An Introduction to Statistical Learning, Springer, 2ᵉ éd., 2021, chap. 3 (régression linéaire simple et multiple, interprétation des coefficients, colinéarité, extensions non linéaires).
Hastie, T., Tibshirani, R. et Friedman, J., The Elements of Statistical Learning, Springer, 2ᵉ éd., 2009, chap. 3 (moindres carrés, géométrie de la projection, régression multiple par orthogonalisations successives).
Bishop, C. M., Pattern Recognition and Machine Learning, Springer, 2006, chap. 3 (modèles linéaires en leurs paramètres, fonctions de base, maximum de vraisemblance sous bruit gaussien).
Murphy, K. P., Probabilistic Machine Learning: An Introduction, MIT Press, 2022, chap. 11 (régression linéaire, lecture probabiliste, solution par décomposition QR et SVD).
Azencott, C.-A., Introduction au Machine Learning, Dunod, 2ᵉ éd., 2022 (le vocabulaire et la régression linéaire, en français).

Logement	1	2	3	4	5	6	7	8
$\overset{y}{^}_{i}$ (CHF)	1 128,96	1 431,32	1 532,11	1 774,00	1 935,26	2 237,61	2 439,19	2 741,55
$e_{i}$ (CHF)	71,04	118,68	−152,11	66,00	−185,26	112,39	−189,19	158,45

9 m²	−1 711,39	−277,62	295,37	0,9518
7 m²	3 014,55	313,13	−295,37	0,9518

Objectifs du chapitre

Le problème: la meilleure droite

Les moindres carrés en dimension un

Lire les coefficients

La forme matricielle

La lecture géométrique: une projection

Le coefficient de détermination

Régression multiple

Caractéristiques polynomiales

La lecture probabiliste

Quand $\tilde{X}^{⊤} \tilde{X}$ n'est pas inversible

Synthèse

Exercices

Références

La pente

Une prédiction

La somme des résidus

Le prix de la contrainte

Objectifs du chapitre

Le problème: la meilleure droite

Les moindres carrés en dimension un

Lire les coefficients

La forme matricielle

La lecture géométrique: une projection

Le coefficient de détermination

Régression multiple

Caractéristiques polynomiales

La lecture probabiliste

Quand X~⊤X~ n'est pas inversible

Synthèse

Exercices

Références

La pente

Une prédiction

La somme des résidus

Le prix de la contrainte

Quand $\tilde{X}^{⊤} \tilde{X}$ n'est pas inversible