Travaux pratiques - Introduction à l’apprentissage supervisé

_images/jupyter_logo.png

Cahier Jupyter

L’objectif de cette première séance de travaux pratiques est de vous faire découvrir les problèmes de classement et les problèmes de régression avec des modèles linéaires et des perceptrons multi-couches (PMC).

Références externes utiles :

Note

Il est fortement conseillé d’utiliser Python 3 pour ces travaux pratiques, voir cette page d’installation. Python 2 n’est plus maintenu depuis janvier 2020 et scikit-learn a abandonné son support depuis sklearn 0.20. Certaines fonctionnalités plus récentes pourraient ne pas être disponibles sous Python 2.x.

# Import des bibliothèques utiles
import matplotlib.pyplot as plt
import numpy as np
import sklearn

Régression avec modèles linéaires et polynomiaux

Pour commencer ce premier TP, nous examinerons un problème simple de régression avec une seule variable prédictive, similaire à celui illustré dans le support de cours. À partir d’un ensembles d’observations \(\{(x_1, y_1), (x_2, y_2), \dots, (x_n, y_n)\}\), l’objectif est de trouver la fonction \(\hat{f}\) telle que \(\hat{f}(x_i) = y_i\).

Pour cet exemple, commençons par générer des données synthétiques. En pratique, nous n’aurions bien entendu pas à réaliser ce travail : le jeu de données est constitué des observations réelles. Notre jeu de données d’exemple consiste en une portion de sinusoïde sur laquelle un bruit gaussien a été appliqué.

# Fixe l'aléas pour pouvoir reproduire les résultats
np.random.seed(0)
# On tire au hasard 100 points dans l'intervalle [0,2]
X = 2 * np.random.rand(100, 1)
# On calcule la valeur de sin(x) pour chaque point, plus un bruit gaussien aléatoire
y = np.sin(X[:,0]) + 0.15 * np.random.randn(100)

Nous pouvons visualisation toutes les données générées, qui forment notre matrice d’observation pour ce problème de régression :

plt.figure(figsize=(8, 6))
plt.scatter(X, y, s=50)
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.show()

Nous générons maintenant un premier découpage entre données d’apprentissage et données de test. Ce découpage est aléatoire et s’effectue à l’aide de la fonction train_test_split de scikit-learn. Le paramètre test_size permet de spécifier le pourcentage du jeu de données qui sera contenu dans le jeu de test.

from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33, random_state=0)
plt.figure(figsize=(8, 6))
plt.scatter(X_train, y_train, s=50, edgecolors='none', label="Exemples d'entraînement")
plt.scatter(X_test, y_test, c='none', s=50, edgecolors='blue', label="Exemples d'évaluation")
plt.legend()
plt.show()

Nous cherchons d’abord un modèle linéaire pour ce problème de régression. Regardez la classe LinearRegression de scikit-learn (ainsi que ces explications). Dans le cas présent, scikit-learn va chercher la fonction affine \(f: x \rightarrow y = ax + b\) en déterminant les paramètres \(a\) et \(b\) à l’aide de la méthode des moindres carrés.

Les résultats sont évalués ici à travers le coefficient de détermination, qui est le rapport entre la variance expliquée par le modèle et la variance totale (de la variable expliquée), voir les explications.

Commençons par instancier un modèle de régression linéaire.

from sklearn import linear_model
reg = linear_model.LinearRegression()

Tous les modèles implémentés dans scikit-learn suivent la même interface et présentent au moins les trois méthodes suivantes :

  • .fit() permet d’apprendre le modèle en estimant ses paramètres à partir d’un jeu de données,

  • .predict() ou .transform() permet d’appliquer le modèle à de nouvelles données,

  • .score() permet d’évaluer le modèle selon un critère prédéfini (taux de bonne classification, coefficient de détermination, etc.) sur un jeu de test.

Par exemple, dans notre cas, la méthode .fit(X,y) permet de déterminer les paramètres de la régression linéaire à partir de nos observations X et de notre vérité terrain y.

