Introduction à l'apprentissage automatique: TP1 - Exercice 1


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 et de sa documentation (utilisez la recherche en haut à droite sur cette dernière page). Tous les modèles d'apprentissage de scikit-learn s'utilisent selon le même principe: le TP aujourd'hui est représentatif de notre travail dans les prochaines séances.

Nous allons nous intéresser à un problème de régression sous l'angle des problématiques de l'apprentissage automatique. Rappelons que les problèmes de régression font partie de l'apprentissage supervisé. Vous avez étudié la régression en cours d'analyse de données, nous introduisons à présent la notion de régularisation qui interviendra à plusieurs reprises dans la suite de ce cours.


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 Noyau 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". Nous vous encourageons à écrire vos commentaires en couleur, comme dans l'exemple de la cellule suivante.

Voici un commentaire en couleur, et le code associé (à taper dans une cellule de type "Markdown", cf onglet dans la barre en haut):

<font color=green>
Un texte en couleur: il faut utiliser des balises HTML.
</font>

</font>

On commence par charger les bibliothèques utiles

In [ ]:
# import des bibliothèques Python utiles:
import numpy as np
import sklearn.linear_model as lm
import sklearn.metrics as metrics 
import matplotlib.pyplot as plt

# "magic function" Jupyter pour l'affichage des graphiques dans le carnet:
%matplotlib inline

Données: distance de freinage en fonction de la vitesse d'un véhicule


L'objectif de cet exercice est de prédire la distance d'arrêt d'un véhicule en fonction de sa vitesse, dans une approche "basée données".

Commencez par sauvegarder le fichier freinage.txt dans votre espace de travail.

Dans la cellule ci-dessous, on charge les données (les observations ): Y_data représente la distance de freinage d'un véhicule en fonction de la vitesse X_data du véhicule.

Remarque: pour des raisons pédagogiques, il s'agit de données générées selon un modèle physique, et pas de données expérimentales réelles.

In [ ]:
data = np.loadtxt('freinage.txt')  # on lit les données dans un fichier txt
print("Les observations :\n")
print(data)

# les modules sklearn demandent que les données soient représentées par des vecteurs colonnes
X_data = data[:,0].reshape(len(data),1)  # première colonne (sans reshape, X_data serait un vecteur ligne)
Y_data = data[:,1].reshape(len(data),1)  # deuxième colonne
print("\nnombre d'observations : %d" %len(X_data))


plt.figure(figsize=(10,6))
plt.scatter(X_data, Y_data)
plt.xlabel("X: vitesse (km/h)")
plt.ylabel("Y: distance d'arrêt (m)")
plt.grid()
plt.xlim(0, 130)
plt.ylim(0, 100);  # le ; permet d'éviter un affichage intempestif dans le carnet


Nous allons à présent créer des modèles permettant de prédire les valeurs $y$ correspondant à des valeurs de $x$ entre 0 et 130 km/h. Ces modèles basés sur les données seront créés à partir des 20 valeurs $(x_{data},y_{data})$.


Régression linéaire

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} = a_0 + a_1 x$$

Les valeurs de $a_0$ et $a_1$ sont estimées par la méthode des moindres carrés: on cherche les paramètres $a_0$ et $a_1$ qui minimisent $$\sum_{i=1}^n \left|y_{data}[i] - a_0 - a_1 x_{data}[i]\right|^2$$ sur l'ensemble des $n=20$ observations $(x_{data}[i],y_{data}[i])_{1\leq i\leq n}$.

In [ ]:
# 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 
# voir les détails dans la documentation de LinearRegression: on ne se servira que de quelques uns d'entre eux

# On estime les paramètres a_0 et a_1: 
# (remarque: sklearn attend des données sous forme de vecteurs colonnes)
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 ici dans les attributs suivants:
print(lr.intercept_)  
print(lr.coef_)

Question 1. Quelles sont les valeurs de $a_0$ et $a_1$ ? (voir la documentation)

La cellule suivante représente graphiquement la prédiction par régression linéaire.

