Statistique

L’algorithme EM

Deviner ce qui manque, puis réestimer comme si on le savait, et recommencer : l’algorithme « espérance–maximisation » de Dempster, Laird et Rubin, pourquoi il monte toujours, à quelle vitesse, à quoi il sert et où il se trompe.

L’idée en une phrase. Quand une partie de l’information manque (une étiquette perdue, une valeur non relevée, un groupe invisible), le calcul du maximum de vraisemblance devient difficile ; mais il serait facile si l’on connaissait ce qui manque. L’algorithme EM alterne alors deux étapes simples : deviner ce qui manque avec les paramètres actuels, puis réestimer les paramètres comme si ces devinettes étaient vraies.

Le mot vraisemblance désigne la probabilité des données observées, vue comme fonction du paramètre θ ; l’estimateur du maximum de vraisemblance θ̂ est la valeur qui la rend maximale. On travaille avec son logarithme ℓ(θ), la log-vraisemblance. Les données complètes sont les données observées augmentées de ce qui manque ; leur log-vraisemblance est souvent si simple qu’on la maximise de tête.

Étape	Nom	Ce qu’on fait
E	Espérance (Expectation)	avec les paramètres actuels θₜ, on calcule la loi de ce qui manque sachant ce qu’on a vu, et l’on remplace chaque quantité inconnue par sa valeur moyenne (une « responsabilité », un effectif attendu, une valeur prédite). On obtient la fonction Q(θ \| θₜ), espérance de la log-vraisemblance complète.
M	Maximisation	on maximise Q(θ \| θₜ) comme on le ferait avec des données complètes : souvent une simple moyenne ou une fréquence pondérée. Le résultat est θₜ₊₁.
Garantie	Montée	ℓ(θₜ₊₁) ≥ ℓ(θₜ) à chaque tour, sans réglage ni dérivée à calculer (planche 2).
Prix	Lenteur et départ	la convergence est lente quand beaucoup d’information manque (planche 3), et l’on peut finir sur un sommet secondaire (planche 5).

Deux pièces, et des étiquettes perdues

On dispose de deux pièces A et B, de probabilités de pile θ_A et θ_B inconnues. Cinq fois, on prend l’une des deux au hasard et on la lance dix fois. On note les piles… mais on a oublié de noter quelle pièce a servi. Si l’on savait, on estimerait θ_A par la fréquence de pile des séries faites avec A. Sans le savoir, EM procède ainsi : étape E, pour chaque série, calculer la probabilité qu’elle vienne de A avec les valeurs actuelles (sa responsabilité) ; étape M, recalculer θ_A comme une fréquence où chaque série compte pour sa responsabilité, et θ_B de même avec le complément. C’est l’exemple de Do et Batzoglou (2008).

À gauche : les cinq séries (pile en bleu foncé, face en clair) et, à droite de chaque série, la part attribuée à A (bleu) et à B (orange) par l’étape E. À droite : le plan des paramètres (θ_A, θ_B), les lignes de niveau de la log-vraisemblance et le chemin suivi par EM (rouge). Le point vert marque le maximum trouvé par une recherche indépendante.

Départ θ_A Départ θ_B

Itération et paramètres—

Log-vraisemblance ℓ(θₜ)—

Contrôle : maximum par recherche directe—

Si les étiquettes étaient connues0,80 ; 0,45vraies pièces des séries : B, A, A, B, A → 24 piles sur 30 et 9 sur 20

En quelques pas, EM retrouve θ_A ≈ 0,80 et θ_B ≈ 0,52, exactement le maximum de vraisemblance obtenu par une recherche directe sur une grille. Il ne retrouve pas les fréquences 0,80 et 0,45 qu’on aurait avec les étiquettes : l’information perdue a un coût. Placez le départ avec θ_A < θ_B : EM trouve le même sommet avec les noms des pièces échangés — un premier signe que le résultat dépend du point de départ.

Pourquoi ça monte toujours : une courbe qui minore et qui touche

