Rappels sur scikit-learn et le machine learning (correction)

Links: notebook, html, PDF, python, slides, GitHub

Quelques exercices simples sur scikit-learn. Le notebook est long pour ceux qui débutent en machine learning et sans doute sans suspens pour ceux qui en ont déjà fait.

from jyquickhelper import add_notebook_menu
add_notebook_menu()
%matplotlib inline

Des données synthétiques

On simule un jeu de données aléatoires.

from numpy import random
n = 1000
X = random.rand(n, 2)
X[:5]
array([[0.33151303, 0.70686719],
       [0.13039027, 0.58941167],
       [0.612744  , 0.37799233],
       [0.20215973, 0.11095186],
       [0.56857961, 0.10783821]])
y = X[:, 0] * 3 - 2 * X[:, 1] ** 2 + random.rand(n)
y[:5]
array([0.10121972, 0.49414321, 2.19975264, 0.74372472, 2.27103021])

Exercice 1 : diviser en base d’apprentissage et de test

Simple train_test_split.

from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y)

Exercice 2 : caler une régression linéaire

Et calculer le coefficient R^2. Pour ceux qui ne savent pas se servir d’un moteur de recherche : LinearRegression, r2_score.

from sklearn.linear_model import LinearRegression
reg = LinearRegression()
reg.fit(X_train, y_train)
LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)
from sklearn.metrics import r2_score
score = r2_score(y_test, reg.predict(X_test))
score
0.908988490753245

Exercice 3 : améliorer le modèle en appliquant une transformation bien choisie

Le modèle de départ est : Y = 3 X_1 - 2 X_2^2 + \epsilon. Il suffit de rajouter des featues polynômiales avec PolynomialFeatures.

from sklearn.preprocessing import PolynomialFeatures
poly = PolynomialFeatures()
poly.fit(X_train)
X_train2 = poly.transform(X_train)
reg2 = LinearRegression()
reg2.fit(X_train2, y_train)
LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)
score2 = r2_score(y_test, reg2.predict(poly.transform(X_test)))
score2
0.9362394073926681

Le coefficient R^2 est plus élevé car on utilise les mêmes variables que le modèle. Il n’est théoriquement pas possible d’aller au delà.

Exercice 4 : caler une forêt aléatoire

from sklearn.ensemble import RandomForestRegressor
rf = RandomForestRegressor()
rf.fit(X_train, y_train)
RandomForestRegressor(bootstrap=True, criterion='mse', max_depth=None,
           max_features='auto', max_leaf_nodes=None,
           min_impurity_decrease=0.0, min_impurity_split=None,
           min_samples_leaf=1, min_samples_split=2,
           min_weight_fraction_leaf=0.0, n_estimators=10, n_jobs=1,
           oob_score=False, random_state=None, verbose=0, warm_start=False)
r2_score(y_test, rf.predict(X_test))
0.9153506166386053

Le modèle linéaire est le meilleur modèle dans notre cas puisque les données ont été construites de la sorte. Il est attendu que le R^2 ne soit pas plus élevé tout du moins pas significativement plus élevé. On regarde avec les features polynômiales…

rf2 = RandomForestRegressor()
rf2.fit(X_train2, y_train)
r2_score(y_test, rf2.predict(poly.transform(X_test)))
0.9119002367619022

Avant de tirer des conclusions hâtives, il faudrait recommencer plusieurs fois l’expérience avant de dire que la performance est plus ou moins élevée avec ces features ce que ce notebook ne fera pas puisque la réponse théorique est connue dans ce cas.

Exercice 5 : un peu de math

Comparer les deux modèles sur les données suivantes ? Que remarquez-vous ? Expliquez pourquoi ?

X_test2 = random.rand(n, 2) + 0.5
y_test2 = X_test2[:, 0] * 3 - 2 * X_test2[:, 1] ** 2 + random.rand(n)
res = []
for model in [reg, reg2, rf, rf2]:
    name = model.__class__.__name__
    try:
        pred = model.predict(X_test)
        pred2 = model.predict(X_test2)
    except Exception:
        pred = model.predict(poly.transform(X_test))
        pred2 = model.predict(poly.transform(X_test2))
        name += " + X^2"
    res.append(dict(name=name, r2=r2_score(y_test, pred),
                    r2_jeu2=r2_score(y_test2, pred2)))

import pandas
df = pandas.DataFrame(res)
df
name r2 r2_jeu2
0 LinearRegression 0.908988 0.682467
1 LinearRegression + X^2 0.936239 0.948110
2 RandomForestRegressor 0.915351 0.493273
3 RandomForestRegressor + X^2 0.911900 0.517105

