Travaux pratiques - Estimation de densité

(correspond à 2 séances de TP)

Références externes utiles :

L’estimation de densités dans Scikit-learn

Pour estimer des densités, Scikit-learn propose les histogrammes et les noyaux comme méthodes non paramétriques, et les modèles de mélange comme méthodes paramétriques. Manquent, en revanche, des outils de choix du paramètre unique (largeur de « fenêtre ») pour les méthodes non paramétriques et des outils de comparaison entre densités (même lorsque celles-ci sont issues de mélanges de lois normales).

Les méthodes non paramétriques proposées sont présentées dans http://scikit-learn.org/stable/modules/density.html et l’utilisation de modèles de mélange dans http://scikit-learn.org/stable/modules/mixture.html.

Estimation par noyaux

Afin de permettre la visualisation des résultats et sachant que des outils de comparaison entre densités ne sont pas directement accessibles, nous examinerons uniquement des données unidimensionnelles et bidimensionnelles.

La description de l’implémentation de l’estimation de densités par noyaux se trouve dans http://scikit-learn.org/stable/modules/generated/sklearn.neighbors.KernelDensity.html.

Le type de noyau est spécifié dans le paramètre kernel est il est possible d’utiliser comme valeurs 'gaussian' (choix par défaut), 'tophat', 'epanechnikov', 'exponential', 'linear' ou 'cosine'. Une fois le type de noyau choisi, tous les noyaux (ou « fenêtres ») employés présentent une même valeur pour l’unique paramètre de largeur de « fenêtre » appelé ici bandwidth.

Il faut noter que le coût du calcul de la densité en un point donné peut être élevé si l’ensemble d’observations à partir duquel l’estimation de densité a été réalisée est très volumineux (en effet, cela impose de calculer la valeur du noyau entre le point donné et chaqune des observations de départ !). Pour réduire le coût de calcul, on observe que les noyaux sont « localisés », c’est à dire leur valeur est numériquement non nulle seulement pour une valeur relativement faible de la norme de la différence entre les arguments ; la valeur du noyau calculé entre un point et une observation très éloignée est ainsi assimilable à 0 et il n’est donc pas nécessaire de faire le calcul. Il est alors possible d’utiliser un index spatial (ici de type KD tree ou Ball tree, algorithm ['kd_tree'|'ball_tree'|'auto']) pour éviter de calculer le noyau entre données trop éloignées les unes des autres.

Les méthodes qui peuvent être employées :

  • fit(X, y=None) : calcul du modèle à partir des observations qui sont les lignes de X.
  • sample([n_samples, random_state]) : générer un échantillon (composé de n_samples données) à partir du modèle (pour le moment, accessible seulement pour les noyaux 'gaussian' et 'tophat').
  • score(X, y=None) : retourne la log-vraisemblance totale des données de X par rapport au modèle.
  • score_samples(X) : retourne le logarithme de la densité calculée pour chaque donnée de X.
  • get_params([deep]) : lire les valeurs des paramètres de l’estimateur employé.
  • set_params(**params) : donner des valeurs aux paramètres de l’estimateur employé.

Il n’y a pas d’attributs publics, pourquoi ?

Estimation à partir de données générées

Cette partie de la séance a pour objectif de faire comprendre l’impact des différents aspects qui interviennent dans l’estimation par noyaux : le volume de données, le choix du noyau, le paramètre d’échelle (bandwidth, correspondant en général à la variance du noyau). Pour cela, nous générons des données à partir d’une densité connue, estimons la densité à partir de ces données et comparons le résultat de l’estimation à la densité qui a généré les données.

Une première utilisation sur des données unidimensionnelles, sur lesquelles il est possible d’afficher à la fois la vraie densité et la densité estimée :

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.stats import norm
>>> from sklearn.neighbors import KernelDensity

# générer l'échantillon à partir de deux lois normales
>>> N = 100
>>> X = np.concatenate((np.random.normal(0, 1, int(0.3 * N)),
                        np.random.normal(5, 1, int(0.7 * N))))[:, np.newaxis]

# préparer les points où on calculera la densité
>>> X_plot = np.linspace(-5, 10, 1000)[:, np.newaxis]