# optimisation du modèle: détermination des paramètres de la régression linéaire par la méthode des moindres carrés
reg.fit(X_train, y_train)

Nous pouvons évaluer à quel point le modèle « colle » aux données à l’aide du coefficient de détermination calculé sur les exemples du jeu d’apprentissage :

# attention, score() ici ne renvoie pas l'erreur mais la valeur du coefficient de détermination R² !
coeff_train = reg.score(X_train, y_train)
print(f"Coefficient de détermination R² en train : {coeff_train:.2f}")

Question:

Calculer le coefficient de détermination R² sur le jeu de test. Que constatez-vous ?

Correction

coeff_test = reg.score(X_test, y_test)
print(f"Coefficient de détermination R² en test : {coeff_test:.2f}")

Nous pouvons ensuite, à partir des coefficients qui ont été estimés, tracer notre modèle de régression :

plt.figure(figsize=(8, 6))
plt.scatter(X_train,y_train, s=50, edgecolors='none', label="Exemples d'apprentissage")
plt.scatter(X_test,y_test, c='none', s=50, edgecolors='blue', label="Exemples d'évaluation")
plt.xlabel("$x$")
plt.ylabel("$y$")

x_min, x_max = plt.xlim()
nx = 100
xx = np.linspace(x_min, x_max, nx).reshape(-1,1)
plt.plot(xx,reg.predict(xx), color='k', label="Régression linéaire")
plt.title("Régression linéaire unidimensionnelle")
plt.legend()
plt.show()

Question

Calculez l’erreur quadratique moyenne du modèle sur les données d’apprentissage et ensuite sur les données de test. Vous pouvez l’implémenter vous-même à l’aide de NumPy ou bien consulter la documentation du module sklearn.metrics.

Correction

from sklearn.metrics import mean_squared_error
y_pred_train = reg.predict(X_train)
y_pred_test = reg.predict(X_test)

mse_train = mean_squared_error(y_train, y_pred_train)
mse_test = mean_squared_error(y_test, y_pred_test)
print(f"MSE = {mse_train:.3f} (train)")
print(f"MSE = {mse_test:.3f} (test)")

Régression polynomiale

La régression linéaire est un outil pratique mais quelque peu limité : les relations entre les variables explicatives et la variable à prédire sont rarement linéaires en pratique ! Nous pourrions chercher dans une famille paramétrique plus grande, par exemple les polynômes à une seule variable de degré \(d\) : \(P : x \rightarrow a_0 + a_1 x + a_2 x^2 + \dots + a_d x^d\).

def polynomial_fit(degree, show=False):
      pipe = Pipeline([('poly', PolynomialFeatures(degree=degree)),
                      ('reg', linear_model.LinearRegression())])

      pipe.fit(X_train, y_train)
      train_score = pipe.score(X_train, y_train)
      test_score = pipe.score(X_test, y_test)

      if show:
          plt.figure(figsize=(8, 6))
          plt.scatter(X_train, y_train, s=50, edgecolors='none', label="Exemples d'apprentissage")
          plt.scatter(X_test, y_test, c='none', s=50, edgecolors='blue', label="Exemples de test")
          x_min, x_max = plt.xlim()
          xx = np.linspace(x_min, x_max, nx).reshape(-1,1)
          plt.plot(xx,pipe.predict(xx),color='black', label="Régression polynomiale")
          plt.legend()
          plt.show()

          print(f"Coefficient R² (train): {train_score:.3f}")
          print(f"Coefficient R² (test): {test_score:.3f}")
      return train_score, test_score

polynomial_fit(2, show=True)

Question

Que se passe-t-il lorsque l’on modifie le degré du polynôme ? Quelles observations pouvez-vous faire concernant le sous-apprentissage ? Le sur-apprentissage ?

Correction