In [ ]:
# On "prédit" les valeurs de y pour 30 valeurs de x comprises entre 0 et 130
X = np.linspace(0,130,num=30).reshape(30,1)
Y_pred_lr = lr.predict(X)  
# la méthode predict permet de prédire les valeurs y pour les valeurs de x passées en argument
# (à l'aide du modèle lr)

plt.figure(figsize=(10,6))
plt.plot(X_data, Y_data,'o')
plt.plot(X, Y_pred_lr, '-g')
plt.xlim(0, 130)
plt.ylim(-10, 100)
plt.xlabel("X: vitesse (km/h)")
plt.ylabel("Y: distance d'arrêt (m)")
plt.grid()
plt.title('observations et régression linéaire')
plt.legend(["observations","droite de régression"]);

Comme on le voit, la droite de régression passe globalement entre les observations.

Question 2. Ce modèle vous semble-t-il réaliste?

La cellule suivante permet d'afficher l'erreur quadratique moyenne de prédiction sur les observations.

Question 3. Comment cet indicateur est-il défini? (cherchez dans la documentation)

In [ ]:
print("régression linéaire: MSE = "+str(metrics.mean_squared_error(Y_data,lr.predict(X_data))))


Nous allons à présent construire des modèles prédictifs polynomiaux: $$ y_{pred} = a_0 + \sum_{i=1}^d a_d x^d$$ pour des entiers $d>1$.

Les coefficients $a_0, \dots, a_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 comme dans le cours d'analyse de données: au lieu d'utiliser comme régresseur la seule variable $x$, on utilise également les variables $x^2, x^3,\dots, x^d$. Tout se passe comme si les observations sont les vecteurs $(x,x^2,x^3,\dots, x^d)$.

In [ ]:
# création des vecteurs "puissance"
from sklearn.preprocessing import PolynomialFeatures
poly2 = PolynomialFeatures(degree=2,include_bias=False)
X_data2 = poly2.fit_transform(X_data)
poly6 = PolynomialFeatures(degree=6,include_bias=False)
X_data6 = poly6.fit_transform(X_data)
print("les cinq premières observations dans X_data6: \n(les colonnes contiennent les puissances successives de la première)")
print(X_data6[:5,:])

# modèle polynomial avec d=2
lrp2 = lm.LinearRegression()
lrp2.fit(X_data2,Y_data)
print("\nrégression polynomiale degré 2")
print(lrp2.intercept_)
print(lrp2.coef_)

# modèle polynomial avec d=6
lrp6 = lm.LinearRegression()
lrp6.fit(X_data6,Y_data)
print("\nrégresion polynomiale degré 6")
print(lrp6.intercept_)
print(lrp6.coef_)

Question 4. Comment s'écrivent les deux modèles de prédiction de $y$ en fonction de $x$? (ne reportez que quelques chiffres significatifs)

La représentation graphique des différents modèle est donnée par la cellule suivante.

In [ ]:
X2=poly2.fit_transform(X)
X6=poly6.fit_transform(X)
Y_pred_lrp2=lrp2.predict(X2)  # prédiction du modèle sur les valeurs X définies plus haut, pour le graphique
Y_pred_lrp6=lrp6.predict(X6)

plt.figure(figsize=(10,6))
plt.plot(X_data, Y_data,'o')
plt.plot(X, Y_pred_lr, '-g')
plt.plot(X, Y_pred_lrp2, '-b')
plt.plot(X, Y_pred_lrp6, '-c')
plt.xlim(0, 130)
plt.ylim(-10, 100)
plt.xlabel("X: vitesse (km/h)")
plt.ylabel("Y: distance d'arrêt (m)")
plt.grid()
plt.title('régression')
plt.legend(["observations","régression linéaire","modèle degré 2","modèle degré 6"]);

Question 5. Comment interprêter ces résultats, ainsi que les calculs de MSE ci-dessous, dans le vocabulaire de l'apprentissage automatique?

In [ ]:
print("régression linéaire: MSE = %.3f" %metrics.mean_squared_error(Y_data,lr.predict(X_data)))

print("régression polynomiale degré 2, MSE = %.3f" %metrics.mean_squared_error(Y_data,lrp2.predict(X_data2)))

