"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).\n",
"\n",
"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:\n",
"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).\n",
"\n",
"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$.\n",
# Introduction à l'apprentissage automatique: TP1 - Exercice 1
<br>
Le but de ce premier exercice est de se familiariser avec l'environnement __[Jupyter Notebook](https://jupyter-notebook.readthedocs.io/en/stable/)__ et l'utilisation de quelques fonctionnalités de la bibliothèque Python __[scikit-learn](http://scikit-learn.org)__ et de sa __[documentation](https://scikit-learn.org/stable/user_guide.html)__ (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 linéaire sous l'angle de l'apprentissage automatique. Rappelons que les problèmes de régression font partie de l'apprentissage supervisé. Vous avez étudié la régression linéaire 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.
<br>
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` ou `Kernel` 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](https://www.ibm.com/support/knowledgecenter/en/SSGNPV_2.0.0/dsx/markd-jupyter.html)__. Vous pouvez aussi examiner le code source des cellules de cette page en "double cliquant". Nous vous encourageons à écrire vos commentaires en couleur, comme dans l'exemple de la cellule suivante.
%% Cell type:markdown id: tags:
<fontcolor=green>
Voici un commentaire en couleur, et le code associé (à taper dans une cellule de type "Markdown", cf menu déroulant dans la barre en haut):
```
<font color=green>
Un texte en couleur: il faut utiliser des balises HTML.
</font>
```
</font>
<br>
%% Cell type:markdown id: tags:
## 1. On commence par charger les bibliothèques utiles
%% Cell type:code id: tags:
``` python
# import des bibliothèques Python utiles:
importnumpyasnp
importsklearn.linear_modelaslm
importsklearn.metricsasmetrics
importsklearn.preprocessingaspreprocessing
importmatplotlib.pyplotasplt
# "magic function" Jupyter pour l'affichage des graphiques dans le carnet de manière interactive (zoom possible par exemple):
# (voir plus loin)
%matplotlibnotebook
# si on veut supprimer les interactions (si elles ralentissent trop l'affichage par exemple), il faut faire:
# si on veut supprimer les interactions (si elles ralentissent trop l'affichage par exemple, ou si elles ne fonctionnent pas), il faut faire:
# %matplotlib inline
```
%% Cell type:markdown id: tags:
## 2. Données: distance de freinage en fonction de la vitesse d'un véhicule
<br>
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 le même dossier que ce carnet.
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.
%% Cell type:code id: tags:
``` python
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))
# représentation graphique
plt.figure(figsize=(9,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 ; final permet d'éviter un affichage intempestif dans le carnet
plt.ylim(0,100)
plt.show()
```
%% Cell type:markdown id: tags:
<br>
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})$.
<br>
## 3. Régression linéaire (rappel)
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
sur l'ensemble des $n=20$ observations $(x_{data}[i],y_{data}[i])_{1\leq i\leq n}$.
%% Cell type:code id: tags:
``` python
# 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:
# On estime les paramètres a_0 et a_1 à l'aide de la méthode fit:
# (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, estimation par la méthodes des moindres carrés)
# Ils sont stockés dans les attributs suivants de l'objet lr:
# Les paramètres sont stockés dans les attributs suivants de l'objet lr:
print(lr.intercept_)
print(lr.coef_)
```
%% Cell type:markdown id: tags:
__Question 1__. Quelles sont les valeurs de $a_0$ et $a_1$ ? (voir la [documentation](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html))
%% Cell type:markdown id: tags:
<fontcolor=red>
_Votre réponse._
</font>
%% Cell type:markdown id: tags:
La cellule suivante représente graphiquement la prédiction par régression linéaire.
%% Cell type:code id: tags:
``` python
# On "prédit" les valeurs de y pour 30 valeurs de x réparties régulièrement 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 défini précédemment)
# représentation graphique
plt.figure(figsize=(9,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"]);
plt.legend(["observations","droite de régression"])
plt.show()
```
%% Cell type:markdown id: tags:
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?
%% Cell type:markdown id: tags:
<fontcolor="red">
_Votre réponse._
</font>
%% Cell type:markdown id: tags:
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)
Nous allons à présent construire des modèles prédictifs polynomiaux:
$$ y_{pred} = a_0 + \sum_{j=1}^d a_j x^j$$
pour des entiers $d>1$. On se limitera à deux modèles: $d=2$ et $d=6$.
Autrement dit, on prédit la distance d'arrêt comme une fonction polynomiale de la vitesse.
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 la $i$-ème observation est le vecteur $x[i],x[i]^2,x[i]^3,\dots, x[i]^d$.
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:
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$.
Rappelons deux points discutés dans le polycopié:
- le paramètre $a_0$ n'est pas inclus dans la régularisation;
- comme les paramètres $a_1,\dots, a_j$ ont le même poids 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. Dans le cas contraire, la régularisation n'aurait pas le même effet sur chaque caractéristique.
Cette méthode est la régression _ridge_.
<br>
Dans un premier temps, on normalise les observations en divisant chaque vecteur de caractéristique par son écart-type empirique.
%% Cell type:code id: tags:
``` python
# on définit un objet "scaler6" sur les données de X_data6
# les coefficients suivants sont les quantités par lesquelles sont divisées les 6 caractéristiques:
print(scaler6.scale_)
# on constate que ces coefficients sont bien l'écart-type de la caractéristique:
print(np.std(X_data6,axis=0))
# la procédure suivante permet de normaliser les caractéristiques de X_data6 et X6 en les divisant par leur écart-type:
X_data6_n=scaler6.transform(X_data6)
X6_n=scaler6.transform(X6)
# remarquons que l'on ne normalise pas les Y
```
%% Cell type:markdown id: tags:
Nous allons estimer les paramètres $a'_j$ de modèles de régression sur les données normalisées (où la $j$-ème caractéristique du dataset de départ est divisée par $\sigma_j$). Si $x$ est une vitesse en km/h, il faut donc diviser $x^j$ par $\sigma_j$ pour prédire la distance d'arrêt, selon:
La prédiction s'exprime donc en fonction de la vitesse en km/h de la manière suivante:
$$ y_{pred} = a_0 + \sum_{j=1}^d a_j x^j$$
où $a_0=a'_0$, et $\;\forall 1\leqslant j \leqslant d, \, a_j=a'_j/\sigma_j$.
Dans la suite, ce sont ces coefficients $a_j$ que l'on affichera pour pouvoir comparer aux coefficients obtenus avec la régression simple (section 3).
<br>
L'influence du paramètre $\alpha$ est illustrée par la cellule suivante, dans le cas de la régression polynomiale de degré 6: on effectue des regressions ridges pour différentes valeurs de $\alpha$.
Notez que vous retrouvez bien les valeurs de la régression linéaire classique pour $\alpha=0$.
%% Cell type:code id: tags:
``` python
ridgealpha0=lm.Ridge(alpha=0)
ridgealpha0.fit(X_data6_n,Y_data)
print("ridge regression alpha=0")
print(ridgealpha0.intercept_)
print(ridgealpha0.coef_/scaler6.scale_)# on affiche bien les a_j, pas les a'_j
Y_pred_ridgealpha0=ridgealpha0.predict(X6_n)# on prédit avec les a'_j, donc sur les données normalisées
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$ le plus adapté au problème.
Un _hyperparamètre_ est un paramètre du modèle qui n'est pas déterminé par apprentissage mais qui doit être fixé par l'utilisateur ou de manière à optimiser un critère.
<br>
La méthode `RidgeCV` (à utiliser à la place de `Ridge`) permet de faire une régression ridge en sélectionnant un hyperparamètre par _validation croisée_. **Dans cette question, nous l'utilisons en "boîte noire"**.
Une partie de la prochaine séance sera consacrée à cette méthode de sélection de modèle.
<br>
On commence par normaliser les datasets pour la régression linéaire et polynomiale de degré 2, comme précédemment pour la régression polynomiale de degré 6.
Par exemple, pour la régression ridge dont la valeur de $\alpha$ est choisie par validation croisée parmi 11 valeurs de la forme $10^{-i}$ avec $i$ prenant 11 valeurs régulièrement espacées entre entre -5 et 5:
%% Cell type:code id: tags:
``` python
ridge1=lm.RidgeCV(alphas=np.logspace(-5,5,11))
ridge1.fit(X_data_n,Y_data)
print("ridge regression, polynome degré 1")
print(ridge1.intercept_)
print(ridge1.coef_/scaler1.scale_)
print("alpha sélectionné: %.5f"%ridge1.alpha_)
```
%% Cell type:markdown id: tags:
Avec des données normalisées, il est d'usage de chercher $\alpha$ autour de la valeur 1 comme ici.
<br>
La cellule suivante réalise la régression _ridge_ avec sélection de l'hyperparamètre $\alpha$ par validation croisée, dans les cas $d=1$, $d=2$, et $d=6$.
__Question 6__. 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.
%% Cell type:markdown id: tags:
<fontcolor=red>
_Votre réponse._
</font>
%% Cell type:markdown id: tags:
# 5. Lasso
Expérimentez le _Lasso_ , décrit __[dans la documentation](https://scikit-learn.org/stable/modules/linear_model.html#lasso)__.
__Question 7__. Quelle est la différence essentielle avec la régression ridge? On utilisera __[LassoCV](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LassoCV.html#sklearn.linear_model.LassoCV)__ qui sélectionne l'hyperparamètre. Comment comprenez-vous la remarque _"As the Lasso regression yields sparse models, it can thus be used to perform feature selection"_ dans la documentation?
_Remarque_ :
observez dans la documentation que la stratégie pour fixer l'hyperparamètre $\alpha$ n'est pas la même que pour `RidgeCV`
%% Cell type:markdown id: tags:
<fontcolor=red>
_Votre réponse._
</font>
%% Cell type:code id: tags:
``` python
lasso1=lm.LassoCV()
lasso1.fit(X_data_n,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_/scaler1.scale_)
print("alpha sélectionné: %.5f"%lasso1.alpha_)
lasso2=lm.LassoCV()
lasso2.fit(X_data2_n,np.ravel(Y_data))
print("\nlasso regression, polynome degré 2")
print(lasso2.intercept_)
print(lasso2.coef_/scaler2.scale_)
print("alpha sélectionné: %.5f"%lasso2.alpha_)
lasso6=lm.LassoCV()
lasso6.fit(X_data6_n,np.ravel(Y_data))
print("\nlasso regression, polynome degré 6")
print(lasso6.intercept_)
print(lasso6.coef_/scaler6.scale_)
print("alpha sélectionné: %.5f"%lasso6.alpha_)
```
%% Cell type:markdown id: tags:
__Question 8__. Que dire des coefficients calculés?
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](https://fr.wikipedia.org/wiki/Distance_d%27arr%C3%AAt)__.
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:
%% Cell type:code id: tags:
``` python
print("a_0 = 0")
print("a_1 = %.5f"%(1/3.6))
print("a_2 = %.5f"%(1/(3.6**2*2*10)))
```
%% Cell type:markdown id: tags:
__Question 9__. Comparez aux modèles basés sur les données obtenus précédemment.
