La plupart des mesures que l'on obtient des expériences en biologie suivent approximativement une distribution dite "normale", ou "gaussienne", dont la densité a la forme d'une cloche, symétrique avec un unique sommet au milieu. C'est aussi l'hypothèse d'un grand nombre d'outils d'analyse statistique. Mais que faire quand on observe deux sommets ou plus ? Le plus probable, c'est qu'on observe alors un mélange de plusieurs composantes normales, qu'on voudrait séparer. Une fois décomposées, on peut travailler sur chaque composante séparément. Le tout s'appelle un "modèle de mélange gaussien", ou "gaussian mixture model" en anglais pour la littérature.
Les exemples d'utilisation sont fréquents. D'ailleurs, si j'ai choisi d'écrire sur ce sujet, c'est que j'ai reçu il y a quelques jours le message suivant :
Salut,
J'ai un service à te demander ! J'ai une série de 100 points qui représentent des distances d'expression de gène entre deux cellules après division. J'obtiens une distribution qui n'est pas une gaussienne, et j'aimerais savoir s'il y a moyen de trouver une double gaussienne qui pourrait modéliser de façon adéquate ma distribution. Le but étant de fixer un seuil pour les paires de cellules avec une distance élevée. Est-ce que tu pourrais me montrer comment faire ça ? C'est assez urgent…
A.
Pour déterminer quel individu appartient à quelle composante, la méthode usuelle est d'utiliser un algorithme EM (pour Expectation-Maximization). Le résultat de l'algorithme est un ensemble de paires (moyenne, variance) — une pour chaque composante normale — , et une probabilité "a posteriori" pour chaque observation d'appartenir à telle ou telle composante.
Sans entrer dans les détails, l'algorithme fonctionne en alternant deux phases : on attribue — arbitrairement, la première fois — chaque point à l'une des composantes (phase E), dont on peut alors estimer la densité de probabilité (phase M). Puis on considère fixes les densités, et on redistribue au mieux les points entre elles (phase E), par un critère de maximum de vraisemblance. Et on recommence : on recalcule les densités (phase M), etc. On peut montrer que chaque itération améliore la vraisemblance, et donc qu'on converge toujours vers un maximum local.
Voici comment faire en utilisant le langage R. Pour ne pas s'encombrer avec l'import de données réelles, générons des données aléatoires, à partir de deux distributions normales, mélangées dans un même vecteur :
Mais comme les points sont difficilement visibles, on préfère en principe dessiner un histogramme ou une courbe de densité, comme ça :
Histogramme et courbe de densité (non-paramétrique) des données.
On devine qu'il doit y avoir deux composantes (deux "bosses"). Il sera important de préciser à l'algorithme combien de composantes il y a, sans quoi le résultat sera n'importe quoi.
Ensuite il existe de nombreux packages à même de faire le travail. Voilà un exemple avec la librairie R "mixtools" :
En rouge la densité jointe (non-paramétrique) ; en bleu/vert les courbes de densité des deux composantes normales, d'après les estimations de mixtools.
Et voilà, on a séparé les données en deux composantes normales. On peut ajouter à volonté au graphe la position des moyennes, les écarts-types, les points colorés selon leur appartenance à l'une ou l'autre composante, etc.
Alternativement, en utilisant la librairie "mclust" et en utilisant d'autres possibilités graphiques :
On a ajouté à l'histogramme des données la densité jointe (non-paramétrique, en rouge), la moyenne (bleu plein) et la variance (bleu trait-tillé) de chaque composante normale.
Généralisation
Notez que si dans notre exemple nous avons deux composantes normales à une dimension, on peut facilement généraliser à X composantes en N dimensions (comme dans l'exemple ci-dessous) mais aussi avec d'autres sortes de distributions de probabilité. Dans le cas de deux dimensions, la ligne de points de notre toute première figure ci-dessus devient un nuage de points dans un plan.
A gauche les données initiales ; on devine trois sous-groupes. A droite, la décomposition par mixtools, où on a coloré les points selon la composante à laquelle ils appartiennent le plus probablement.
Remarque : le calcul a tout de même pris une bonne dizaine de minutes sur mon vieux Mac.
Il existe dans à peu près tous les langages des librairies de ce genre. Un exemple notable est le tout prometteur langage Julia, encore en développement, qui possède une librairie "MixtureModels", par exemple :
r.mixture.components#liste des composantes normales
r.mixture.pi#proportions
r.Q#probabilitésaposteriori
Parfois c'est plus compliqué : l'algorithme EM peut être implémenté dans sa forme générale (il a beaucoup d'autres usages), et c'est à vous de spécifier à quoi l'appliquer. Mais c'est un autre sujet qui mérite un livre à lui tout seul.
Bioinformaticien @ CHUV / UNIL
Diplômé EPFL en mathématiques appliquées.
Développement software, analyse de données génomiques.
Pour continuer la lecture :
Commentaires
34 réponses à “Les mélanges gaussiens”
prince
Bonjour Mr. Julien Delafontaine,
je trouve ce tutoriel très enrichissant et aimerais vous poser une question par rapport à ma thèse de master que j'ecris en ce moment.
Le thème est le suivant : Recherche sur le em algorithme.
je suis ainsi à la recherche des exemples d'applications du modèle.
j'ai deja fait un exemple pour les valeurs manquantes avec R et aimerais en faire 2 autres exemples.
auriez vous d'autres cas d'apllications ?
Bonjour,
Oui, il y a beaucoup d'autres applications de l'algo EM.
Très souvent c'est en probabilités, pour un calcul de maximum de vraisemblance (voir en particulier Monte-Carlo EM).
Pour un usage plus spécifique en bioinfo, il y a l'étude de l'épissage alternatif (splicing) pour le RNA-seq : des programmes comme RSEM ou Cufflinks en implémentent un.
Enfin dans un autre domaine, le Filtre de Kalman est une application très importante (GPS par exemple).
merci pour cette reponse rapide et les exemples mentionnes.
en image segmentation on pourrait aussi apparenment lutiliser mais je narrive pas a comprendre cela.
En plus mon profeseur souherait avoir des codes dont il pourrait aussi tester a la maison.
auriez vous des exemples , codes y compris qui me permettront de mieux comprendre.
merci
Je n'en ai pas directement, mais on peut trouver du code sur le net, comme en cherchant "kalman filter EM example", on tombe sur du code R qui fait ca : http://www.econ.umn.edu/~karib003/help/kalman_example2.htm
(SSM = State Space Model)
Bonjour Monsieur Julien Delafontaine,
Merci pour ce tutoriel qui ma beaucoup aidé dans mon travail, je souhaiterais vous posez quelques questions sur la performance de cette approche lorsque le problème se complique, je veux dire lorsque les composantes deviennent plus imbriqués. Comment peut-on étudier par Monte-Carlo la précision de EM (variance des estimateurs)?
En exécutant votre programme, j'ai aussi remarqué que dans le cas multivarié l’algorithme est lent, je voudrais savoir à quoi cela est du ? Est ce que la taille de l'échantillon y est pour quelques chose ou bien c'est le nombre des composants (ici 3)?
Je souhaite faire une comparaison entre l'algorithme EM et k‑mean mais je viens de débuter en R (je travaille avec SAS) je souhaite savoir quelle est la commande(ou bibliothèque) qui me permettra de faire cela?(en passant je ne comprend pas pourquoi la bibliothèque(mclust) ne marche pas chez moi,il n'existe pas dans mon aide).
Infiniment merci pour ce tutoriel.
cordialement.
Pour étudier la précision de la méthode, on peut générer soi-même des données comme je l'ai fait (Monte-Carlo c'est juste une méthode pour le sampling d'une distribution donnée), par exemple un certain nombre de composantes normales. Comme on connaît par construction l'appartenance de chaque point a l'un des générateurs, on peut la comparer à ce que retrouve l'algorithme EM sans cette information.
Je ne sais pas exactement pourquoi il est lent, mais ça peut être parce que
— c'est écrit en R
— il lui faut beaucoup d'itérations pour converger
— (?) l'algorithme est intrinsèquement lent (demande beaucoup d'opérations au CPU)
Si c'est pour du clustering, les k‑means seront plus efficaces, d'autant que dans les deux cas il faut deviner en avance le nombre de composantes. En R, c'est la fonction "kmeans" qui est présente sans package particulier.
Mclust au contraire doit être installé au préalable : install.packages("mclust").
En considérant les données bivariees (eruptions, waiting) de la table faithful, comment peut-on appliquer la méthode ci-dessus ? En appliquant les codes j'ai l'impression que la fonction mvnormalmixEM marchent de façon infinie.
Je ne connais pas ces tables. Si c'est lent c'est normal. Soit le jeu de donnees est tres grand, soit l'algorithme met beaucoup de temps a converger.
dans le cadre de travail, j'aimerais aussi faire une comparaison entre le k‑means et Lalgorithme EM. Jusqu'à présent mes essais sont encore loin detre concluants. Aurais tu deja des meilleures idées a propos.
Mr. Julien Delafontaine votre aide sera aussi la bienvenue.
je compte en fait generer des donnees dans R et y appliquer les 2 methodes et comparer les resultats.
Merci d'avance
La différence c'est que l'EM va retourner les paramètres des gaussiennes, alors que le K‑means va seulement donner des groupes. Par ailleurs le modèle de mélange gaussien est appliquable comme son nom l'indique uniquement aux données supposées gaussiennes.
Le plus simple en effet est de "generer des donnees dans R et y appliquer les 2 methodes et comparer les resultats". Dans l'article on montre comment générer les données.
Comme a dit ci haut Mr Julien Delafontaine,constate que même si le k‑means souffre d’un cadre théorique rigoureux pour le choix des classes, il donne des résultats plus satisfaisant que le EM. Avec un grand nombre de données, cet algorithme est rapide, performant et efficace. Comparativement à EM, quand le nombre de classe devient plus imbriqué (m=3), ont attend au moins huit minutes pour qu’il converge et un essai avec quatre mélange ne fait que confirme cette lenteur. Toute fois, elle permet de donner un cadre théorique et d’être ainsi une approche qui viendra renforcer le k‑means. La lenteur de EM peut s'expliquer par le fait que la maximisation de la log-vraisemblance qui intervient dans le calcul des paramètres se fait linéairement.Ainsi lorsque le EM atteint un maximun local il peut converger alors, ce maximum n'est pas le maximum global recherché.
Bonjour,
Je suis très intéressée par votre méthodologie pour étudier ces mélanges gaussiens. J'aurais aimé savoir s'il y existe selon vous un moyen pour identifier le point d'iso-probabilité, c'est à dire le point délimitant les deux composants du mélange (entre les deux composants)?
Re !
En fait voilà j'aimerais effectuer une analyse de mélanges sur la distribution d'une variable. J'aimerais voir s'il est possible d'identifier au sein de cette distribution deux composants. Lorsque j'applique le script R présenté ci-dessus via le package mixtool j'obtiens les proportions de chaque composant, les moyennes…Mais j'aimerais pouvoir savoir :
— à quel groupe chaque observation appartient
— quel est le point à partir duquel on passe d'un composant à l'autre
— quel est le taux d'erreur (soit le taux indiquant la probabilité de placer une observation dans le mauvais composant)
Bonjour,
Ce sont les probabilités a posteriori, "mixmodel$posterior" dans l'exemple. On ne peut pas attribuer avec certitude une observation à une distribution, par définition. Par contre on peut choisir, par exemple dans le cas d'un mélange de deux gaussiennes, d'attribuer à chacune les observations qui ont une probabilité > 0.5 de puis appartenir. C'est ce qui a été fait pour mettre les couleurs dans le graphique illustrant le cas à deux dimensions. Le point d'iso-probabilité est l'endroit où les deux courbes de densité se croisent.
Bonjour,
Merci beaucoup ! Oui c'est ce que je pensais avoir compris concernant le point d'iso-probabilité mais je me demandais s'il y avait un moyen de le calculer directement au lieu de l'identifier graphiquement.
J'ai une autre question concernant l'allure du graphique obtenu…En fait sur la représentation graphique obtenue l'une des deux courbes du mélange est vraiment "plate". Dans ce cas est-ce qu'il est juste de dire qu'on a pu identifier deux distributions gaussiennes mais que l'une des deux distribution regroupe un très faible nombre d'observations ? Ou bien la forme plate de cette deuxième composante empêche de conclure sur l'existence réelle de deux composantes et il vaut mieux retester le modèle en prenant une autre valeur de k ?
Par avance merci
Il faut utiliser cette méthode si on sait 1. qu'on a des populations normales, et 2. que d'après l'histogramme il n'y en a pas qu'une (même si c'est "plat"). Si c'est "plat" pour une autre raison, c'est qu'on n'a pas une distribution normale, et le modèle de mixture gaussienne est inapproprié.
Ces deux critères découlent du contexte du problème, et donc c'est entièrement la responsabilité du chercheur de les deviner. Dans l'exemple j'ai généré des points à partir de deux lois normales, donc il est certain que le modèle est approprié.
On peut calculer le point d'iso-probabilité en connaissant les paramètres des deux normales, car on a une formule pour ces deux courbes (voir Wiki "loi normale", celle avec l'exponentielle). On peut écrire un système à deux équations et le résoudre.
Merci beaucoup pour votre blog très instructeur.
Une question ici sur le test d'ajustement qui viendrait si on veut avoir une idée de la précision de l'approximation gaussienne ; le test du chi-deux serait j'imagine le fleuron de l'arsenal classique pour ceci, mais que devient-il si nous travaillons sur des lois bivariées par exemple ?
Et avec les données issues de la réalité l'approximation (même avec un mélange gaussien) est souvent très brute & la probabilité issue du chi-deux est trop faible pour être "montrable" ; auriez-vous idée d'autres tests d'ajustements qui pourraient parvenir à quantifier la précision de l'approximation ? Si vous avez des liens vers d'autres articles je suis bien sûr preneur.
Merci encore & bonne continuation !
Cordialement,
A mon avis si le test retourne des résultats insatisfaisants, le modèle est mauvais et il faut réfléchir à autre chose, pas changer le test jusqu'à obtenir une bonne p‑valeur. Le test a une signification (la somme des carrés des erreurs, par exemple, est tout ce qu'il y a de concret).
On peut comparer différents modèles par leur critère d'AIC (qui découle du rapport des vraisemblances).
Bonjour Julien et bravo pour ce tuto.
J'aimerai savoir s'il y a un moyen de savoir la meilleur valeur de k (modalités).
J'ai vu un article qui parlait de la " MLE error" mais je ne sais pas trop s'il y a d'autres moyens.
Bonjour, Il y a des moyens mathématiques, mais je ne le conseillerais pas. La meilleure méthode c'est de produire un graphique des données pour vérifier si déjà à l'oeil on observe une séparation, comme la courbe de densité ci-dessus. Si on ne voit rien de probant, il ne faut même pas faire de stats du tout.
Si on y tient quand même, ou si on a des normales multi-dimensionnelles, il faut d'abord observer que plus il y a de composantes mélangées, moins le modèle sera performant, donc la plupart du temps on peut juste tester 1,2,3,4 composantes. Enfin, la librairie utilisée devrait retourner un estimateur de vraisemblance, et on peut l'utiliser pour calculer des critères d'accord du modèle aux données (p.ex. Akaike, AIC) ou faire un test des rapports de vraisemblance (likelihood ratio test). Ces critères indiqueront quel est le modèle qui correspond le mieux aux données.
Il n'y a pas de façon par contre de déterminer le nombre idéal de composantes sans appliquer le modèle dans chaque cas et comparer les résultats.
J'ai utiliser la fonction boot.comp sur mon jeu de données pour trouver un nombre de composante qui retranscrit au "mieux la réalité". Par la suite j'utilise normalmixEM pour une simulation sur mon jeu de données. Ma question est la suivante, en fonction du nombre d'opérations, le nb de k sélectionnés peut varier, de plus avec normalmixEM pour la même ligne de code les valeurs mu lambda et sigma peuvent changer également. Comment choisir alors le "meilleur modèle". J'aimerais avoir l'AIC ou le BIC qui existe dans mixtools avec les mixture et autre modele (comparaison des modeles en fonction du nombre de composante) mais pas avec la fonction normalmixEM. Je ne vois pas comment faire, je pensais le calculer manuellement mais comment avoir accès aux nb de résidus, a la somme des carrés des res.…DSL j'espère que je ne suis pas a côté de la plaque.
Je ne comprends pas. Avec normalmixEM(), on choisit soi-meme le k. Si on n'arrive pas à le choisir en voyant l'histogramme, c'est probablement que le modèle de mélange gaussien est un mauvais choix au départ.
En ce qui concerne l'implémentation, je ne me souviens plus précisément et ne me substitue pas à la documentation du package.
Je suis une étudiante qui est en train de travailler sur les mélanges gaussienne.J'avais utilisé votre code pour le réaliser, mais chaque fois quand je lance le programme, je pourrais obtenir des résultats différents. En plus, je trouve que des valeur qu'il a trouvé c'est pas tout du le résultat que je veux. est ce que vous pouvez m'aidez ?
Pour vous aider, il faudrait que vous précisiez ces termes : "le programme", "résultats différents", "des valeur qu'il a trouvé" et "le résultat que je veux".
Bonjour,
Pour la valeur de k, je fais référence aux commentaires ci-dessus.
"Not convergent" veut dire que la condition de termination de l'algorithme n'a pas été atteinte après un grand nombre d'itérations.
Si l'algorithme ne converge pas, c'est probablement que les données ne se prêtent pas à ce modèle (p.ex. avec distribution hors de la famille exponentielle dont la normale fait partie). S'il ne converge pas, le résultat n'est pas utilisable.
Merci pour la reponse. J'ai une autre question à vous poser aussi. Si le k=3 par exemple, j'aimerais bien de trouver un bias de ma système, donc il faut que je concentré sur quelle gaussienne ? Laquelle a une plus grosse surface ?
merci d'avance.
Je suis à la recherche d'une solution pour présenter des résultats issus d'expériences où on mesure les pptés nanomécaniques, notamment le module élastique des bactéries en contact des nanoparticules oxydes. Généralement, on obtient une distribution en double gaussienne. Je suis tombée sur votre tutoriel, mais je n'ai pas réussi à m'en sortir.
Pourriez-vous m'indiquer la marche à suivre. Mes données sont sous forme d'histogrammes (Fréquence en fct de l'élasticité).
Merci par avance
Angelina
Bonjour Julien,
Merci pour le tutoriel.
Je voulais savoir comment faire pour avoir les différents paramètres du modèle.
En fait moi dans le cadre de mes travaux, J'aimerais pourvoir générer les paramètres des modèles car elles ont une signification biologique.
paramètres : optimums, largeurs des courbes …
Pour insérer du code dans vos commentaires, utilisez les balises <code> et <\code>.
Gérer le consentement aux cookies
Pour offrir les meilleures expériences, nous utilisons des technologies telles que les cookies pour stocker et/ou accéder aux informations des appareils. Le fait de consentir à ces technologies nous permettra de traiter des données telles que le comportement de navigation ou les ID uniques sur ce site. Le fait de ne pas consentir ou de retirer son consentement peut avoir un effet négatif sur certaines caractéristiques et fonctions.
Fonctionnel
Toujours activé
Le stockage ou l’accès technique est strictement nécessaire dans la finalité d’intérêt légitime de permettre l’utilisation d’un service spécifique explicitement demandé par l’abonné ou l’internaute, ou dans le seul but d’effectuer la transmission d’une communication sur un réseau de communications électroniques.
Préférences
Le stockage ou l’accès technique est nécessaire dans la finalité d’intérêt légitime de stocker des préférences qui ne sont pas demandées par l’abonné ou la personne utilisant le service.
Statistiques
Le stockage ou l’accès technique qui est utilisé exclusivement à des fins statistiques.Le stockage ou l’accès technique qui est utilisé exclusivement dans des finalités statistiques anonymes. En l’absence d’une assignation à comparaître, d’une conformité volontaire de la part de votre fournisseur d’accès à internet ou d’enregistrements supplémentaires provenant d’une tierce partie, les informations stockées ou extraites à cette seule fin ne peuvent généralement pas être utilisées pour vous identifier.
Marketing
Le stockage ou l’accès technique est nécessaire pour créer des profils d’internautes afin d’envoyer des publicités, ou pour suivre l’internaute sur un site web ou sur plusieurs sites web ayant des finalités marketing similaires.
Laisser un commentaire