print("régression polynomiale degré 6, MSE = %.3f" %metrics.mean_squared_error(Y_data,lrp6.predict(X_data6)))

Régression ridge

Une manière classique de contrer le surapprentissage est de contraindre les paramètres du modèle prédictif à prendre des valeurs "pas trop grandes" (on dit qu'on régularise le modèle).

Dans le cadre de la régression, au lieu d'estimer les paramètres $a_i$ par minimisation des moindres carrés, on peut chercher à minimiser: $$\sum_{i=1}^n \left|y_{data}[i] - a_0 - \sum_{j=1}^d a_j x_{data}[i]^j\right|^2 + \alpha \sum_{j=1}^d a_j^2$$ où $\alpha$ est un paramètre positif, fixé a priori par l'utilisateur (on parle d'hyperparamètre car $\alpha$ 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 (mesurée par la MSE, premier terme de l'expression) et la valeur des paramètres $a_j$ (qui interviennent par le carré de la norme euclidienne du vecteur $(a_1,\dots,a_d) $). On remarque que la régression linéaire classique correspond au cas particulier $\alpha=0$.

Notons deux points discutés dans le polycopié:

  • le coefficient $a_0$ n'est pas inclus dans la régularisation;
  • il vaut mieux normaliser les caractéristiques (les composantes des observations) de manière à ce qu'elles varient dans le même intervalle, sinon la régularisation n'aurait pas le même effet sur chacune d'entre elles.

Cette méthode est la régression ridge.


L'influence du paramètre $\alpha$ est illustrée par la cellule suivante, dans le cas d'une régression polynomiale.

Remarquez normalize=True, et notez que vous retrouvez bien les valeurs de la régression linéaire classique pour $\alpha=0$.

In [ ]:
ridgealpha0 = lm.Ridge(normalize=True,alpha=0)
ridgealpha0.fit(X_data6,Y_data)
print("ridge regression alpha=0")
print(ridgealpha0.intercept_)
print(ridgealpha0.coef_)
Y_pred_ridgealpha0 = ridgealpha0.predict(X6)

ridgealpha01 = lm.Ridge(normalize=True,alpha=0.1)
ridgealpha01.fit(X_data6,Y_data)
print("\nridge regression alpha=0.1")
print(ridgealpha01.intercept_)
print(ridgealpha01.coef_)
Y_pred_ridgealpha01 = ridgealpha01.predict(X6)

ridgealpha1 = lm.Ridge(normalize=True,alpha=1)
ridgealpha1.fit(X_data6,Y_data)
print("\nridge regression alpha=1")
print(ridgealpha1.intercept_)
print(ridgealpha1.coef_)
Y_pred_ridgealpha1 = ridgealpha1.predict(X6)

ridgealpha10 = lm.Ridge(normalize=True,alpha=10)
ridgealpha10.fit(X_data6,Y_data)
print("\nridge regression alpha=10")
print(ridgealpha10.intercept_)
print(ridgealpha10.coef_)
Y_pred_ridgealpha10 = ridgealpha10.predict(X6)

ridgealpha100 = lm.Ridge(normalize=True,alpha=100)
ridgealpha100.fit(X_data6,Y_data)
print("\nridge regression alpha=100")
print(ridgealpha100.intercept_)
print(ridgealpha100.coef_)
Y_pred_ridgealpha100 = ridgealpha100.predict(X6)

plt.figure(figsize=(10,6))
plt.plot(X_data, Y_data,'o')
plt.plot(X, Y_pred_ridgealpha0, '-g')
plt.plot(X, Y_pred_ridgealpha01, '-b')
plt.plot(X, Y_pred_ridgealpha1, '-c')
plt.plot(X, Y_pred_ridgealpha10, '-r')
plt.plot(X, Y_pred_ridgealpha100, '-k')
plt.xlim(0, 130)
plt.ylim(-10, 80)
plt.xlabel("X: vitesse (km/h)")
plt.ylabel("Y: distance d'arrêt (m)")
plt.grid()
plt.title('régression ridge, d=6')
plt.legend(["observations","alpha=0","alpha=0.1","alpha=1","alpha=10","alpha=100"]);

Notez que les coefficients $a_j$ ($1\leq j \leq 6$) diminuent et que $a_0$ tend vers la moyenne des $Y_\text{data}[i]$ lorsque $\alpha$ augmente. Il faut maintenant choisir l'hyperparamètre $\alpha$.


Question 6. que fait RidgeCV lorsqu'on l'appelle par lm.RidgeCV()? (voir la documentation: on ne vous demande que le principe général, pas le rôle de chaque paramètre)

La cellule suivante réalise la régression ridge avec sélection de l'hyperparamètre $\alpha$ par validation croisée. Ici, on cherche dans une gamme de valeurs de $\alpha$ plus large que celle utilisée par défaut: quelle est-elle exactement? Notez qu'on normalise les caractéristiques et qu'on utilise une validation croisée à 5 plis. Après normalisation, il est d'usage de chercher $\alpha$ autour de la valeur 1.

In [ ]:
ridge1 = lm.RidgeCV(normalize=True, alphas=np.logspace(-5, 5, 20), cv=5)
ridge1.fit(X_data,Y_data)
print("ridge regression, polynome degré 1")
print(ridge1.intercept_)
print(ridge1.coef_)
print("alpha sélectionné: %.5f" %ridge1.alpha_)

ridge2 = lm.RidgeCV(normalize=True, alphas=np.logspace(-5, 5, 20), cv=5)
ridge2.fit(X_data2,Y_data)
print("\nridge regression, polynome degré 2")
print(ridge2.intercept_)
print(ridge2.coef_)
print("alpha sélectionné: %.5f" %ridge2.alpha_)

ridge6 = lm.RidgeCV(normalize=True, alphas=np.logspace(-5, 5, 20), cv=5)
ridge6.fit(X_data6,Y_data)
print("\nridge regression, polynome degré 6")
print(ridge6.intercept_)
print(ridge6.coef_)
print("alpha sélectionné: %.5f" %ridge6.alpha_)
In [ ]:
Y_pred_lrr1=ridge1.predict(X)
Y_pred_lrr2=ridge2.predict(X2)
Y_pred_lrr6=ridge6.predict(X6)

plt.figure(figsize=(10,6))
plt.plot(X_data, Y_data,'o')
plt.plot(X, Y_pred_lrr1, '-g')
plt.plot(X, Y_pred_lrr2, '-b')
plt.plot(X, Y_pred_lrr6, '-c')
plt.xlim(0, 130)
plt.ylim(-10, 80)
plt.xlabel("X: vitesse (km/h)")
plt.ylabel("Y: distance d'arrêt (m)")
plt.grid()
plt.title('régression ridge')
plt.legend(["observations","modèle degré 1","modèle degré 2","modèle degré 6"]);
In [ ]:
print("régression ridge polynomiale degré 1, MSE = %.2f" %metrics.mean_squared_error(Y_data,ridge1.predict(X_data)))

print("régression ridge polynomiale degré 2, MSE = %.2f" %metrics.mean_squared_error(Y_data,ridge2.predict(X_data2)))

print("régression ridge polynomiale degré 6, MSE = %.2f" %metrics.mean_squared_error(Y_data,ridge6.predict(X_data6)))

Question 7. Comment évoluent les valeurs des paramètres par rapport à la régression classique? Constatez que la régularisation permet de définir des modèles moins affectés par le surapprentissage.

Lasso

Expérimentez le Lasso , décrit dans la documentation.

Question 8. Quelle est la différence essentielle avec la régression ridge? On utilisera LassoCV. Comment comprenez-vous la remarque "As the Lasso regression yields sparse models, it can thus be used to perform feature selection" dans la documentation?

Remarques :

  • observez dans la documentation que la stratégie pour fixer l'hyperparamètre $\alpha$ n'est pas la même que pour RidgeCV
  • on change la valeur de max_iter pour éviter un problème de convergence dans le dernier modèle (essayez sans cette option)
In [ ]:
lasso1 = lm.LassoCV(normalize=True)
lasso1.fit(X_data,np.ravel(Y_data))  # le lasso attend un 1D array pour y, d'où "np_ravel"
print("lasso regression, fonction affine")
print(lasso1.intercept_)
print(lasso1.coef_)
print("alpha sélectionné: %.5f" %lasso1.alpha_)

lasso2 = lm.LassoCV(normalize=True)
lasso2.fit(X_data2,np.ravel(Y_data))  
print("\nlasso regression, polynome degré 2")
print(lasso2.intercept_)
print(lasso2.coef_)
print("alpha sélectionné: %.5f" %lasso2.alpha_)

lasso6 = lm.LassoCV(normalize=True, max_iter=10000)  # (max_iter pour éviter des problèmes numériques ensuite)  
lasso6.fit(X_data6,np.ravel(Y_data))
print("\nlasso regression, polynome degré 6")
print(lasso6.intercept_)
print(lasso6.coef_)
print("alpha sélectionné: %.5f" %lasso6.alpha_)

Question 9. Que dire des coefficients calculés?

In [ ]:
Y_pred_lasso1=lasso1.predict(X)
Y_pred_lasso2=lasso2.predict(X2)
Y_pred_lasso6=lasso6.predict(X6)

plt.figure(figsize=(10,6))
plt.plot(X_data, Y_data,'o')
plt.plot(X, Y_pred_lr, '-g')
plt.plot(X, Y_pred_lasso2, '-b')
plt.plot(X, Y_pred_lasso6, '-c')
plt.xlim(0, 130)
plt.ylim(-10, 80)
plt.xlabel("X: vitesse (km/h)")
plt.ylabel("Y: distance d'arrêt (m)")
plt.grid()
plt.title('Lasso')
plt.legend(["observations","régression linéaire","Lasso modèle degré 2","Lasso modèle degré 6"]);

Les physiciens nous apprennent qu'un modèle réaliste de la distance d'arrêt $D$ en fonction de la vitesse initiale $v_0$ est donné par la relation: $$ D = v_0 t_r + \frac{v_0^2}{2a}$$

où $a$ est la décélération (supposée constante), et $t_r$ le temps de réaction.

Voir cette page wikipedia.

Les données de cet exercice ont été générées en supposant un temps de réaction de une seconde et une décélération de $10 m.s^{-2}$, ce qui aboutit à cette équation: $$Y=X/3.6+(X/3.6)^2/(2*10)$$

Un bruit gaussien d'écart-type grandissant avec X a été ensuite ajouté à Y.

On devrait donc obtenir un modèle de degré 2 avec pour valeurs de paramètres:

In [ ]:
print("a_0 = 0")
print("a_1 = %.5f" %(1/3.6))
print("a_2 = %.5f" %(1/(3.6**2*2*10)))

Question 10. Comparez aux modèles basés sur les données obtenus précédemment.

Question 11. En utilisant la fonction cross_val_score expliquée sur cette page et en vous inspirant de l'exemple ci-dessous pour la régression linéaire, déterminez quel est le meilleur modèle parmi ceux que vous avez envisagés (c'est-à-dire lr, lrp2, lrp6, ridge1, ridge2, ridge6, lasso1, lasso2, lasso6). Utilisez la même syntaxe que dans fit pour spécifier les données en argument de cross_val_score. Par défaut, le score calculé pour tous ces modèles est le coefficient de détermination $R^2$: plus il est proche de 1, meilleur est le modèle sur le pli de test. On voit d'après la définition sur la page wikipedia que si l'écart entre la vraie valeur et la valeur prédite sur le pli de test est grand, alors $R^2$ peut être négatif. C'est le cas en particulier s'il y a surapprentissage.

Remarque : on utilise une validation croisée à 5 plis (valeur par défaut dans la dernière version de scikit-learn) par cohérence avec ce que l'on a utilisé pour sélectionner l'hyperparamètre de la régression ridge et du lasso. Néanmoins, le critère est ici le coefficient de détermination, ce qui n'est pas le cas dans ̀RidgeCV() et LassoCV().

In [ ]:
# votre code ici:

from sklearn.model_selection import cross_val_score

scores_lr = cross_val_score(lr, X_data, Y_data, cv=5)
print("lr - Accuracy: %0.2f (+/- %0.2f)" % (scores_lr.mean(), scores_lr.std() * 2))
In [ ]: