Introduction à l'apprentissage automatique: TP2 - Exercice 3


Mélanges de gaussiennes et jeu de données "Old Faithful", prédiction par Gaussian Mixture Regression


On va étudier un jeu de données relatif au geyser "Old Faithful" dans le parc du Yellowstone aux Etats-Unis.

Chaque observation est constituée de deux caractéristiques: la durée d'une éruption et la durée avant l'éruption suivante.

Les données sont décrites ici: http://www.stat.cmu.edu/~larry/all-of-statistics/=data/faithful.dat
et doivent être téléchargées sur Arche ou sur la page web du cours.


On commence par charger les bibliothèques qui nous seront utiles:

In [ ]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from sklearn import mixture

%matplotlib inline 

puis on charge les données du fichier old_faithful.txt:

In [ ]:
# attention: le notebook doit être ouvert dans le même répertoire que le fichier old_faithful.txt
data=np.loadtxt('old_faithful.txt')
print(data)  # pour vérifier que les données sont bien chargées

On commence par représenter les données par un nuage de point.

In [ ]:
plt.figure(figsize=[12,10]);
plt.scatter(data[:,1], data[:,2]);
plt.grid()
plt.xlabel('durée éruption')
plt.ylabel('temps avant la prochaine éruption')
plt.title('Old Faithful');

Modélisation d'une densité de probabilité par mélange de gaussiennes

Les observations sont des couples $(x,y)$ où $x$ est la durée d'éruption et $y$ le temps avant la prochaine éruption. Nous allons modéliser la loi jointe du couple $(x,y)$ à l'aide d'un mélange de gaussiennes.

Question 1. Observez l'évolution du critère BIC dans l'ajustement d'un mélange de gaussiennes à $M$ composantes, pour $M$ variant entre 1 et 10 (inspirez-vous du code ci-dessous, et utilisez la documentation). Vous afficherez les lignes d'iso-valeur du logarithme de la densité de probabilité estimée pour chaque valeur de $M$ en utilisant la fonction show_mixture ci-dessous (inspirée du code de la documentation scikit-learn).

Constatez que l'évolution du critère BIC semble justifier que ce geyser présente deux "modes de fonctionnement".


Remarque 1 : $M$ est appelé un hyperparamètre du modèle car il doit être fixé par l'utilisateur, à l'aide d'indicateurs comme BIC ici. Les paramètres du modèle sont les poids, moyennes, et matrices de covariance: ils sont estimés à partir des données par la fonction fit ci-dessous (qui, ici, met en oeuvre l'algorithme EM).

Remarque 2 : par définition, la log-vraisemblance calculée à partir d'une observation (ce qui est retourné par score_samples ci-dessous) est le logarithme de la valeur de la densité de probabilité en cette observation.

In [ ]:
def show_mixture(mixture,X_train):
    # mixture: objet GaussianMixture ajusté à des données
    # X_train: pour le scatter-plot, observations sous forme d'un tableau (observations en lignes, caractéristiques en colonnes)
    x = np.linspace(1, 5.5)
    y = np.linspace(20., 120.)
    X, Y = np.meshgrid(x, y)
    XX = np.array([X.ravel(), Y.ravel()]).T
    Z = -mixture.score_samples(XX)  # score_samples est le log de la densité (log-vraisemblance) en chaque point  
    Z = Z.reshape(X.shape)
    plt.figure(figsize=[8,6])
    CS = plt.contour(X, Y, Z, norm=LogNorm(vmin=1.0, vmax=100.0),
                    levels=np.logspace(0, 3, 10))   # lignes d'iso-valeur
    CB = plt.colorbar(CS, shrink=0.8, extend='both')
    plt.scatter(X_train[:, 0], X_train[:, 1])  # nuage des observations
    plt.scatter(mixture.means_[:,0],mixture.means_[:,1],c='red')  # moyennes des composantes en rouge
    plt.title('Level lines of negative log-likelihood predicted by a GMM, K='+str(mixture.n_components))
    plt.axis('tight')
    plt.show();
In [ ]:
# votre code ici