L’exemple de l’article de 1977 vient de la génétique (données de Rao) : 197 animaux répartis en quatre catégories, d’effectifs 125, 18, 20 et 34, de probabilités ½ + θ/4, (1 − θ)/4, (1 − θ)/4 et θ/4. La première case mélange deux sous-cases invisibles, de probabilités ½ et θ/4. Étape E : on répartit les 125 animaux entre elles au prorata, soit 125·θₜ/(2 + θₜ) dans la sous-case θ/4. Étape M : θₜ₊₁ est alors une simple fréquence. Géométriquement, l’étape E construit une courbe (orange) qui reste sous la log-vraisemblance et la touche en θₜ ; l’étape M saute au sommet de cette courbe. Comme la vraie courbe est au-dessus, elle a forcément monté.

À gauche : log-vraisemblance ℓ(θ) (bleu), courbe minorante ℓ(θₜ) + Q(θ | θₜ) − Q(θₜ | θₜ) (orange), point de contact θₜ et sommet de la minorante θₜ₊₁ (rouge). À droite : distance au maximum |θₜ − θ̂| en échelle logarithmique, itération par itération.

Départ θ₀

Itération, θₜ et ℓ(θₜ)—

Maximum exact—racine positive de 197θ² − 15θ − 68 = 0, obtenue en annulant la dérivée de ℓ

La minorante reste sous ℓ—

Rapport des erreurs successives—

L’écart entre la log-vraisemblance et sa minorante n’est pas quelconque : c’est une divergence de Kullback-Leibler entre deux lois de la donnée cachée, toujours positive ou nulle (inégalité de Jensen) et nulle au point de contact. Ce seul fait garantit la montée, quel que soit le modèle. Et l’erreur est multipliée à chaque tour par un facteur presque constant, environ 0,133 ici : c’est la convergence linéaire, dont la planche 3 donne la raison.

La vitesse : la part d’information manquante

Dempster, Laird et Rubin ont montré que, près du sommet, l’erreur est multipliée à chaque tour par la fraction d’information manquante : 1 − I_obs/I_comp, où I_obs mesure la courbure de la vraie log-vraisemblance et I_comp celle qu’on aurait avec les données complètes. Expérience : 500 mesures viennent de deux groupes de lois normales N(−d, 1) et N(d, 1), et l’on cherche seulement la proportion w du premier groupe (vraie valeur 0,3). Quand les groupes se chevauchent, l’appartenance de chaque mesure est presque inconnue : l’information manque, et EM rampe. On le compare à la méthode de Newton, qui utilise la dérivée seconde.

À gauche : histogramme des 500 mesures et les deux groupes estimés (bleu et orange). À droite : erreur |wₜ − ŵ| en échelle logarithmique, pour EM (bleu) et pour Newton (orange) ; la droite pointillée a la pente prédite par la théorie.

Écart entre les groupes d

Proportion estimée ŵ—

Facteur de réduction de l’erreur, mesuré / théorique—

Itérations pour une erreur 10⁻¹⁰ : EM / Newton—

Mesuré et prédit coïncident à quatre décimales. Avec d = 3, les groupes sont bien séparés, presque rien ne manque et EM converge en quelques tours ; avec d = 0,3, il manque plus de 90 % de l’information et il faut près de 300 tours. Newton, lui, converge quadratiquement en quelques pas, mais chaque pas exige des dérivées secondes et peut sortir de l’intervalle permis. EM ne sort jamais de [0 ; 1] : chaque étape M donne une vraie proportion.

Trois usages : nuages de points, groupes sanguins, valeurs manquantes

Trois groupes de points dans le plan, chacun suivant une loi normale à deux dimensions, avec sa propre forme : on ignore à quel groupe appartient chaque point. EM pour un mélange gaussien attribue à chaque point une responsabilité envers chaque groupe (étape E, couleur mélangée), puis recalcule pour chaque groupe son poids, son centre et son ellipse de dispersion (étape M). L’algorithme des k-means (k moyennes) en est la version « dure » : chaque point est attribué entièrement au centre le plus proche, et tous les groupes sont supposés ronds et de même taille.

À gauche : EM, couleur de chaque point = mélange de ses trois responsabilités, ellipses à deux écarts-types. À droite : k-means partant des mêmes centres, avec ses frontières de Voronoï implicites.

Allongement des groupes

Itération EM et log-vraisemblance—

Points mal classés : EM / k-means—

Le groupe sanguin dépend de deux allèles parmi A, B et O ; O est récessif. On voit le phénotype (A, B, AB ou O), pas le génotype : une personne du groupe A peut être AA ou AO. On veut les fréquences p, q, r des trois allèles. Étape E : sous l’équilibre de Hardy-Weinberg, les personnes du groupe A se répartissent entre AA et AO dans le rapport p² : 2pr, et de même pour B. Étape M : on compte les allèles dans ces génotypes attendus. C’est le comptage de gènes de Cedric Smith (1955), un EM avant la lettre. Échantillon fictif de 1 000 donneurs, aux proportions voisines de celles observées en France.

À gauche : les quatre phénotypes et, pour A et B, la part attendue d’homozygotes (foncé) et d’hétérozygotes (clair). À droite : p, q et r au fil des itérations (bleu, orange, vert).

Effectif du groupe AB

Fréquences p ; q ; r par EM—

Contrôle : maximisation directe—

Formules de Bernstein (1925)—

On relève deux notes par élève, en mathématiques (x) et en physique (y). Les élèves faibles en mathématiques ont souvent évité l’épreuve de physique : leur note y manque. La moyenne des seules notes présentes surestime alors la moyenne de physique. EM traite y comme une donnée manquante : étape E, remplacer chaque note absente par sa prédiction sachant x, avec la droite de régression actuelle (et ajouter la variance résiduelle dans les sommes de carrés) ; étape M, recalculer moyennes, variances et covariance. Le manque dépend de x, qui est observé : c’est le cas manquant au hasard (MAR) où la vraisemblance reste valable. Comme seule y peut manquer (motif « monotone »), une seule itération atteint déjà la formule explicite d’Anderson ; EM devient indispensable quand les deux notes peuvent manquer, ou pour des modèles plus riches.

Points bleus : élèves avec les deux notes. Points orange : notes de physique manquantes, placées à leur valeur prédite par l’étape E en cours. Traits horizontaux : vraie moyenne (vert), moyenne des notes présentes (rouge), estimation EM (bleu).

Part des notes manquantes

Vraie moyenne de y (population)12valeur du modèle de simulation

Moyenne des notes présentes—

Estimation EM—

Contrôle : formule d’Anderson (1957)—

Dans les trois cas, l’étape M est la statistique la plus banale (moyennes, comptages, covariances) appliquée à des données complétées par l’étape E. Les nuages allongés montrent la différence entre une attribution souple qui apprend la forme des groupes et une attribution dure qui suppose des groupes ronds. Les groupes sanguins montrent qu’EM respecte automatiquement p + q + r = 1, ce que les formules de Bernstein ne font pas. Les notes manquantes montrent qu’il corrige le biais de l’analyse des seuls cas complets, et qu’il retombe exactement sur la formule explicite connue dans ce cas.

Signature — le bon sommet dépend du départ

EM monte toujours, mais il monte vers le sommet le plus proche, qui n’est pas forcément le plus haut. La log-vraisemblance d’un mélange a souvent plusieurs sommets ; chacun attire une région de points de départ, son bassin d’attraction.

200 mesures viennent d’un mélange : 30 % autour de μ₁ = 0 et 70 % autour de μ₂ = 3 (écarts-types 1, proportions connues). On ne cherche que les deux moyennes. Chaque point du plan (μ₁, μ₂) est un départ possible : on y lance EM et l’on colore le point selon le sommet atteint. Touchez ou cliquez la carte pour lancer EM depuis le point choisi et suivre son chemin.

Lignes de niveau de la log-vraisemblance. Fond bleu : départs qui mènent au maximum global ; fond orange : départs qui mènent au maximum secondaire, où les deux groupes sont intervertis ; la diagonale μ₁ = μ₂ (gris) mène à un point-selle. En rouge : le chemin d’EM depuis votre départ.

Maximum global—

Maximum secondaire—

Point-selle de la diagonale—

Votre départtouchez la carte

Parts des départs de la grille—

Un mélange de cinq lois normales sur une droite (poids, moyennes et écarts-types tous inconnus, 250 mesures). On lance EM 120 fois, chaque fois depuis cinq mesures tirées au hasard comme moyennes de départ. Chaque essai finit sur un sommet ; si une fraction p des essais atteint le meilleur, la probabilité de le trouver au moins une fois en k essais vaut 1 − (1 − p)ᵏ.