# préparation de l'affichage de la vraie densité, qui est celle à partir
#  de laquelle les données ont été générées (voir plus haut)
# la pondération des lois dans la somme est la pondération des lois
#  dans l'échantillon généré (voir plus haut)
>>> true_density = (0.3*norm(0,1).pdf(X_plot[:,0]) + 0.7*norm(5,1).pdf(X_plot[:,0]))

# estimation de densité par noyaux gaussiens
>>> kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(X)

# calcul de la densité pour les données de X_plot
>>> density = np.exp(kde.score_samples(X_plot))

# affichage : vraie densité et estimation
>>> fig, ax = plt.subplots()
>>> ax.fill(X_plot[:,0], true_density, fc='b', alpha=0.2, label='Vraie densité')
>>> ax.plot(X_plot[:,0], density, '-', label="Estimation")
>>> ax.plot(X[:, 0], -0.005 - 0.01 * np.random.random(X.shape[0]), '+k')
>>> ax.legend(loc='upper left')
>>> plt.show()

Question :

Faites varier la taille de l’échantillon et examinez les résultats.

Question :

Faites varier le nombre de lois normales lors de la génération de l’échantillon et examinez les résultats.

Question :

Faites varier la valeur de bandwidth et examinez les résultats.

Question :

Faites varier le type de noyau et examinez les résultats.

Il est possible de générer de nouvelles données suivant cette densité estimée et de visualiser le résultat d’une nouvelle estimation faite à partir de ces dernières données :

>>> Xg = kde.sample(N)
>>> kde2 = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(Xg)
>>> fig, ax = plt.subplots()
>>> ax.fill(X_plot[:,0], true_density, fc='b', alpha=0.2, label='Vraie densité')
>>> ax.plot(X_plot[:,0], density, '-', label="Estimation")
>>> ax.plot(X_plot[:,0], np.exp(kde2.score_samples(X_plot)), 'r-', label="Estimation2")
>>> ax.legend(loc='upper left')
>>> plt.show()

Une deuxième utilisation de l’estimation de densité par noyaux sur des données bidimensionnelles, sur lesquelles il n’est malheureusement plus possible de bien visualiser à la fois la vraie densité et la densité estimée :

# générer l'échantillon
>>> N = 20
>>> kd = np.random.rand(N, 2)

# définir la grille pour la visualisation
>>> grid_size = 100
>>> Gx = np.arange(0, 1, 1/grid_size)
>>> Gy = np.arange(0, 1, 1/grid_size)
>>> Gx, Gy = np.meshgrid(Gx, Gy)

# définir la largeur de bande pour le noyau
>>> bw = 0.02

# estimation, puis calcul densité sur la grille
>>> kde3 = KernelDensity(kernel='gaussian', bandwidth=bw).fit(kd)
>>> Z = np.exp(kde3.score_samples(np.hstack(((Gx.reshape(grid_size*grid_size))[:,np.newaxis],
        (Gy.reshape(grid_size*grid_size)[:,np.newaxis])))))

# affichage
>>> from mpl_toolkits.mplot3d import Axes3D
>>> from matplotlib import cm
>>> fig = plt.figure()
>>> ax = fig.gca(projection='3d')
>>> ax.plot_surface(Gx, Gy, Z.reshape(grid_size,grid_size), rstride=1,
                    cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=True)
>>> ax.scatter(kd[:,0], kd[:,1], -10)
>>> plt.show()

Question :

Faites varier la valeur de bandwidth et examinez les résultats.

Question :

Générez l’échantillon à partir d’une somme de 2 lois normales différentes (avec N plus grand) et examinez les résultats.

Pour mieux comprendre ce qu’est une estimation de densité par noyaux, il peut être utile de construire le modèle et visualiser les résultats directement, sans passer par la classe spécifique de Scikit-learn. En effet, dans l’estimation de densité par noyaux il n’y a pas de « modèle synthétique » obtenu à partir des N données (ou observations) de départ, le modèle de densité est constitué par les N noyaux centrés sur ces N données.

