Le but de ce premier exercice est de se familiariser avec l'environnement Jupyter Notebook et l'utilisation de quelques fonctionalités de la bibliothèque Python scikit-learn. Nous allons nous intéresser à un problème de régression classique, sous l'angle des problématiques de l'apprentissage automatique. Le carnet est adapté de ce tutoriel.
Exécutez les cellules de ce carnet Jupyter les unes après les autres. En cas de problème d'exécution du code Python, vous pouvez redémarrer le noyau / kernel (onglet dans la barre du carnet Jupyter en haut). Pour ajouter vos commentaires personnels, créez une nouvelle cellule (onglet Insert) de type Markdown (à sélectionner dans le menu déroulant en haut du carnet). Un aide-mémoire pour Markdown est disponible ici. Vous pouvez aussi examiner le code source des cellules de cette page à l'aide du "double clic".
# import des bibliothèques Python utiles:
import numpy as np
import sklearn.linear_model as lm
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt
# "magic function" Jupyter pour l'affichage des graphiques dans le carnet:
%matplotlib inline
Dans la cellule ci-dessous, on génère des données y selon un modèle dépendant de $x\in[0,1.5]$. Quel est ce modèle, et quel type de bruit affecte les données?
# modèle génératif:
# les x_mod et y_mod serviront uniquement pour les représentations graphiques
x_mod = np.linspace(0., 1.5, 150)
print("taille de x_mod avant reshape: "+str(x_mod.shape))
x_mod = x_mod.reshape(len(x_mod),1) # scikit-learn exige une ligne par observation
print("taille de x_mod après reshape: "+str(x_mod.shape)) # que fait reshape?
y_mod = np.exp(3*x_mod)
# données:
# les y correspondant à certaines valeurs de x sont générés selon le modèle, et on ajoute un bruit
x_data = np.array([0, .1, .2, .5, .8, .9, 1])
x_data = x_data.reshape(len(x_data),1)
y_data = np.exp(3*x_data) + 2.0 * np.random.randn(len(x_data),1)
print("données x_data:")
print(x_data)
print("données y_data:")
print(y_data)
La cellule suivante permet une représentation graphique des données (on parle aussi des observations) et du modèle génératif.
# représentation graphique:
plt.figure(figsize=(10,6))
plt.plot(x_mod, y_mod, '--k')
plt.plot(x_data, y_data, 'or', ms=10)
plt.xlim(0, 1.5)
plt.xlabel("x")
plt.ylim(-10, 80)
plt.ylabel("y")
plt.title('observations et modèle génératif')
plt.legend(["modèle génératif","observations"]); # rem: le ";" évite des messages intempestifs sous Jupyter
Nous allons à présent essayer de créer des modèles permettant de prédire des valeurs $y$ correspondant à des valeurs de $x$ entre 0 et 1.5. Bien entendu, ces modèles seront créés à partir des données $(x_{data},y_{data})$. Nous confronterons nos prédictions aux valeurs obtenues par le modèle génératif des données, mais, bien entendu, dans une application réelle, le modèle génératif serait inconnu.
Remarque: attention, le mot "modèle" est utilisé dans deux sens différents...
Nous commençons par un modèle de régression linéaire. Cela consiste à prédire les valeurs de $y$ par une fonction affine de $x$: $$ y_{pred} = \alpha_0 + \alpha_1 x$$
Les valeurs de $\alpha_0$ et $\alpha_1$ sont estimées par la méthode des moindres carrés: on cherche les paramètres $\alpha_0$ et $\alpha_1$ qui minimisent $$\sum_{i=1}^n |y_{data}[i] - \alpha_0 - \alpha_1 x_{data}[i]|^2$$ sur l'ensemble des $n$ observations $(x_{data}[i],y_{data}[i])_{1\leq i\leq N}$ (ici, $n=7$).
# On crée un objet scikit-learn pour la régression linéaire:
lr = lm.LinearRegression()
# lorsqu'on crée un objet scikit-learn, on dispose de méthodes et attributs
# cf la documentation: on ne se servira que de quelques uns d'entre eux
# On estime les paramètres alpha_0 et alpha_1:
lr.fit(x_data, y_data)
# de manière générale, la méthode fit permet l'apprentissage des paramètres du modèle
# (ici par la méthodes des moindres carrés), qui sont stockés dans les attributs suivants:
print(lr.intercept_)
print(lr.coef_)
# On prédit des valeurs de y pour les x entre 0 et 1.5
y_pred_lr = lr.predict(x_mod)
# la méthode predict permet de prédire les valeurs y pour les x passés en argument
Quelles sont les valeurs de $\alpha_0$ et $\alpha_1$ ? Comparez aux valeurs obtenues par votre voisin.
La cellule suivante permet la représentation graphique de la prédiction par régression linéaire.
# représentation graphique:
plt.figure(figsize=(10,6))
plt.plot(x_mod, y_mod, '--k')
plt.plot(x_data, y_data, 'or', ms=10)
plt.plot(x_mod, y_pred_lr, '-g')
plt.xlim(0, 1.5)
plt.xlabel("x")
plt.ylim(-10, 80)
plt.ylabel("y")
plt.title('observations, modèle, et régression linéaire')
plt.legend(["modèle","observations","droite de régression"]);
Comme on le voit, la droite de régression passe globalement entre les points rouges, mais est en fait assez peu représentative du modèle génératif. Cela était prévisible, car ce modèle n'est pas linéaire.
La cellule suivante permet d'afficher l'erreur quadratique moyenne de prédiction sur les observations. Comment est défini cet indicateur?
print("régression linéaire, MSE = "+str(mean_squared_error(y_data,lr.predict(x_data))))
Nous allons à présent construire un modèle prédictif polynomial (donc non linéaire): $$ y_{pred} = \alpha_0 + \sum_{i=1}^d \alpha_d x^d$$ pour un entier $d>1$.
Les coefficients $\alpha_0, \dots, \alpha_d$ sont toujours estimés par la méthode des moindres carrés.
Pour ce faire, on utilise une régression linéaire multivariée: au lieu de la seule variable $x$, on introduit les variables $x, x^2, x^3,\dots, x^d$.
# création des vecteurs "puissances"
x_data2=np.power(x_data,2)
x_data3=np.power(x_data,3)
x_data4=np.power(x_data,4)
x_data5=np.power(x_data,5)
x_mod2=np.power(x_mod,2)
x_mod3=np.power(x_mod,3)
x_mod4=np.power(x_mod,4)
x_mod5=np.power(x_mod,5)
# modèle polynomial avec $d=2$
lrp2 = lm.LinearRegression()
lrp2.fit(np.hstack((x_data,x_data2)),y_data)
# hstack permet de créer un tableau de 7 lignes (observations) et 2 colonnes (attributs x et x^2)
print("régression: polynome degré 2")
print(lrp2.intercept_)
print(lrp2.coef_)
y_pred_lrp2=lrp2.predict(np.hstack((x_mod,x_mod2)))
# modèle polynomial avec $d=5$
lrp5 = lm.LinearRegression()
lrp5.fit(np.hstack((x_data,x_data2,x_data3,x_data4,x_data5)),y_data)
print("régresion: polynome degré 5")
print(lrp5.intercept_)
print(lrp5.coef_)
y_pred_lrp5=lrp5.predict(np.hstack((x_mod,x_mod2,x_mod3,x_mod4,x_mod5)))
Notez les valeurs des coefficients des deux modèles.
Représentation graphique:
plt.figure(figsize=(10,6))
plt.plot(x_mod, y_mod, '--k')
plt.plot(x_data, y_data, 'or', ms=10)
plt.plot(x_mod, y_pred_lr, '-g')
plt.plot(x_mod, y_pred_lrp2, '-b')
plt.plot(x_mod, y_pred_lrp5, '-c')
plt.xlim(0, 1.5)
plt.xlabel("x")
plt.ylim(-10, 80)
plt.ylabel("y")
plt.title('régression polynomiale')
plt.legend(["modèle","observations","régression linéaire","modèle degré 2","modèle degré 5"]);
print("régression polynomiale degré 2, MSE = %.3f" %mean_squared_error(y_data,lrp2.predict(np.hstack((x_data,x_data2)))))
print("régression polynomiale degré 5, MSE = %.3f" %mean_squared_error(y_data,lrp5.predict(np.hstack((x_data,x_data2,x_data3,x_data4,x_data5)))))
On voit que le modèle de degré 5 "colle" le mieux aux données mais est très mauvais pour prédire le modèle génératif. On dit qu'il y a surapprentissage, le modèle prédictif a de mauvaises capacités de généralisation. Par exemple, quelle valeur de $y$ est-elle prédite pour $x=1.2$?
A l'inverse, le modèle de la régression linéaire simple est dans une situation de sous-apprentissage, la MSE étant plutôt élevée dans ce cas.
Ces notions seront l'objet des cours suivants. On voit que le surapprentissage semble venir de la complexité du modèle de degré 5: il dépend de 6 paramètres donc il n'est pas étonnant qu'il puisse passer très près des 7 observations. Néanmoins, il ne faut pas oublier qu'à cause du bruit, il n'est sans doute pas judicieux d'utiliser un modèle prédictif si bien adapté aux observations.
Une manière classique de contrer le surapprentissage est de contraindre les paramètres du modèle prédictif (on dit qu'on régularise le modèle).
Dans le cadre de la régression, au lieu d'estimer les paramètres $\alpha_i$ par minimisation des moindres carrés, on peut chercher à minimiser: $$\sum_{i=1}^n \left|y_{data}[i] - \alpha_0 - \sum_{j=1}^d \alpha_j x_{data}[i]^j\right|^2 + C \sum_{j=0}^d \alpha_j^2$$ où $C$ est un paramètre positif (on parle d'hyperparamètre car il ne fait pas partie des paramètres estimés par minimisation de la fonction précédente).
On voit apparaître un compromis entre l'adéquation aux données et la valeur des paramètres $\alpha_i$. La régression linéaire classique correspond au cas particulier $C=0$.
Cette approche est la régression ridge.
Les cellules suivantes réalisent la régression ridge pour les modèles polynomiaux de degrés 2 et 5. On utilise les objets scikit_learn RidgeCV qui permettent d'estimer automatiquement une valeur optimale pour l'hyperparamètre $C$ par validation croisée (notion que l'on verra dans une prochaine séance).
ridge2 = lm.RidgeCV()
ridge2.fit(np.hstack((x_data,x_data2)),y_data)
print("ridge regression, polynome degré 2")
print(ridge2.intercept_)
print(ridge2.coef_)
y_pred_lrr2=ridge2.predict(np.hstack((x_mod,x_mod2)))
ridge5 = lm.RidgeCV()
ridge5.fit(np.hstack((x_data,x_data2,x_data3,x_data4,x_data5)),y_data)
print("ridge regression, polynome degré 5")
print(ridge5.intercept_)
print(ridge5.coef_)
y_pred_lrr5=ridge5.predict(np.hstack((x_mod,x_mod2,x_mod3,x_mod4,x_mod5)))
Comment évoluent les valeurs des paramètres?
Recherchez dans la documentation comment afficher la valeur de $C$ estimée automatiquement par RidgeCV.
plt.figure(figsize=(10,6))
plt.plot(x_mod, y_mod, '--k')
plt.plot(x_data, y_data, 'or', ms=10)
plt.plot(x_mod, y_pred_lr, '-g')
plt.plot(x_mod, y_pred_lrr2, '-b')
plt.plot(x_mod, y_pred_lrr5, '-c')
plt.xlim(0, 1.5)
plt.xlabel("x")
plt.ylim(-10, 80)
plt.ylabel("y")
plt.title('Ridge régression')
plt.legend(["modèle","observations","régression linéaire","régression ridge degré 2","régression ridge degré 5"])
plt.show()
print("régression ridge polynomiale degré 2, MSE = "+str(mean_squared_error(y_data,ridge2.predict(np.hstack((x_data,x_data2))))))
print("régression ridge polynomiale degré 5, MSE = "+str(mean_squared_error(y_data,ridge5.predict(np.hstack((x_data,x_data2,x_data3,x_data4,x_data5))))))
Que peut-on dire de l'évolution de la MSE, et de la capacité de généralisation?
Expérimentez le Lasso, décrit dans la documentation. Quelle est la différence essentielle avec la régression ridge? On utilisera LassoCV. Comment comprenez-vous la remarque "the Lasso regression yields sparse models, it can thus be used to perform feature selection" dans la documentation?