Lorsque le degré du polynôme est faible (0 ou 1), la correspondance entre la courbe réelle et le modèle de régression est mauvaise (score R² faible, voire négatif). Nous sommes en sous-apprentissage: l’erreur d’approximation est élevée car le modèle que l’on choisit (constante ou régression linéaire) est une mauvaise hypothèse par rapport aux données réelles.

Lorsque l’on augmente le degré, le coefficient R² se rapproche de 1 sur les exemples d’apprentissage, mais diminue sur les exemples de test. Visuellement, la variance des prédictions augmente: nous sommes en sur-apprentissage. Le modèle dispose de trop de paramètres libres qui viennent « coller » au bruit statistique dans les données, mais cela ne généralise pas sur les exemples de test.

Question

(pour aller plus loin) Tracer les courbes du coefficient R² sur le jeu d’apprentissage et sur le jeu de test en fonction du degré du polynôme (entre 1 et 10). Vous pouvez utiliser la fonction polynomial_fit pour ce faire qui renvoie un tuple contenant les deux scores. Que constatez-vous ?

Correction

degrees = np.arange(1, 11, 1)
train_scores, test_scores = zip(*[polynomial_fit(d, show=False) for d in degrees])

fig = plt.figure(figsize=(8, 6))
plt.plot(degrees, train_scores, label="Train")
plt.plot(degrees, test_scores, label="Test")
plt.xlabel("Degré du polynôme")
plt.ylabel("R²")
plt.xlim(1, 10)
plt.show()

Classement avec modèles linéaires et perceptrons multicouches

Voyons maintenant comment les notions de sur-apprentissage et sous-apprentissage peuvent se manifester dans le cas d’un problème de classement (ou classification). Nous examinerons un problème simple de discrimination entre deux classes, similaire à celui illustré dans le support de cours.

Nous générons des données à partir de lois normales bidimensionnelles. Pour la première classe nous employons une seule loi avec des variances différentes et une rotation qui introduit des covariances. La seconde classe est un mélange de 4 lois normales avec des centres différents.

import numpy as np
import matplotlib.pyplot as plt
# fixer la graine aléatoire de numpy permet d'obtenir systématiquement les mêmes résultats
np.random.seed(150)

from sklearn import datasets

dataset = "gaussian"

if dataset == "gaussian":
    # définir matrices de rotation et de dilatation
    rot = np.array([[0.94, -0.34], [0.34, 0.94]])
    sca = np.array([[3.4, 0], [0, 2]])
    # générer données classe 1
    c1d = (np.random.randn(100,2)).dot(sca).dot(rot)

    # générer données classe 2
    c2d1 = np.random.randn(25,2)+[-10, 2]
    c2d2 = np.random.randn(25,2)+[-7, -2]
    c2d3 = np.random.randn(25,2)+[-2, -6]
    c2d4 = np.random.randn(25,2)+[5, -7]

    data = np.concatenate((c1d, c2d1, c2d2, c2d3, c2d4))

    # générer étiquettes de classe
    l1c = np.ones(100, dtype=int)
    l2c = np.zeros(100, dtype=int)
    labels = np.concatenate((l1c, l2c))
elif dataset == "iris":
    iris_X, iris_y = datasets.load_iris(return_X_y=True)
    data = iris_X[:, :2]
    labels = iris_y
elif dataset == "wine":
    wine_X, wine_y = datasets.load_wine(return_X_y=True)
    data = wine_X[:, (0,10)]
    labels = wine_y

Visualisation de toutes les données :

# les échantillons du premier groupe sont en rouge 'r', ceux du deuxième groupe en vert ("green") 'g'
cmp = np.array(['r','g', 'b'])
plt.figure(figsize=(8, 6))
plt.scatter(data[:,0],data[:,1], c=cmp[labels], s=50, edgecolors='none')
plt.show()

Nous générons maintenant un premier découpage entre données d’apprentissage et données de test. Les données de test sont affichées avec cercles vides (c='none'), les données d’apprentissage avec des cercles remplis.

from sklearn.model_selection import train_test_split

# découpage des données en train et test
X_train, X_test, y_train, y_test = train_test_split(data, labels, test_size=0.33, random_state=0)