Pour la visualisation directe nous nous servirons de la fonction de calcul d’une loi normale bidimensionelle en (x,y), bivariate_normal(x,y,sigmax,sigmay,mux,muy,sigmaxy), où sigmax est la variance de X, sigmay la variance de Y, sigmaxy la covariance entre X et Y, mux la moyenne de X et muy la moyenne de Y.

# générer l'échantillon
>>> N = 10
>>> kd = np.random.rand(N, 2)

# définir la grille pour la visualisation
>>> grid_size = 100
>>> Gx = np.arange(0, 1, 1/grid_size)
>>> Gy = np.arange(0, 1, 1/grid_size)
>>> Gx, Gy = np.meshgrid(Gx, Gy)
>>> bw = 0.01

# calcul direct de la densité sur la grille
>>> from matplotlib.mlab import bivariate_normal
>>> Z = np.zeros((grid_size,grid_size))
>>> for i in range(10):
        Z = Z + bivariate_normal(Gx,Gy,bw,bw,kd[i,0],kd[i,1],0)

# affichage
>>> from mpl_toolkits.mplot3d import Axes3D
>>> from matplotlib import cm
>>> fig = plt.figure()
>>> ax = fig.gca(projection='3d')
>>> ax.plot_surface(Gx, Gy, Z, rstride=1, cstride=1, cmap=cm.coolwarm,
                    linewidth=0, antialiased=True)
>>> ax.scatter(kd[:,0], kd[:,1], -500)
>>> plt.show()

Question :

Faites varier la taille de l’échantillon et examinez les résultats. Quelles sont les conséquences d’une augmentation de N ?

Estimation à partir des données « textures »

Pour rappel, ces données correspondent à 5500 observations décrites par 40 variables. Chaque observation appartient à une des 11 classes de textures ; chaque classe est représentée par 500 observations. Les données sont issues de https://www.elen.ucl.ac.be/neural-nets/Research/Projects/ELENA/databases/REAL/texture/. Nous appliquerons l’estimation de densité par noyaux à ces données.

Si les données ne sont pas dans le répertoire de travail il est nécessaire de les chercher d’abord :

[nom@machine ~]$ wget http://cedric.cnam.fr/~crucianm/src/texture.dat

Pour visualiser facilement les résultats, il faudra travailler sur des projections unidimensionnelles ou bidimensionnelles de ces données.

Question :

Appliquez l’analyse en composantes principales à ces données et estimez la densité des données projetées sur le premier axe principal (ou les deux premiers axes principaux). Faites varier la largeur de « fenêtre » (bandwidth). Visualisez les résultats.

Question :

Appliquez l’analyse discriminante à ces données et estimez la densité des données projetées sur le premier axe discriminant (ou les deux premiers axes discriminants). Faites varier la largeur de « fenêtre » (bandwidth). Visualisez les résultats.

Modèles de mélange gaussiens

La présentation des outils disponibles pour les modèles de mélange se trouve dans http://scikit-learn.org/stable/modules/mixture.html.

La description de l’implémentation de l’estimation de densités par noyaux se trouve dans http://scikit-learn.org/stable/modules/generated/sklearn.mixture.GaussianMixture.html.

La sélection d’un modèle (le choix du nombre de composantes du mélange, mais aussi du type de matrice de covariances entre 'full', 'tied', 'diag', 'spherical', voir plus bas) peut être faite grâce au calcul de Akaike information criterion (AIC) ou de Bayes information criterion (BIC).

Parmi les paramètres, arrêtons-nous aux suivants (regarder la description de l’implémentation pour les autres) :

  • n_components : le nombre de composantes du mélange (1 par défaut).
  • covariance_type : type de la matrice de covariance employée, choix entre 'full', 'tied', 'diag' et 'spherical', par défaut 'full' ; 'full' = matrice quelconque, 'tied' = matrice quelconque mais identique entre les différentes composantes du mélange, 'diag' = matrice diagonale (variances quelconques mais covariances nulles), 'spherical' = chaque composante a sa propre variance (mais cette valeur est partagée par toutes les variables pour cette composante) et covariances nulles.
  • init_params : méthode d’initialisation des paramètres, choix entre 'kmeans' (classification automatique des observations avec K-means, ensuite utilisation de chaque centre de groupe comme moyenne d’une composante du mélange et de la matrice de covariances empiriques du groupe comme matrice de variances-covariances pour cette composante) et 'random' ; par défaut 'kmeans'.
  • n_init : nombre d’initialisations (suivies d’exécutions de EM) effectuées ; les meilleurs résultats sont conservés.

