Régression par moindres carrés - Exemples

L’idée

À partir d’observations $(x_i,y_i)$, la régression linéaire simple cherche une droite $y=mx+b$. Un résidu est l’écart vertical $y_i-(mx_i+b)$ entre une observation et sa prédiction.

La méthode des moindres carrés choisit $m$ et $b$ en minimisant la somme des carrés des résidus :

$E(m,b)=\sum_{i=1}^{n}\left(y_i-(mx_i+b)\right)^2.$

Le carré empêche les résidus positifs et négatifs de s’annuler et pénalise davantage les grandes erreurs. Lorsque les entrées ne sont pas constantes, cet objectif est convexe et possède un minimum unique.

Une variable explicative

Pour une régression simple, dériver $E(m,b)$ par rapport à $m$ et $b$ donne

$\hat m=\frac{\sum_{i=1}^n(x_i-\bar x)(y_i-\bar y)}{\sum_{i=1}^n(x_i-\bar x)^2},\qquad \hat b=\bar y-\hat m\bar x.$

Le dénominateur est non nul lorsque les valeurs observées $x_i$ ne sont pas toutes identiques.

Plusieurs variables explicatives

Pour $n$ observations et $p$ variables explicatives, rassemblons les observations dans une matrice de conception $X\in\mathbb R^{n\times(p+1)}$. Sa première colonne contient des 1 pour l’ordonnée à l’origine ; chaque autre colonne contient une variable. Les coefficients inconnus forment $\beta\in\mathbb R^{p+1}$ et les valeurs observées forment $y\in\mathbb R^n$ :

$y=X\beta+\varepsilon.$

L’estimateur des moindres carrés est

$\hat\beta=\operatorname*{arg\,min}_{\beta}\|X\beta-y\|_2^2.$

Son gradient vaut

$\nabla_{\beta}\|X\beta-y\|_2^2=2X^T(X\beta-y).$

Au minimum, ce gradient est nul, ce qui donne les équations normales

$X^TX\hat\beta=X^Ty.$

Si $X$ est de rang colonne plein, $X^TX$ est inversible et

$\hat\beta=(X^TX)^{-1}X^Ty.$

Sans cette hypothèse, le minimiseur peut ne pas être unique. La solution de norme minimale est $\hat\beta=X^+y$, où $X^+$ est la pseudo-inverse de Moore–Penrose. En calcul numérique, on préfère généralement une décomposition QR ou en valeurs singulières au calcul explicite de $(X^TX)^{-1}$.

Exemple en Python

Cet exemple génère des observations bruitées autour de $y=2x+3$ et résout la régression avec numpy.linalg.lstsq :

import numpy as np

rng = np.random.default_rng(42)
x = np.linspace(0, 10, 100)
y = 2 * x + 3 + rng.normal(0, 1, len(x))

X = np.column_stack((x, np.ones(len(x))))
coefficients, residuals, rank, singular_values = np.linalg.lstsq(
    X, y, rcond=None
)

pente, ordonnee_origine = coefficients
predictions = X @ coefficients
print(pente, ordonnee_origine)

La construction reste identique avec davantage de variables : il suffit d’ajouter une colonne à $X$ pour chaque variable. L’objet ajusté devient alors un hyperplan dans l’espace des variables ; il ne s’obtient pas en intersectant plusieurs hyperplans séparés.