plt.figure(figsize=(8, 6))
plt.scatter(X_train[:,0],X_train[:,1],c=cmp[y_train], s=50, edgecolors='none', label="Exemples d'apprentissage")
plt.scatter(X_test[:,0], X_test[:,1], c='none' ,s=50, edgecolors=cmp[y_test], label="Exemples de test")
plt.legend()
plt.show()

Question

Combien d’échantillons le jeu de données d’apprentissage contient-il ? La documentation de scikit-learn peut vous aider.

Correction

Le paramètre test_size=0.33 signifie que l’on conserve 33% des données pour le jeu de test. Le jeu d’apprentissage contient donc 67% des observations, soit 134 échantillons (le jeu de données gaussien contient 200 exemples).

Modèle linéaire : AFD

Nous cherchons d’abord un modèle linéaire pour ce problème de discrimination entre deux classes et utilisons pour cela l’étape décisionnelle de l’analyse factorielle discriminante (AFD). En résumé, cette méthode cherche la frontière de décision linéaire qui sépare au mieux les données en maximisant la variance entre les barycentres des deux groupes. Pour plus de détails, vous pouvez consulter la documentation de l’AFD dans scikit-learn ainsi que le chapitre et les travaux pratiques de l”UE RCP208) portant sur les méthodes factorielles.

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
lda = LinearDiscriminantAnalysis()

# évaluation et affichage sur split1
lda.fit(X_train, y_train)
print("Le score sur le jeu d'apprentissage est de : {:.3f}".format(lda.score(X_train, y_train)))

print("Le score sur le jeu de test est de : {:.3f}".format(lda.score(X_test, y_test)))

Nous pouvons examiner visuellement le modèle trouvé (la frontière de discrimination linéaire, c’est-à-dire ici une droite dans le plan) :

# on créé une nouvelle figure sur laquelle on affiche les points
plt.figure(figsize=(8, 6))
plt.scatter(X_train[:,0], X_train[:,1], c=cmp[y_train], s=50, edgecolors='none', label="Train")
plt.scatter(X_test[:,0],  X_test[:,1], c='none', s=50, edgecolors=cmp[y_test], label="Test")

# on calcule pour chaque point du plan sa probabilité d'appartenir à chaque classe
nx, ny = 400, 400
x_min, x_max = plt.xlim()
y_min, y_max = plt.ylim()
# meshgrid permet d'échantillonner tous les points du plan (entre x_min et x_max)
xx, yy = np.meshgrid(np.linspace(x_min, x_max, nx),np.linspace(y_min, y_max, ny))
# .predict_proba permet de prédire le score de la LDA pour un ensemble d'observations
Z = lda.predict_proba(np.c_[xx.ravel(), yy.ravel()])
for cls_idx in range(max(labels)):
    zz = Z[:, cls_idx].reshape(xx.shape)
    # on dessine la frontière correspond à un score de 0,5
    # les scores < 0,5 correspondent à la classe 0
    # les scores > 0,5 correspondent à la classe 1
    plt.contour(xx, yy, zz, [0.5])
plt.legend()
plt.show()

Modèle non-linéaire: perceptron multicouche

Comme vous avez pu le constater, la frontière de séparation idéale entre les deux groupes n’est pas linéaire. Il nous faut donc changer de famille paramétrique pour en choisir une de plus grande capacité. Nous allons donc utiliser un réseau de neurones simple (un perceptron multicouche) pour réaliser la discrimination binaire. Nous utilisons pour cela la classe MLPClassifier de scikit-learn (voir aussi ces explications concernant le fonctionnement général de cet objet). Il est très instructif d’examiner les valeurs par défaut des paramètres.

Les perceptrons multicouches ont été présentés dans le dernier chapitre de l’UE RCP208 portant sur les réseaux de neurones. Il n’est pas indispensable de connaître le fonctionnement précis des réseaux de neurones pour ce TP, ces modèles seront revus par la suite dans la partie du cours portant sur l’apprentissage profond.