À gauche : histogramme des mesures, meilleur ajustement trouvé (bleu) et ajustement le plus fréquent parmi les autres (orange). À droite : log-vraisemblance finale des essais, rangés par ordre croissant ; la marche verte la plus haute est le meilleur sommet.

Essais terminés—

Part p des essais qui trouvent le meilleur sommet—

Essais nécessaires pour le trouver à 99 %—

Contrôle : tirages parmi les essais réels, k = 5—

La carte est coupée exactement en deux par la diagonale : un départ « du mauvais côté » (μ₁ > μ₂) mène toujours au sommet où les groupes sont intervertis, de log-vraisemblance nettement plus basse. Ce n’est pas un hasard : ici, l’étape M ne change jamais l’ordre des deux moyennes (la responsabilité du premier groupe décroît quand la mesure grandit, donc sa moyenne pondérée reste la plus petite), si bien qu’EM ne peut pas franchir la diagonale. Rien dans l’algorithme ne le signale : la montée est garantie, le sommet ne l’est pas. D’où les deux réflexes des praticiens : partir d’un point raisonnable (par exemple le résultat de k-means) et multiplier les départs en gardant le meilleur. La garantie de Wu (1983) est exactement celle-ci : convergence vers un point stationnaire, pas vers le maximum global.

Repères historiques

1886Simon Newcomb, astronome, estime une mesure en supposant les observations issues d’un mélange de lois normales, par des repondérations successives : une première forme du principe.
1926Anderson McKendrick applique une idée semblable à des données d’épidémie dont une partie n’est pas observable.
1955Ceppellini, Siniscalco et Smith publient le comptage de gènes pour les groupes sanguins (planche 4).
1958Herman Hartley propose une méthode générale d’estimation à partir de données incomplètes de comptage.
1970Leonard Baum et ses collègues publient l’algorithme d’apprentissage des chaînes de Markov cachées (Baum-Welch), un EM pour la reconnaissance de la parole.
1974Rolf Sundberg étudie l’algorithme pour les familles exponentielles.
1977Arthur Dempster, Nan Laird et Donald Rubin unifient tous ces cas, nomment l’algorithme EM et prouvent la montée, dans le Journal of the Royal Statistical Society (série B).
1983C. F. Jeff Wu corrige la preuve de convergence et en précise les conditions.
2008Chuong Do et Serafim Batzoglou popularisent l’exemple des deux pièces (planche 1).

Ce qu’il faut retenir

Le principe

Compléter les données par leur valeur attendue (E), puis estimer comme si elles étaient complètes (M). Chaque étape est simple ; l’alternance résout un problème difficile.

La garantie

La log-vraisemblance ne baisse jamais : l’étape E construit une minorante qui touche la courbe, l’étape M en prend le sommet. Les contraintes (proportions, variances positives) sont respectées d’elles-mêmes.

Les limites

Convergence lente, au rythme de la fraction d’information manquante ; arrivée possible sur un sommet secondaire ou un point-selle. Remèdes : bon départ, départs multiples, accélérations.

Pour aller plus loin

La décomposition clé. Si z désigne ce qui manque, ℓ(θ) = Q(θ | θₜ) − H(θ | θₜ), avec H(θ | θₜ) = E[ln k(z | y, θ)], l’espérance étant prise sous la loi de z sachant y et θₜ. L’inégalité de Jensen donne H(θ | θₜ) ≤ H(θₜ | θₜ) : augmenter Q suffit pour augmenter ℓ. C’est pourquoi un EM « généralisé », qui ne fait qu’augmenter Q sans la maximiser, garde la propriété de montée.

Lien avec k-means. Quand les variances des groupes tendent vers zéro, les responsabilités deviennent 0 ou 1 et EM pour un mélange gaussien de groupes ronds devient exactement l’algorithme de Lloyd, alias k-means.

Variantes. ECM (maximisations partielles successives), SEM (EM supplémenté, qui récupère les erreurs-types à partir de la vitesse de convergence), EM stochastique et Monte-Carlo (étape E simulée), EM variationnel, accélérations par Aitken ou par extrapolation. La maximisation a posteriori (MAP) se traite de la même façon en ajoutant le logarithme de la loi a priori à Q.