Parmi les attributs accessibles nous pouvons mentionner (regardez la documentation pour voir la totalité des attributs) :

  • weights_ : les coefficients de mélange (pondérations des différentes composantes du mélange).
  • means_ : les moyennes des composantes du mélange.
  • covariances_ : les matrices de covariance des composantes du mélange.
  • converged_ : True si EM a convergé (dans .fit()), False sinon.
  • lower_bound_ : log-vraisemblance atteinte à la fin des itérations par la meilleure exécution de EM.

Les méthodes qui peuvent être employées :

  • fit(X, y=None) : calcul du modèle à partir des observations qui sont les lignes de X.
  • aic(X) : valeur du critère d’information de Akaike sur les données de X pour le modèle courant.
  • bic(X) : valeur du critère d’information de Bayes sur les données de X pour le modèle courant.
  • predict_proba(X) : probabilités a posteriori par rapport à chaque composante du modèle courant pour chaque donnée de X.
  • predict(X, y=None) : étiquettes de groupe (classification automatique issue du modèle de mélange grâce aux probabilités a posteriori) obtenues avec le modèle courant pour les données de X.
  • sample([n_samples, random_state]) : générer un échantillon (composé de n_samples données) à partir du modèle (pour le moment, accessible seulement pour les noyaux 'gaussian' et 'tophat').
  • score(X, y=None) : retourne la log-vraisemblance totale des données de X par rapport au modèle.
  • score_samples(X) : retourne le logarithme de la densité calculée pour chaque donnée de X.
  • get_params([deep]) : lire les valeurs des paramètres de l’estimateur employé.
  • set_params(**params) : donner des valeurs aux paramètres de l’estimateur employé.

Estimation à partir de données générées

Une première utilisation sur des données unidimensionnelles similaires à celles employées pour l’estimation par noyaux (attention, la classe GaussianMixture est disponible à partir de la version 0.18 de Scikit-learn, pour des versions antérieures il faut se servir plutôt de la classe GMM) :

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.stats import norm
>>> from sklearn.mixture import GaussianMixture

# générer l'échantillon
>>> N = 100
>>> X = np.concatenate((np.random.normal(0, 1, int(0.3 * N)),
           np.random.normal(5, 1, 0.7 * N)))[:, np.newaxis]

# préparer les données où on calculera la densité
>>> X_plot = np.linspace(-5, 10, 1000)[:, np.newaxis]

>>> true_density = (0.3*norm(0,1).pdf(X_plot[:,0]) + 0.7*norm(5,1).pdf(X_plot[:,0]))

# estimation par mélange gaussien, avec le « bon » nombre de composantes
>>> gmm = GaussianMixture(n_components=2,n_init=3).fit(X)
>>> gmm.converged_
True
>>> gmm.lower_bound_
-1.8839063293411455

# calcul de la densité pour les données de X_plot
>>> density = np.exp(gmm.score_samples(X_plot))

# affichage : vraie densité et estimation
>>> fig, ax = plt.subplots()
>>> ax.fill(X_plot[:,0], true_density, fc='b', alpha=0.2, label='Vraie densité')
>>> ax.plot(X_plot[:,0], density, '-', label="Estimation")
>>> ax.plot(X[:, 0], -0.005 - 0.01 * np.random.random(X.shape[0]), '+k')
>>> ax.legend(loc='upper left')
>>> plt.show()

Question :

(\rightarrow Compte-rendu) Faites varier le nombre de composantes et examinez visuellement les résultats. Regardez les pondérations des composantes et les moyennes de ces composantes. Examinez de quelle manière évolue la valeur finale atteinte par la log-vraisemblance avec le nombre de composantes.