Les trois paramètres importants que nous allons manipuler sont le coefficient d’oubli, le nombre de neurones par couche et le nombre de couches. Par défaut, scikit-learn instancie un réseau de neurones à une seule couche cachée contenant 100 neurones.

Nous utilisons d’abord un coefficient « d’oubli » (weight decay) alpha = 1e-5. Ce terme correspond à l’intensité de la régularisation L2 appliquée sur le perceptron. Pour rappel, l’objectif de la régularisation est de réduire la capacité du modèle.

from sklearn.neural_network import MLPClassifier
# nous utilisons ici l'algorithme L-BFGS pour optimiser le perceptron
clf = MLPClassifier(solver='lbfgs', alpha=1e-5, max_iter=1000, random_state=42)

# évaluation et affichage sur split1
clf.fit(X_train, y_train)
train_score = clf.score(X_train, y_train)
print(f"Le score en train est {train_score:.3f}")

test_score = clf.score(X_test, y_test)
print(f"Le score en test est {test_score:.3f}")

Question

À l’aide de la documentation de scikit-learn, déterminez à quoi correspond la valeur renvoyée par la méthode .score().

Correction

.score() calcule la métrique habituelle de classification: l”accuracy (en français, le taux de bonnes prédictions), c’est-à-dire le pourcentage d’exemples dont la classe a été correctement prédite par le modèle.

Comme précédemment pour l’analyse factorielle discriminante, nous pouvons désormais tracer la frontière de décision du perceptron multicouche que nous venons d’optimiser.

# créer une nouvelle figure
plt.figure()
# afficher les nuages de points apprentissage (remplis) et de test (vides)
plt.scatter(X_train[:,0],X_train[:,1],c=cmp[y_train],s=50, edgecolors='none', label="Train")
plt.scatter(X_test[:,0], X_test[:,1], c='none', s=50, edgecolors=cmp[y_test], label="Test")

# calculer la probabilité de sortie du perceptron pour tous les points du plan
nx, ny = 200, 200
x_min, x_max = plt.xlim()
y_min, y_max = plt.ylim()
xx, yy = np.meshgrid(np.linspace(x_min, x_max, nx),np.linspace(y_min, y_max, ny))
Z = clf.predict_proba(np.c_[xx.ravel(), yy.ravel()])
for cls_idx in range(max(labels)):
    zz = Z[:, cls_idx].reshape(xx.shape)
    # on dessine la frontière correspond à un score de 0,5
    # les scores < 0,5 correspondent à la classe 0
    # les scores > 0,5 correspondent à la classe 1
    plt.contour(xx, yy, zz, [0.5])
plt.legend()
plt.show()

Question

Refaites l’expérience avec alpha = 1. Dans quel cas la régularisation est plus forte ? Quelle est la conséquence sur les résultats ?

Correction

La régularisation est plus élevée lorsque alpha=1. La frontière est alors nettement moins irrégulière (on réduit la variance en diminuant la capacité du modèle).

La fonctions mlp_fit ci-dessous permet de manipuler simultanément les trois paramètres importants que nous avons décrit précédemment : coefficient d’oubli, nombre de couches, nombre de neurones par couche. Dans le perceptron, les paramètres à optimiser sont les poids des connexions entre les neurones. Plus il y a de couches et plus il y a de neurones dans une couche, plus le nombre de connexions est élevé (et par conséquent plus il y a de paramètres).

