%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$.
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
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}$$
x
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)]
erreur_donnees = np.sum((Y - np.sin(x))**2)
fig = plt.figure(figsize=(10,4))
ax1 = plt.subplot(131)
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(132)
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()
ax3 = fig.add_subplot(133)
ax3.plot(np.arange(n), erreur_test/erreur_donnees, label = "$R^2$")
ax3.set_xlabel("p")
ax3.legend()
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()
Pour voir le notebook en version interactive : https://mybinder.org/v2/gh/fortierq/ML/HEAD?filepath=regression_lineaire.ipynb
Sans terme de régularisation :
Avec régularisation ($\lambda = 1$), pour les mêmes données :
from sklearn.linear_model import LinearRegression
reg = LinearRegression().fit(X[p].T, Y.T)
X0 = (X_plot.reshape(-1, 1)**np.arange(p+1))
%matplotlib inline
plt.plot(X_plot, reg.predict(X0)[:, 0])
plt.plot(X_plot, sinX)
plt.show()