TP3. Cross-validation, boosting, hyperparameter tuning¶
Dans ce TP, on va voir des aspects plus "méta" du machine learning. On tentera de répondre aux questions suivantes :
- Comment évaluer de manière robuste les modèles, sans que cela ne dépende totalement de la séparation train/test? C'est ce que la cross validation nous aidera à faire dans un premier temps.
- Comment aggréger plein de modèles faibles pour en faire un modèle fort? On fera du boosting.
- Comment trouver les meilleurs hyperparamètres d'un modèle? On verra l'algorithme GridSearch
Dans ce notebook, les éléments explicatifs seront écrits au fil des questions, pour suivre un enchaînement plus naturel.
Consignes¶
Cross-validation
- Charger les datasets
irisetdiabetesde scikit-learn - Séparer le dataset
irisen deux train/test sets différents,(X_train1, y_train1, X_test1, y_test1), (X_train2, y_train2, X_test2, y_test2). Afficher l'histogramme des espèces pour chacun des train/test. - Faire de même pour
diabetes. Que remarquez-vous? - Entraîner une forêt aléatoire
rf1surX_train1, y_train1, et une autre forêt aléatoirerf2surX_train2, y_train2. Evaluer. Comparer. - Faire une cross-validation, pour
irisetdiabetes.
Boosting
- Evaluer les performances d'adaboost, avec le cross-val score, sur
irisetdiabetes - Les performances changent-elles en fonction du
base_estimatorchoisi?
Hyperparameter tuning
- Avec GridSearch
Cross-validation¶
1. Charger les datasets iris et diabetes de scikit-learn¶
Toujours avec les fonctions de load de scikit-learn:
from sklearn.datasets import load_iris, load_diabetes
iris = load_iris()
X_iris, y_iris = iris.data, iris.target
diabetes = load_diabetes()
X_diabetes, y_diabetes = diabetes.data, diabetes.target
2. Séparer le dataset iris en deux train/test sets différents, (X_train1, y_train1, X_test1, y_test1), (X_train2, y_train2, X_test2, y_test2). Afficher l'histogramme des espèces pour chacun des train/test.¶
On reprend, comme dans les sujets précédents, la fonction train_test_split de scikit-learn
from sklearn.model_selection import train_test_split
X_train1, X_test1, y_train1, y_test1 = train_test_split(X_iris, y_iris)
X_train2, X_test2, y_train2, y_test2 = train_test_split(X_iris, y_iris)
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
plt.figure()
plt.hist(y_test1, align='left', label="split 1")
plt.hist(y_test2, align='right', label="split 2")
plt.legend()
plt.show()
3. Faire de même pour diabetes. Que remarquez-vous?¶
Idem en régression :
X_train1, X_test1, y_train1, y_test1 = train_test_split(X_diabetes, y_diabetes)
X_train2, X_test2, y_train2, y_test2 = train_test_split(X_diabetes, y_diabetes)
plt.figure()
plt.hist(y_test1, align='left', label="split 1")
plt.hist(y_test2, align='right', label="split 2")
plt.legend()
plt.show()
On remarque que selon le train-test split, les différents X_train, X_test peuvent changer, à cause du tirage aléatoire. Cela a pour conséquence que pour deux train-test split différents, on a des datasets d'entraînement/test complètement différents, même si l'on ajoute l'option stratify de scikit-learn (voir TP1).
En évaluation, on aura des répartitions des classes/valeurs différentes pour deux y_test : mécaniquement, on aura des métriques d'évaluation différentes selon le tirage aléatoire dans les train-test split.
On le voit clairement dans l'expérimentation suivante :
4. Entraîner une forêt aléatoire rf1sur X_train1, y_train1, et une autre forêt aléatoire rf2sur X_train2, y_train2. Evaluer. Comparer.¶
from sklearn.ensemble import RandomForestClassifier
X_train1, X_test1, y_train1, y_test1 = train_test_split(X_iris, y_iris)
X_train2, X_test2, y_train2, y_test2 = train_test_split(X_iris, y_iris)
rf1 = RandomForestClassifier()
rf2 = RandomForestClassifier()
rf1.fit(X_train1, y_train1) # Notre première forêt
rf2.fit(X_train2, y_train2) # La deuxième
from sklearn.metrics import classification_report
print(classification_report(y_test1, rf1.predict(X_test1)))
print(classification_report(y_test2, rf1.predict(X_test2)))
Sur le dataset iris (classification), on voit que le même modèle (une forêt aléatoire en classification), entraîné et évalué sur deux splits différents, peut passer d'une accuracy de 97% à une accuracy de 100%.
Regardons ce qu'il en est sur la régression:
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import r2_score
X_train1, X_test1, y_train1, y_test1 = train_test_split(X_diabetes, y_diabetes)
X_train2, X_test2, y_train2, y_test2 = train_test_split(X_diabetes, y_diabetes)
rf1 = RandomForestRegressor()
rf2 = RandomForestRegressor()
rf1.fit(X_train1, y_train1)
rf2.fit(X_train2, y_train2)
print(r2_score(y_test1, rf1.predict(X_test1)))
print(r2_score(y_test2, rf1.predict(X_test2)))
Ici, c'est encore pire : le R2 (qui varie entre 0 et 1 et qui est proche de 1 pour les bons modèles et proche de 0 pour les mauvais) change quasiment du simple au double.
On veut avoir une estimation plus robuste des performances du modèle, et c'est là que la cross-validation intervient.
On appelle une validation croisée "k-fold" l'algorithme suivant :
Diviser le dataset en $k$ sous-datasets de taille égale ($\sim \frac{n}{k}$ lignes dans chaque subset). On appelle $$g : \{1, ..., N \} \rightarrow \{1, ..., k\}$$ la fonction qui à un $x_i$ lui associe sa partition (1, 2, ... ou k)
Pour i=1, ...,k, on entraîne un modèle sur toutes les données sauf la partition $i$, qu'on appelle $$\hat{f}^{-i}(x) $$
Une fois que l'on a entraîné nos $k$ modèles, on peut estimer l'erreur par: $$ \text{CV}(\hat{f}) = \frac{1}{N} \sum_{i=1}^N L\left(y_i, \hat{f}^{-g(i)}(x_i) \right)$$
Où $L$ est la fonction de perte/d'évaluation choisie.
Dans la suite, on utilisera uniquement cette estimation consolidée pour évaluer nos algorithmes. Regardons plus en détail comment la cross validation fonctionne avec sklearn.
5. Faire une cross-validation, pour iris et diabetes.¶
La fonction cross_validate renvoie les scores de cross-validation, mais aussi des informations sur le temps nécessaire pour fitter et évaluer chacun des modèles.
Dans la pratique, ces informations ne sont pas toujours nécessaires, on utilise donc plutôt en général la fonction cross_val_score.
from sklearn.model_selection import cross_validate
cross_validate(RandomForestClassifier(), X_iris, y_iris )
from sklearn.model_selection import cross_val_score
# cross_val_score renvoie un array qui contient les scores de chacun des k=5 modèles
# On calcule nous même la moyenne avec le .mean():
cross_val_score(RandomForestClassifier(), X_iris, y_iris, cv=5).mean()
Maintenant qu'on a vu comment utiliser les deux fonctions, regardons la fluctuation du cross val score (il fluctue encore légèrement : le tirage est aléatoire), en fonction du paramètre $k$.
Ecrire une fonction qui fait $100$ fois une cross validation pour une forêt aléatoire sur le dataset iris, avec comme paramètre $k$.
Tracer des boxplot des séries de $100$ estimations pour $k=5$ et pour $k=10$.
Tracer la courbe d'évolution de la variance des séries de $100$ estimations pour $k=1, ..., 100$.
Que remarquez-vous?
def see_fluc(k):
fluctuation = []
for i in range(100):
fluctuation.append(cross_val_score(RandomForestClassifier(), X_iris, y_iris, cv=k).mean())
return fluctuation
fluctuation1 = see_fluc(5)
fluctuation2 = see_fluc(10)
plt.figure()
plt.boxplot([fluctuation1, fluctuation2])
plt.show()
Les boxplots ne nous permettent pas de très bien voir la différence de fluctuation.
Traçons la courbe de l'évolution de la variance des séries en fonction de $k$: (un peu long, on entraîne beaucoup de modèles)
k = [2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 25]
v = [np.array(see_fluc(i)).std() for i in k]
plt.figure()
plt.plot(k, v)
plt.xlabel('k')
plt.ylabel('Variance pour 100 k-fold CV')
plt.figure()
Boosting¶
Le boosting est une technique d'apprentissage qui consiste à aggréger plusieurs prédicteurs faibles, pour construire un prédicteur plus performant.
Par soucis de clarté des notations, on considèrera dans la suite un problème de classification binaire tel que les $n$ éléments du dataset appartiennent à $\mathcal{X}\in \mathbb{R}^p$ ($p$ features numériques), et les labels appartiennent à $\mathcal{Y}\in \{+1, -1\}$.
On définit en outre $\mathcal{G}$ notre ensemble de classifieurs faibles : on peut par exemple prendre $\mathcal{G}$ ensemble des arbres de décision de profondeur $2$, etc...
On définit l'ensemble suivant:
$$\mathcal{H}_{\lambda} = \{ \lambda \times \text{Conv}(\mathcal{G}) \} $$$$=\left \{ \sum_{j=1}^{m} \alpha_j g_j, m \in \mathbb{N}^*, g_1, ..., g_m \in \mathcal{G}, \alpha_1, ..., \alpha_m \geq 0, \sum_{i=1}^m \alpha_j \leq \lambda \right\} $$
Dans ce contexte, le predicteur de boosting $\hat{h}^{\text{boost}}$ est le suivant :
$$\hat{h}^{\text{boost}} \in \arg \underset{h\in \mathcal{H}_{\lambda}}{\min} \frac{1}{n} \sum_{i=1}^n e^{-y_i h(x_i)}$$Les raisons pour lesquelles on pose le prédicteur de boosting de cette manière seront détaillées dans les cours plus théoriques de machine learning au S2 de 2A et en 3A. L'important est de remarquer qu'on est placé dans un problème d'optimisation convexe, qui peut donc nous garantir un minimum global de la fonction de perte; mais qui a le désavantage d'être défini sur un espace de dimension infinie.
On ne va en général pas pouvoir résoudre de manière exacte ce programme d'optimisation. À la place, on approxime la solution avec un algorithme. Ici, on présente adaboost:
- Initialisation : On définit $n$ poids initiaux (pour chacuns des exemples). Pour ne mettre aucun à priori, on donne à chaque élément le même poids : $$w_1 = \frac{1}{n}, ..., w_n = \frac{1}{n} $$
- Pour m = 1, ..., M :
- a) On définit le prédicteur $\hat{g}_m \in \arg \underset{g\in \mathcal{G}}{\min} \sum_{i=1}^n w_i 1_{\{y_i \neq g(x_i)\}}$, c'est à dire que l'on minimise la fonction de perte pondérée sur notre ensemble de prédicteurs $\mathcal{G}$. En suivant l'exemple des arbres de profondeur 2, $\hat{g}_m$ est l'arbre de profondeur 2 qui minimise la perte pondérée par les poids au temps $m$.
- b) Ensuite, on évalue l'erreur commise au temps $m$ : $\hat{e}_m = \frac{1}{\sum_{i=1}^n w_i} \left ( \sum_{i=1}^n w_i 1_{\{y_i \neq \hat{g}(x_i)\}}\right)$.
- c) On définit un poids qui sera associé au prédicteur $\hat{g}_m$ : $\hat{a}_m = \frac{1}{2}\ln \left( \frac{1-\hat{e}_m}{\hat{e}_m}\right)$
- d) Enfin, on met à jour les poids associés à chacuns des exemples du dataset pour mettre l'emphase sur les éléments mal classifiés, afin qu'à chaque étape l'algorithme apprenne à classifier les éléments les plus difficiles à classifier : $$ w_i := w_i \exp \{ \hat{a}_m 1_{\{y_i \neq \hat{g}_m\}}\}$$
- Pour m = 1, ..., M :
- Enfin, on définit le classifieur final : $$\hat{g}(x) = \text{signe} \left( \sum_{m=1}^M \hat{a}_m \hat{g}_m \right)$$
Cet algorithme d'approximation a historiquement fourni de bons résultats sur les benchmarks classiques en machine learning. Voyons ce qu'il donne, en l'évaluant évidemment avec le cross val score.
6. Evaluer les performances d'adaboost, avec le cross-val score, sur iris et diabetes¶
from sklearn.ensemble import AdaBoostClassifier
clf = AdaBoostClassifier() #DecisionTreeClassifier(max_depth=1)
print(cross_val_score(clf, X_iris, y_iris, cv=10).mean())
On est à pein meilleur que sur les algorithmes précédents. Regardons l'évolution en fonction du base estimator (l'ensemble $\mathcal{G}$ de prédicteurs faibles) choisi.
7. Les performances changent-elles en fonction du base_estimator choisi?¶
Estimons les performances d'Adaboost pour des ensembles $\mathcal{G}$ d'arbres de profondeur croissante :
from sklearn.tree import DecisionTreeClassifier
perf_dt = []
for i in range(1,20):
clf = AdaBoostClassifier(DecisionTreeClassifier(max_depth=i))
perf_dt.append(cross_val_score(clf, X_iris, y_iris, cv=10).mean())
plt.figure()
plt.plot(np.arange(1, 20), perf_dt)
plt.xlabel('Profondeur des arbres')
plt.ylabel('Crossval score')
plt.show()
Regardons les performances de ce qui semble être empiriquement la profondeur optimale, pour ce dataset:
clf = AdaBoostClassifier(DecisionTreeClassifier(max_depth=3))
print(cross_val_score(clf, X_iris, y_iris, cv=10).mean())
Passons maintenant au problème de régression, avec :
- Des arbres de régression
- Des régressions linéraires
from sklearn.ensemble import AdaBoostRegressor
rg = AdaBoostRegressor() #DecisionTreeClassifier(max_depth=1) - default parameter
print(cross_val_score(rg, X_diabetes, y_diabetes, cv=10, scoring='r2').mean())
from sklearn.tree import DecisionTreeRegressor
perf_dt = []
for i in range(1,20):
rg = AdaBoostRegressor(DecisionTreeRegressor(max_depth=i))
perf_dt.append(cross_val_score(rg, X_diabetes, y_diabetes, cv=10, scoring='r2').mean())
plt.figure()
plt.plot(np.arange(1, 20), perf_dt)
plt.xlabel('Profondeur des arbres')
plt.ylabel('Crossval score - R2')
plt.show()
from sklearn.linear_model import LinearRegression
rg = AdaBoostRegressor(LinearRegression())
print(cross_val_score(rg, X_diabetes, y_diabetes, cv=10, scoring='r2').mean())
Hyperparameter tuning¶
On voit qu'une fois que l'on sait entraîner un algorithme de machine learning, et l'évaluer de manière robuste grâce au cross-val score, il reste une dernière chose sur laquelle on peut jouer : les hyperparamètres.
Il existe de nombreux algorithmes d'optimisation d'hyperparamètres mais on n'utilisere ici que GridSearch, qui est assez simple. Comme son nom l'indique, on fournit à l'algorithme d'optimisation une grille d'hyperparamètres à tester. L'algorithme teste toutes les combinaisons d'hyperparamètres et retourne le meilleur modèle.
8. Avec GridSearch, trouver les meilleurs paramètres en boosting sur la problème de régression de diabetes¶
AdaBoostRegressor().get_params()
from sklearn.linear_model import Lasso
param_grid = {'base_estimator': [DecisionTreeClassifier(), LinearRegression(), Lasso()],
'learning_rate': [0.001, 0.01, 0.1, 1.0],
'n_estimators': [50, 100, 150, 200]}
from sklearn.model_selection import GridSearchCV
rg = GridSearchCV(AdaBoostRegressor(), param_grid=param_grid, scoring='r2')
rg.fit(X_diabetes, y_diabetes)
import pandas as pd
df = pd.DataFrame(rg.cv_results_)
df.sort_values('rank_test_score')