Une utilisation sur des données bidimensionnelles :

# générer l'échantillon
>>> md1 = 1.5 * np.random.randn(200,2) + [3,3]
>>> md2 = np.random.randn(100,2).dot([[2, 0],[0, 0.8]]) + [-3, 3]
>>> md3 = np.random.randn(100,2) + [3, -3]
>>> md = np.concatenate((md1, md2, md3))

# préparer les données où on calculera la densité
>>> grid_size = 100
>>> Gx = np.arange(-10, 10, 20/grid_size)
>>> Gx.shape
(100,)
>>> Gy = np.arange(-10, 10, 20/grid_size)
>>> Gx, Gy = np.meshgrid(Gx, Gy)
>>> Gx.shape
(100, 100)

# estimation par mélange gaussien
>>> gmm = GaussianMixture(n_components=3,n_init=3).fit(md)
>>> gmm.converged_
True
>>> gmm.lower_bound_
-4.3399960707653937

# calcul de la densité pour les données de la grille
>>> density = np.exp(gmm.score_samples(np.hstack(((Gx.reshape(grid_size*grid_size))[:,np.newaxis],
        (Gy.reshape(grid_size*grid_size)[:,np.newaxis])))))

# affichage : données et estimation
>>> from mpl_toolkits.mplot3d import Axes3D
>>> from matplotlib import cm
>>> fig = plt.figure()
>>> ax = fig.gca(projection='3d')
>>> ax.plot_surface(Gx, Gy, density.reshape(grid_size,grid_size), rstride=1,
                    cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=True)
>>> ax.scatter(md[:,0], md[:,1], -0.025)
>>> plt.show()

Question :

Faites varier le nombre de composantes et examinez visuellement les résultats. Examinez de quelle manière évolue la valeur finale atteinte par la log-vraisemblance avec le nombre de composantes.

Question :

Générez des données bidimensionnelles suivant une distribution uniforme dans [0, 1]^2 et estimez un mélange gaussien avec 3 composantes. Visualisez les résultats. Utilisez la méthode predict de GaussianMixture pour obtenir des étiquettes de groupe pour les données. Appliquez plusieurs fois de suite la modélisation suivie d’affectation d’étiquettes de groupe et examinez la stabilité des partitionnements en utilisant l’indice de Rand ajusté, comme dans les TP de classification automatique. Que constatez-vous ?

Choix du nombre de composantes et du type de matrice de covariances pour les données générées

Pour choisir le nombre de composantes du mélange, nous comparerons d’abord les critères AIC et BIC en utilisant des matrices de covariance 'full' (par défaut). Sur les données bidimensionnelles générées :

>>> n_max = 8    # nombre de valeurs pour n_components
>>> n_components_range = np.arange(n_max)+1
>>> aic = []
>>> bic = []

# construction des modèles et calcul des critères
>>> for n_comp in n_components_range:
>>>     gmm = GaussianMixture(n_components=n_comp).fit(md)
>>>     aic.append(gmm.aic(md))
>>>     bic.append(gmm.bic(md))

>>> aic
[3634.7623945384089, 3505.7538502557586, 3512.408840477905, 3521.0648517705217,
3533.2915729309807, 3545.8605525753455, 3545.5083205836022, 3535.7598532705847]

# normalisation des résultats obtenus pour les critères
>>> raic = aic/np.max(aic)
>>> rbic = bic/np.max(bic)

# affichage sous forme de barres
>>> xpos = np.arange(n_max)+1  # localisation des barres
>>> largeur = 0.3              # largeur des barres
>>> plt.ylim([min(np.concatenate((rbic,raic)))-0.1, 1.1])
>>> plt.bar(xpos, raic, largeur, color='r', label="AIC")
>>> plt.bar(xpos+largeur, rbic, largeur, color='b', label="BIC")
>>> plt.legend(loc='upper left')
>>> plt.show()

Pour ces données et avec des matrices de covariances 'full', à la fois AIC et BIC permettent de privilégier l’utilisation de 3 composants, nombre égal à celui de lois normales qui ont permis de générer les données.

Question :

