Régression linéaire

In [1]:
%matplotlib widget
import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact

Modèle

Une regression linéaire consiste à modéliser une variable à expliquer $Y$ par une somme linéaire de variables explicatives $X_0, ..., X_p$.
Dans cet exemple, nous utiliserons $Y$ telle que $\mathbb{E}[Y] = \sin$ et des variables explicatives polynomiales $X_k = x \mapsto x^k$, ce qui revient à chercher à approcher $Y$ par: $$\hat{Y}(x) = X(x) {\alpha} = \alpha_0 + \alpha_1 x^1 + \alpha_2 x^2 + ... + \alpha_p x^p$$ où $X(x) = \begin{bmatrix} 1 & x^{1} & x^2 & \cdots & x^{p} \end{bmatrix}$ et ${\alpha} = \begin{bmatrix} \alpha_{0} \\ \alpha_{1} \\ \vdots \\ \alpha_{p} \end{bmatrix}$

$\hat{Y}(x)$ va être la prédiction de notre modèle pour une valeur $x$.

Méthode des moindres carrés

Etant donné des échantillons $(x_i, y_i)$, notons $X$ la matrice des variables explicatives : $$X = \begin{bmatrix} X(x_0) \\ X(x_1) \\ \vdots \\ X(x_n) \end{bmatrix} = \begin{bmatrix} 1 & x_0^{1} & x_0^2 & \cdots & x_0^{p} \\ \vdots & \vdots & \vdots & \cdots & \vdots \\ 1 & x_n^{1} & x_n^2 & \cdots & x_n^{p} \end{bmatrix}$$ Le vecteur des prédictions devient : $$\hat{Y} = \begin{bmatrix} \hat{Y}(x_0) \\ \hat{Y}(x_1) \\ \vdots \\ \hat{Y}(x_n) \end{bmatrix} = X \alpha$$ Par abus de langage, on nomme encore $Y$ le vecteur des valeurs de la variable aléatoire $Y$ : $$Y = \begin{bmatrix} {Y}(x_0) \\ {Y}(x_1) \\ \vdots \\ {Y}(x_n) \end{bmatrix}$$ Pour savoir si notre prédiction est bonne, il nous faut un critère : on se donne une distance $d$ et on souhaite minimiser la distance $d(Y, \hat{Y})$ de $Y$ à $\hat{Y}$.
Dans la cas de la méthode des moindres carrés, on cherche à minimiser la distance euclidienne.
Remarque : dans l'idéal, $\hat{Y}$ est le projeté orthogonal de $Y$ sur l'espace vectoriel engendré par les $X_k$.

In [2]:
n = 10  # nombre d'échantillons
x = np.linspace(0, 4, n)  # vecteur des x_i (ici choisis uniformément sur [0, 4])
Y = np.random.normal(np.sin(x), 0.5, size = (1, n))  # Y = sin + gaussienne

Équation normale

On cherche donc à minimiser : $$|| Y - \hat Y ||^2 = || Y - X \alpha ||^2 = {}^t(Y - X \alpha)(Y - X \alpha) = {}^tY Y - 2{}^t \alpha {}^tX Y + {}^t \alpha {}^tX X \alpha$$ Cette fonction vectorielle en $\alpha$ est convexe et admet un minimum globale en un point qui annule sa dérivée : $$\boxed{ {}^tX Y - {}^tX X \alpha = 0}$$ Cette équation est appelée équation normale.
Sous réserve que ${}^tX X$ soit inversible, on obtient donc : $$\boxed{\alpha = ({}^tX X)^{-1}~{}^tX Y}$$ Un raffinement de la régression linéaire (ridge regression) consiste à ajouter un terme de régularisation $\lambda$ dans l'objectif de diminuer les coefficients de \alpha et d'éviter l'overfitting : $$\boxed{\alpha = ({}^tX X + \lambda I)^{-1}~{}^tX Y}$$

Code Python

In [3]:
p = n // 2  # degré du polynôme utilisé pour la régression

def normale_eq(Y, X, p):  # renvoie le α optimal
    return np.linalg.inv(X @ X.T + 10*np.eye(p)) @ X @ Y.T

def Y_hat(α, x):  # prédiction
    return np.sum(α[i]*x**i for i in range(len(α)))

X = [np.array([x**k for k in range(p)]) for p in range(1, n+1)]  # matrices X pour différentes valeurs de p
α = [normale_eq(Y, X[p], p+1) for p in range(n)]  # valeurs de α correspondantes

X_plot = np.linspace(0, 4, 100)
sinX = np.sin(X_plot)
Y_plot = [Y_hat(α[p], X_plot) for p in range(n)]
erreur_test = [np.sum((Y_plot[p] - sinX)**2) for p in range(n)]
erreur_train = [np.sum((Y_hat(α[p], x) - Y)**2) for p in range(n)]

fig = plt.figure(figsize=(10,4))
ax1 = plt.subplot(121)
ax1.scatter(x, Y, label="Données")
ax1.plot(X_plot, sinX, linestyle='dotted', color="red", label="$Y$ (sin)")
plot_pred, = ax1.plot(X_plot, Y_plot[p], label="$\hat Y$")
ax1.legend()

ax2 = fig.add_subplot(122)
ax2.plot(np.arange(n), erreur_test, label = "Erreur entre $Y$ et $\hat Y$")
plot_err_test = ax2.scatter([p], [erreur_test[p]], s = 100)
ax2.plot(np.arange(n), erreur_train, label = "Erreur entre données et $\hat Y$")
ax2.set_xlabel("p")
ax2.legend()
plot_err_train = ax2.scatter([p], [erreur_train[p]], s = 100)
plt.tight_layout()

@interact(p = (0, n-1))
def regression(p = n//2):  
    plot_pred.set_ydata(Y_plot[p])
    plot_err_test.set_offsets([[p, erreur_test[p]]])
    plot_err_train.set_offsets([p, erreur_train[p]])
    plt.draw()

Résultats

Sans terme de régularisation :

Avec régularisation ($\lambda = 1$), pour les mêmes données :

Avec scikit-learn

In [4]:
from sklearn.linear_model import LinearRegression
reg = LinearRegression().fit(X[p].T, Y.T)
In [5]:
X0 = (X_plot.reshape(-1, 1)**np.arange(p+1))
In [6]:
%matplotlib inline 
plt.plot(X_plot, reg.predict(X0)[:, 0])
plt.plot(X_plot, sinX)
plt.show()