Le seul modèle qui s’en tire vraiment est la régression linéaire avec les features polynômiales. Comme il équivaut au modèle théorique, il est normal qu’il ne se plante pas trop même si ses coefficients ne sont pas identique au modèle théorique (il faudrait plus de données pour que cela converge).

reg2.coef_, reg2.intercept_
(array([ 0.        ,  2.81692538, -0.29768531,  0.08662761,  0.13367719,
        -1.7515442 ]), 0.5889925538787228)

Pour les autes modèles, voyons déjà visuellement ce qu’il se passe.

Exercice 6 : faire un graphe avec…

Je laisse le code décrire l’approche choisie pour illustrer les carences des modèles précédents. Le commentaire suit le graphique pour les paresseux.

import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 2, figsize=(14, 4))

a, b = 0.9, 1.1
index1 = (X_test2[:, 0] >= a) & (X_test2[:, 0] <= b)
index2 = (X_test2[:, 1] >= a) & (X_test2[:, 1] <= b)
yth = X_test2[:, 0] * 3 - 2 * X_test2[:, 1]


ax[0].set_xlabel("X1")
ax[0].set_ylabel("Y")
ax[0].plot(X_test2[index2, 0], yth[index2], '.', label='Y théorique')
ax[1].set_xlabel("X2")
ax[1].set_ylabel("Y")
ax[1].plot(X_test2[index1, 1], yth[index1], '.', label='Y théorique')

for model in [reg, reg2, rf, rf2]:
    name = model.__class__.__name__
    try:
        pred2 = model.predict(X_test2)
    except Exception:
        pred2 = model.predict(poly.transform(X_test2))
        name += " + X^2"
    ax[0].plot(X_test2[index2, 0], pred2[index2], '.', label=name)
    ax[1].plot(X_test2[index1, 1], pred2[index1], '.', label=name)
ax[0].legend()
ax[1].legend();
../_images/ml_scikit_learn_simple_correction_28_0.png

Le graphe étudie les variables des modèles selon une coordonnées tout en restreignant l’autre dans un intervalle donné. On voit tout de suite que la forêt aléatoire devient constante au delà d’un certain seuil. C’est encore une fois tout à fait normal puisque la base d’apprentissage ne contient des X_1, X_2 que dans l’intervalle [0, 1]. En dehors, chaque arbre de décision produit une valeur constante tout simplement parce que ce sont des fonctions en escalier : une forêt aléatoire est une moyenne de fonctions en escalier, elle est bornée. Quant à la première régression linéaire, elle ne peut saisir les effets du second degré, elle est linéaire par rapport aux variables de départ. Elle s’écarte moins mais elle s’écarte quand même de la variable à prédire.

Cet exercice a pour but d’illustrer qu’un modèle de machine learning est estimé sur un jeu de données qui suit une certaine distribution. Lorsque les données sur lesquelles le modèle est utilisé pour prédire ne suivent plus cette loi, les modèles retournent des réponses qui ont toutes les chances d’être fausses et ce, de manière différente selon les modèles.

C’est pour cela qu’on dit qu’il faut réapprendre régulièrement les modèles de machine learning, surtout s’ils sont appliqués sur des données générées par l’activité humaine et non des données issues de problèmes physiques.

Exercice 7 : illuster l’overfitting avec un arbre de décision

Sur le premier jeu de données.

from sklearn.tree import DecisionTreeRegressor

res = []
for md in range(1, 20):
    tree = DecisionTreeRegressor(max_depth=md)
    tree.fit(X_train, y_train)
    r2_train = r2_score(y_train, tree.predict(X_train))
    r2_test = r2_score(y_test, tree.predict(X_test))
    res.append(dict(profondeur=md, r2_train=r2_train, r2_test=r2_test))

df = pandas.DataFrame(res)
df.head()
profondeur r2_test r2_train
0 1 0.380103 0.446387
1 2 0.636346 0.672284
2 3 0.786778 0.826594
3 4 0.872799 0.892911
4 5 0.896276 0.931297
ax = df.plot(x='profondeur', y=['r2_train', 'r2_test'])
ax.set_title("Evolution du R2 selon la profondeur");
../_images/ml_scikit_learn_simple_correction_32_0.png

Exercice 8 : augmenter le nombre de features et régulariser une régression logistique

L’objectif est de regarder l’impact de la régularisation des coefficients d’une régression logistique lorsque le nombre de features augmentent. On utilise les features polynômiales et une régression Ridge ou Lasso.

from sklearn.linear_model import Ridge, Lasso
import numpy.linalg as nplin
import numpy

def coef_non_nuls(coef):
    return sum(numpy.abs(coef) > 0.001)