(\rightarrow Compte-rendu) Réalisez la même comparaison pour ces données avec des matrices de covariances 'diag'.

Question :

Ajoutez sur le même graphique (adapté) les valeurs normalisées de la log-vraisemblance finale pour chaque valeur de n_components.

Estimation à partir des données « textures »

Nous appliquerons l’estimation de densité par mélange gaussien aux données « textures » projetées sur les deux premiers axes principaux.

# lecture des données et aplication de l'ACP
>>> from sklearn.decomposition import PCA
>>> textures = np.loadtxt('texture.dat')
>>> pca = PCA().fit(textures[:,:40])
>>> texturesp = pca.transform(textures[:,:40])

# construction du modèle de mélange, vérifications
>>> gmm = GaussianMixture(n_components=11,n_init=3).fit(texturesp[:,:2])
>>> gmm.converged_
True
>>> gmm.n_iter_
10
>>> gmm.lower_bound_
-1.7780023505669722

# préparer les données où on calculera la densité
>>> grid_size = 100
>>> xmin = 1.3*np.min(texturesp[:,0])
>>> xmax = 1.3*np.max(texturesp[:,0])
>>> Gx = np.arange(xmin, xmax, (xmax-xmin)/grid_size)
>>> ymin = 1.3*np.min(texturesp[:,1])
>>> ymax = 1.3*np.max(texturesp[:,1])
>>> Gy = np.arange(ymin, ymax, (ymax-ymin)/grid_size)
>>> Gx, Gy = np.meshgrid(Gx, Gy)

# calcul de la densité pour les données de la grille
>>> density = np.exp(gmm.score_samples(np.hstack(((Gx.reshape(grid_size*grid_size))[:,np.newaxis],
        (Gy.reshape(grid_size*grid_size)[:,np.newaxis])))))

# affichage des résultats
>>> from mpl_toolkits.mplot3d import Axes3D
>>> from matplotlib import cm
>>> fig = plt.figure()
>>> ax = fig.gca(projection='3d')
>>> ax.plot_surface(Gx, Gy, density.reshape(grid_size,grid_size), rstride=1,
                    cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=True)
>>> ax.scatter(texturesp[:,0], texturesp[:,1], -0.25)
>>> plt.show()

Question :

Appliquez l’estimation de densité par mélange gaussien sur les données « textures » projetées sur les deux premiers axes discriminants.

Choix du nombre de composantes et du type de matrice de covariances pour les données « textures »

Déterminons d’abord le meilleur nombre de composantes en utilisant les critères AIC et BIC pour les projections de ces données sur les deux premiers axes principaux :

>>> n_max = 8    # nombre de valeurs pour n_components
>>> n_components_range = np.arange(n_max)+6
>>> n_components_range
array([ 6,  7,  8,  9, 10, 11, 12, 13])
>>> aic = []
>>> bic = []

# construction des modèles et calcul des critères
>>> for n_comp in n_components_range:
...     gmm = GaussianMixture(n_components=n_comp).fit(texturesp[:,:2])
...     aic.append(gmm.aic(texturesp[:,:2]))
...     bic.append(gmm.bic(texturesp[:,:2]))
...

# normalisation des résultats obtenus pour les critères
>>> raic = aic/np.max(aic)
>>> rbic = bic/np.max(bic)

# affichage sous forme de barres
>>> xpos = n_components_range  # localisation des barres
>>> largeur = 0.3              # largeur des barres
>>> plt.ylim([min(np.concatenate((rbic,raic)))-0.02, 1.01])
>>> plt.bar(xpos, raic, largeur, color='r', label="AIC")
>>> plt.bar(xpos+largeur, rbic, largeur, color='b', label="BIC")
>>> plt.legend(loc='upper left')
>>> plt.show()

Question :

(\rightarrow Compte-rendu) Comment expliquez-vous que la valeur optimale obtenue avec BIC pour n_components soit inférieure au nombre de classes de textures (qui est de 11) ?

Question :

(\rightarrow Compte-rendu) Réalisez la même étude en utilisant plutôt les projections des données « textures » sur les deux premiers axes discriminants. Comment expliquez-vous la différence par rapport au cas précédent ?