# exemple de code permettant de déterminer les paramètres d'un mélange à 4 composantes gaussiennes:
GMM = mixture.GaussianMixture(n_components=4)
GMM.fit(data[:,1:3])   # estimation des paramètres du mélange  (dans l'écriture Python 1:3, 3 est exclu)
show_mixture(GMM,data[:,1:3])   # pour tracer les courbes de niveau de ce mélange
print(GMM.bic(data[:,1:3]))  # affiche la valeur de BIC sur le jeu de données d'entraînement

Question 2. Ecrivez les paramètres (poids, moyennes, matrices de covariance) du mélange de gaussienne optimal au sens de BIC.

In [ ]:
# votre code ici (affichez les paramètres du "meilleur" mélange)

Question 3. La "méthode du coude" de $K$-Means permet-elle de justifier un classification en deux groupes ?

In [ ]:
# votre code ici
# voir le code de l'exercice 1 

Illustration de la théorie développée au chapitre 4, section 4.1.1: peut être vu comme un complément facultatif cette année

Prédiction du temps avant la prochaine éruption par Gaussian mixture regression

Nous allons à présent prédire le temps d'attente avant la prochaine éruption en fonction de la durée d'une éruption.

Si $x$ désigne la durée d'une éruption et $y$ le temps avant la prochaine éruption, nous venons de modéliser la densité $p(x,y)$ de la loi de probabilité jointe de $x$ et $y$ comme un mélange de deux gaussiennes.

Nous verrons lors de la prochaine séance (polycopié chapitre 4, section 4.1.1) que la meilleure prédiction (en un sens précis) de $y$ en fonction de $x$ est donnée par l'espérance conditionnelle suivante: $$ h(x) = E(y|x) = \int_y y p(y|x) dy$$ Il s'agit de l'espérance de $y$ connaissant $x$. Autrement dit, il s'agit de la valeur moyenne de $y$ à $x$ fixé.

Par définition, $$p(y|x) = \frac{p(x,y)}{p(x)}$$ et la densité marginale $p(x)$ est déterminée par: $$ p(x) = \int_y p(x,y) dy$$ Malheureusement, GaussianMixture de scikit-learn n'implémente pas ces calculs. Nous allons donc déterminer l'espérance conditionnelle $h(x)$ par intégration numérique, à partir du mélange de gaussiennes $p(x,y)$.

C'est ce que fait la fonction GMM_regression dans la cellule suivante: ses arguments sont 1) un objet GaussianMixture représentant un mélange de Gaussiennes, et 2) un tableau de données $x$; elle retourne le tableau des valeurs prédites pour les valeurs de $x$ fournies.

In [ ]:
def GMM_regression(GMM,x):
    y = np.linspace(20., 120.,num=1001)  # y discrétisé entre 20 et 120 par pas de 0.1 (1001 valeurs)
    X, Y = np.meshgrid(x, y)
    XX = np.array([X.ravel(), Y.ravel()]).T
    p = np.exp(GMM.score_samples(XX))   # score_samples fournit la log-vraisemblance d'une observation, p est donc bien la valeur de la densité  
    p = p.reshape(X.shape)  # on construit ainsi la carte des p(x,y)
    condexp = np.sum(Y*p,axis=0)/np.sum(p,axis=0)  # espérance conditionnelle obtenue par intégration numérique (ici, le première axe correspond à y)
    return condexp

Question 4. Complétez le code suivant de manière à tracer le graphique de la prédiction du temps entre éruptions en fonction de la durée d'éruption.

Comparez votre prédiction avec ce qui est proposé sur cette page du National Park Service. (nos données sont un peu ancienne, le "fonctionnement" du geyser évolue au fil des années)

In [ ]:
GMM=mixture.GaussianMixture(n_components=2,covariance_type='full')
GMM.fit(data[:,1:3])
x = np.linspace(1, 5.5)

# votre code ici

Question 5. Tracez la droite de régression linéaire superposée à l'estimation par espérance conditionnelle dans laquelle la loi jointe $p(x,y)$ est supposée gaussienne (il suffit de faire n_components=1 dans GaussianMixture), et constatez que les prédictions sont identiques. La raison est expliquée au chapitre 4 du polycopié, à lire pour la prochaine séance.

In [ ]:
# votre code ici