def mlp_fit(alpha, neurons, layers):
    # nous utilisons ici l'algorithme L-BFGS pour optimiser le perceptron
    layer_sizes = (neurons,) * layers
    clf = MLPClassifier(solver='lbfgs', alpha=alpha, hidden_layer_sizes=layer_sizes, max_iter=1000, tol=1e-3, random_state=0)

    # évaluation et affichage sur split1
    clf.fit(X_train, y_train)
    train_score = clf.score(X_train, y_train)
    print(f"Le score en train est {train_score:.3f}")

    test_score = clf.score(X_test, y_test)
    print(f"Le score en test est {test_score:.3f}")

    # on créé une nouvelle figure sur laquelle on affiche les points
    plt.figure(figsize=(8, 6))
    plt.scatter(X_train[:,0], X_train[:,1], c=cmp[y_train], s=50, edgecolors='none', label="Train")
    plt.scatter(X_test[:,0],  X_test[:,1], c='none', s=50, edgecolors=cmp[y_test], label="Test")

    # on calcule pour chaque point du plan sa probabilité d'appartenir à chaque classe
    nx, ny = 400, 400
    x_min, x_max = plt.xlim()
    y_min, y_max = plt.ylim()
    # meshgrid permet d'échantillonner tous les points du plan (entre x_min et x_max)
    xx, yy = np.meshgrid(np.linspace(x_min, x_max, nx),np.linspace(y_min, y_max, ny))
    # .predict_proba permet de prédire le score de chaque classe pour un ensemble d'observations
    Z = clf.predict_proba(np.c_[xx.ravel(), yy.ravel()])
    for cls_idx in range(max(labels)):
        zz = Z[:, cls_idx].reshape(xx.shape)
        # on dessine la frontière correspond à un score de 0,5
        # les scores < 0,5 correspondent à la classe 0
        # les scores > 0,5 correspondent à la classe 1
        plt.contour(xx, yy, zz, [0.5])
    plt.legend()
    plt.show()

mlp_fit(1e-5, 100, 1)

Question

Expérimentez avec la visualisation ci-dessus. Que se passe-t-il lorsque le nombre de couches augmente ? Lorsque le nombre de neurones augmente ? Comment pouvez-vous interpréter ces résultats du point de vue de la capacité du modèle, du sur-apprentissage et du sous-apprentissage ?

Correction

Augmenter le nombre de couches ou le nombre de neurones revient à augmenter le nombre de paramètres du modèle: la capacité augmente et donc le sur-apprentissage également.

Les résultats changent peu lorsque alpha est suffisamment élevé: la régularisation réduit artificiellement la capacité du modèle, même lorsque le nombre de paramètres augmente. C’est donc une façon de contrôler le sur-apprentissage, même pour un modèle de grande capacité.

Question

Générez d’autres découpages apprentissage/test en réexécutant la fonction train_test_split avec une autre valeur de graine aléatoire (random_state=). Examinez la variabilité des résultats lorsque ce découpage change. Dans quel cas elle est plus forte ?

Correction

La variabilité des résultats selon le découpage apprentissage/test est plus forte lorsque le modèle est en sur-apprentissage (grand nombre de couches/neurones et alpha faible) : le modèle mémorise du bruit statistique lié à l’échantillonnage des données. Les frontières sont irrégulières et modélisent moins bien la frontière réelle.

Pour aller plus loin

Reprenez la deuxième partie du TP sur le problème de discrimination à l’aide d’un modèle linéaire (AFD) et d’un modèle non-linéaire (PMC) mais en changeant le jeu de données. Vous pouvez modifier la variable dataset par "wine" ou "iris" pour utiliser respectivement le Wine Dateset ou les Iris de Fisher.

Le jeu de données Wine Dataset est une collection de 178 vins italiens répartis en 3 catégories, chacune correspondant à un vignoble différent. Treize mesures physico-chimiques sont disponibles mais dans ce TP nous n’utilisons que les deux premières (taux d’alcool et la couleur du vin). Plus de détails ici.

Les Iris de Fisher est un jeu de données célèbre en classification automatique. Il comporte 150 échantillons répartis en trois classes de fleurs (Iris Setosa, Iris Versicolour et Iris Virginica). Pour chaque plante, quatre caractéristiques végétales ont été mesurés : longueur sépale, largeur sépale, longueur de pétale, largeur de pétale. Dans ce TP nous n’utilisons que la longueur et la largeur sépale. Plus d’informations sur Iris ici.