Algorithme EM - Application à des données non simulées

Dans ce TP, nous allons étudier des données issues de l'article suivant :

  • Bachrach LK, Hastie T, Wang M-C, Narasimhan B, Marcus R. Bone Mineral Acquisition in Healthy Asian, Hispanic, Black and Caucasian Youth. A Longitudinal Study. J Clin Endocrinol Metab (1999) 84, 4702-12.

Ces données son disponibles ici. Récupérez les et placez les dans le même répertoire que votre notebook. Elles représentent des mesures relatives de densité minérale osseuse spinale sur des adolescents nord-américains. Chaque valeur est la différence de mesures prises sur deux visites consécutives, divisées par la moyenne.

L'idée est de savoir si la population peut être décrite par un mélange de deux gaussiennes. Le but est donc d'appliquer à ces données l'algorithme EM, vérifier que le $K$ optimal de la méthode de sélection de modèle vue en cours est bien $K=2$ et proposer un clustering de ces observations.

In [1]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as sps
import scipy.linalg as spl

On commence par créer le tableau de donnees "data" à partir du fichier.

In [3]:
data=np.loadtxt("../datasets/densitesOs.txt")
n,d = data.shape

Données

1) De quelle dimension sont les données ? Affichez les données sous la forme d'un nuage de points.

In [4]:
 

Algorithme EM

2) On va maintenant appliquer l'algorithme EM aux données précédentes.

  • Ecrire une fonction EM(K,nbpas,data) prenant en entrée les données 'data', un entier K (nombre de gaussiennes du mélange) et un nombre de pas nbpas.
    • Appliquer la fonction précédente aux données pour les paramètres du GMM sous-jacent pour un nombre de classes K donné.

On fixe le nombre de pas de EM et un intervalle pour les valeurs de K que l'on va tester.

In [5]:
nbpas=int(5e1) # nombre de pas dans EM
rangeK=range(1,5) # valeurs de K testees pour le nombre de classes (dans la selection de modele)

Sélection de modèle

Pour la sélection de modèle, il faut calculer la log-vraisemblance et la complexité du modèle.

Une manière classique de choisir le nombre de composantes $K$ dans un intervalle d'entiers est de minimiser en $K$ $$-\ell(x;\widehat{\theta}_K)+ \frac{\operatorname{dim}_K \times \log n}{2},$$ avec

  • $\widehat{\theta}_K$ l'estimateur de $\theta$ issu de l'algorithme EM dans le modèle de mélange avec $K$ classes,
  • $\ell(x;\hat{\theta}_K)$ la log-vraisemblance de l'échantillon observé $(x_1,\dots , x_n)$
  • $\operatorname{dim}_K$ le nombre de degrés de liberté du modèle de mélange utilisé avec $K$ classes (voir le cours).

Le K optimal pour la sélection de modèle est l'argmin de cette fonction.

3) Ecrire une fonction logLikelyhood(K,alpha,mu,Cov,data) permettant de calculer la vraisemblance du GMM donné par les paramètres (K,alpha,mu,Cov) sur les données data.

4) Implémenter la sélection de modèle et l'appliquer aux données précédentes pour déterminer le K optimal. Commentez le résultat.

5) Afficher les données en les colorant en fonction de leur classe telle qu'estimée par EM.

6) Afficher la fonction $-\ell(x;\widehat{\theta}_K)+ \frac{\operatorname{dim}_K \times \log n}{2},$ en fonction de $K$

In [9]:
 

En utilisant la librarie sklearn

Refaire la même chose en utilisant la fonction sklearn.mixture.GaussianMixture de la bibliothèque scikit-learn. Voir la page suivante

In [11]:
import sklearn.mixture