res = []
for d in range(1, 21):
    poly = PolynomialFeatures(degree=d)
    poly.fit(X_train)
    X_test2 = poly.transform(X_test)

    reg = LinearRegression()
    reg.fit(poly.transform(X_train), y_train)
    r2_reg = r2_score(y_test, reg.predict(X_test2))

    rid = Ridge(alpha=10)
    rid.fit(poly.transform(X_train), y_train)
    r2_rid = r2_score(y_test, rid.predict(X_test2))

    las = Lasso(alpha=0.01)
    las.fit(poly.transform(X_train), y_train)
    r2_las = r2_score(y_test, las.predict(X_test2))

    res.append(dict(degre=d, nb_features=X_test2.shape[1],
                    r2_reg=r2_reg, r2_las=r2_las, r2_rid=r2_rid,
                    norm_reg=nplin.norm(reg.coef_),
                    norm_rid=nplin.norm(rid.coef_),
                    norm_las=nplin.norm(las.coef_),
                    nnul_reg=coef_non_nuls(reg.coef_),
                    nnul_rid=coef_non_nuls(rid.coef_),
                    nnul_las=coef_non_nuls(las.coef_),
                    ))

df = pandas.DataFrame(res)
df
degre nb_features nnul_las nnul_reg nnul_rid norm_las norm_reg norm_rid r2_las r2_reg r2_rid
0 1 3 2 2 2 3.403924 3.568168e+00 3.079896 0.906854 0.908988 0.892232
1 2 6 4 5 5 3.280738 3.334211e+00 2.350989 0.934247 0.936239 0.919758
2 3 10 4 9 9 3.280888 5.090698e+00 2.072340 0.934247 0.935657 0.921656
3 4 15 4 14 14 3.280963 1.180973e+01 1.992448 0.934248 0.934934 0.921532
4 5 21 4 20 20 3.281022 6.657031e+01 1.972762 0.934248 0.935177 0.921656
5 6 28 4 27 26 3.281037 3.242570e+02 1.965189 0.934248 0.934388 0.921823
6 7 36 4 35 35 3.281040 1.428490e+03 1.958965 0.934248 0.933246 0.921927
7 8 45 4 44 43 3.281041 1.441240e+04 1.952981 0.934248 0.931307 0.921976
8 9 55 4 55 53 3.281041 4.893616e+13 1.947737 0.934248 0.616562 0.921997
9 10 66 4 65 65 3.281041 4.659569e+05 1.943600 0.934248 0.929751 0.922013
10 11 78 4 77 75 3.281041 2.360073e+06 1.940623 0.934248 0.925344 0.922035
11 12 91 4 91 88 3.281041 1.708086e+07 1.938666 0.934248 0.921123 0.922067
12 13 105 4 105 99 3.281041 9.955733e+07 1.937518 0.934248 0.874136 0.922108
13 14 120 4 120 111 3.281041 6.854497e+08 1.936961 0.934248 0.927716 0.922158
14 15 136 4 136 127 3.281041 3.786997e+09 1.936808 0.934248 0.842488 0.922212
15 16 153 4 153 145 3.281041 4.467998e+10 1.936913 0.934248 0.664588 0.922269
16 17 171 4 171 161 3.281041 2.361809e+11 1.937165 0.934248 -0.726442 0.922325
17 18 190 4 190 179 3.281041 1.599035e+12 1.937489 0.934248 0.582385 0.922379
18 19 210 4 210 203 3.281041 4.455355e+13 1.937834 0.934248 -25.406536 0.922429
19 20 231 4 231 223 3.281041 2.262080e+13 1.938168 0.934248 -21.684447 0.922475
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
df.plot(x="nb_features", y=["r2_reg", "r2_las", "r2_rid"], ax=ax[0])
ax[0].set_xlabel("Nombre de features")
ax[0].set_ylim([0, 1])
ax[0].set_title("r2")
df.plot(x="nb_features", y=["nnul_reg", "nnul_las", "nnul_rid"], ax=ax[1])
ax[1].set_xlabel("Nombre de features")
ax[1].set_title("Nombre de coefficients non nuls");
../_images/ml_scikit_learn_simple_correction_35_0.png

Numériquement, la régression linéaire devient difficile à estimer lorsque le nombre de features augmente. Théoriquement, il ne devrait pas y avoir de baisse de performances mais le graphe montre des erreurs évidentes. Cela se traduit par une norme des coefficients qui explose. La régularisation parvient à contraindre les modèles. La régression Ridge produira beaucoup de petits coefficients non nuls, la régression Lasso préfèrera concentrer la norme sur quelques coefficients seulement. Cette observation n’est vraie que dans le cas d’une régression linéaire avec une erreur quadratique.