Algorithme Expectation Maximization

**yarf** · 30/04/2009, 20h06

Bonjour,

je suis à la recherche d'une implémentation simple de l'algorithme EM, de préférence en Java, capable de s'appliquer au problème suivant :

Pour une ressource r donnée, sa probabilité est calculée selon trois modèles en utilisant une interpolation linéaire :

P(r) = A P1(r) + B P2(r) + (1 - A - B) P3(r)

Il faudrait donc que les deux paramètres A et B soient calculés selon l'algorithme EM.

Auriez-vous une suggestion ?

**FR119492** · 30/04/2009, 23h09

Salut!
Tu as une explication très détaillée dans Wikipedia.
Jean-Marc Blanc

**yarf** · 01/05/2009, 09h03

Merci, mais j'avais déjà passé une heure à essayer de comprendre l'explication donnée par wikipédia. Cela m'a paru très complexe, et j'ai pensé me tourner vers une solution toute faite.

**pseudocode** · 01/05/2009, 12h09

Tes 3 lois de proba P1,P2,P3 n'ont pas de paramètres indéterminés ?

**yarf** · 01/05/2009, 14h36

Non, seuls A et B le sont.

**pseudocode** · 01/05/2009, 15h19

Envoyé par yarf

Non, seuls A et B le sont.

Dans ce cas tu as juste besoin de calculer les ratios du mélange, comme indiqué sur wikipedia:

- l'étape E se réduit a la formule de Bayes

- l'étape M au calcul des ratios.

Tu boucles sur ces 2 calculs jusqu'a convergence, et tu obtiens les ratios du mélange (pi1,pi2,pi3)

**yarf** · 01/05/2009, 17h33

Merci pseudocode. Cela a l'air plus simple je vais essayer de m'y coller.

Sinon, histoire de ne pas réinventer la roue, y a-t-il des implémentations de cela que je pourrais utiliser ?

**pseudocode** · 01/05/2009, 17h51

Envoyé par yarf

Merci pseudocode. Cela a l'air plus simple je vais essayer de m'y coller.

J'ai testé les formules de wikipedia et l'algo en lui même prend une dizaine de lignes de code. En fait, ce qui m'a pris le le plus de temps/place c'est le code pour créer le test unitaire.

Sinon, histoire de ne pas réinventer la roue, y a-t-il des implémentations de cela que je pourrais utiliser ?

Hum... je ne sais pas si ca existe. Il n'y a pas beaucoup de valeur ajoutée à avoir l'algo EM déja codé. Ce qui est "compliqué" c'est la partie de l'étape M dans laquelle on doit trouver les paramètres des lois. Et cette partie dépend uniquement des lois et pas de l'algo EM.

**yarf** · 01/05/2009, 18h02

OK, merci beaucoup. Bonne fin de journée !

**yarf** · 01/05/2009, 18h58

Je suis encore en train de réfléchir dessus, mais comme je suis plutôt pressé, je me permets de poser encore quelques questions...

J'utilise les P(r) calculées pour trier les ressources, et ensuite je mesure la liste obtenue selon une mesure particulière. Puis-je/dois-je inclure cette mesure dans l'algorithme ? Il me semble que oui, le but final étant de maximiser le résultat obtenu selon cette mesure. Aussi, puis-je/dois-je remplacer le maximum de vraisemblance par cette mesure ?

Je continue de réfléchir...

**Matthieu Brucher** · 01/05/2009, 19h11

Comme le dit pseudocode, il n'y a pas d'algorithme de référence. Ca serait juste :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
loop
{
  Expectation()
  Maximization()
}while (critere)

Si tu veux injecter ta mesure dans l'EM, c'est que tu as une loi de probabilité à maximiser différente. Si tu l'écrit proprement, tu sauras comment l'injecter dans ton algorithme.

**pseudocode** · 01/05/2009, 19h23

Envoyé par yarf

J'utilise les P(r) calculées pour trier les ressources, et ensuite je mesure la liste obtenue selon une mesure particulière. Puis-je/dois-je inclure cette mesure dans l'algorithme ? Il me semble que oui, le but final étant de maximiser le résultat obtenu selon cette mesure. Aussi, puis-je/dois-je remplacer le maximum de vraisemblance par cette mesure ?

L'algo EM utilise la mesure du "maximum de vraisemblance". Si tu as une autre mesure, rien ne dit que l'algo EM soit facilement implémentable, ni qu'il converge.

**yarf** · 01/05/2009, 20h19

Merci pour vos réponses !

**yarf** · 03/05/2009, 19h41

Bonjour,

après quelques heures supplémentaires passées sur le problème, je me suis rendu compte que l'utilisation de la vraisemblance était bien appropriée et j'en suis venu à me diriger vers ce qui semble être un cas particulier de l'algorithme EM, et qui s'appelle Deleted Interpolation.

Si j'ai bien compris, l'algorithme est basé sur le fait que dans le cas d'une combinaison linéaire de deux sous-modèles :

P(r) = A P1(r) + (1 - A) P2(r)

on peut considérer cette combinaison comme un polynôme d'ordre 1 en A, et donc la vraisemblance comme un polynôme d'ordre N (le nombre d'éléments de l'échantillon). Il est montré que ce polynôme a soit un maximum pour A = 0 et la vraisemblance décroit de façon monotone, soit pour A = 1 et la vraisemblance croit de façon monotone, soit entre les deux, de sorte que la vraisemblance croisse entre 0 et le maximum, et décroisse entre le maximum et 1.

L'algorithme consiste donc à étudier la dérivée de la vraisemblance, afin de déterminer dans lequel de ces trois cas on se trouve. Dans le premier cas A = 0, dans le second A = 1. Dans le troisième cas, la valeur optimale est obtenue par encadrements successifs.

Ce principe est transposable pour plus de deux sous-modèles. Je donne cela à titre informatif, et aussi pour me fournir une assurance supplémentaire de ma bonne compréhension. Si vous avez des remarques, n'hésitez pas à m'en faire part.

Il me semble qu'il faut diviser le corpus en un corpus d'apprentissage pour les sous-modèles, et un corpus de validation pour l'algorithme. Mais avec cet algorithme, il me semble qu'il n'y a pas d'ajustement des paramètres sur plusieurs échantillons, comme il en est question dans l'algorithme EM. Du coup, j'ai l'impression que tout ce que fait l'algorithme, c'est des encadrements successifs ! Je ne connais pas encore bien tout cela, qu'en pensez-vous ?

Je vous donnerais bien l'article dont cela provient, mais cela est interdit...

**Matthieu Brucher** · 03/05/2009, 19h44

Envoyé par yarf

Je vous donnerais bien l'article dont cela provient, mais cela est interdit...

Donne la référence tout de même, pour ceux qui y auraient accès.

**pseudocode** · 03/05/2009, 19h57

Envoyé par yarf

L'algorithme consiste donc à étudier la dérivée de la vraisemblance, afin de déterminer dans lequel de ces trois cas on se trouve. Dans le premier cas A = 0, dans le second A = 1. Dans le troisième cas, la valeur optimale est obtenue par encadrements successifs.

Est-ce plus simple/rapide que d'utiliser directement l'algo EM ?

**yarf** · 03/05/2009, 20h13

La voici :

A Fast Algorithm for Deleted Interpolation
L.R. Bahl, P.F. Brown, P.V. de Souza, R.L. Mercer, D. Nahamoo
June 1991
TechReport (Microsoft)

http://research.microsoft.com/apps/p....aspx?id=64740

**yarf** · 03/05/2009, 20h20

Envoyé par pseudocode

Est-ce plus simple/rapide que d'utiliser directement l'algo EM ?

Je suppose que oui, mais je ne sais pas exactement, je suis encore en train de découvrir cela. Il paraît surtout plus simple.

Au passage, l'algorithme ne s'appelle pas Deleted Interpolation, il permet de calculer les coefficients d'une deleted interpolation. Il n'a pas de nom.

**pseudocode** · 03/05/2009, 20h51

Envoyé par yarf

Je suppose que oui, mais je ne sais pas exactement, je suis encore en train de découvrir cela. Il paraît surtout plus simple.

Pourtant la boucle EM dans ton cas est très simple.

Code java :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
 
interface Law {
	double proba(double x);
}
 
double[] mixtureCoefficients(double[] x, Law[] laws) {
	int N=x.length, G=laws.length;
 
	double[] pi = new double[G];
	double[][] t = new double[N][G];
 
	// initial guess for mixture coefficients (uniform)
	for(int k=0;k<G;k++) pi[k]=1.0/G;
 
	// iterative loop (TODO: until convergence of pi)
	for(int loop=0;loop<20;loop++) {
 
		// E Step (Bayes inversion formula)
		for(int i=0;i<N;i++) {
			double denominator = 0;
			for(int l=0;l<G;l++) denominator+=pi[l]*laws[l].proba(x[i]);
			for(int k=0;k<G;k++) {
				double numerator = pi[k]*laws[k].proba(x[i]);
				t[i][k]=numerator/denominator;
			}
		}
 
		// M Step (maximum likelihood estimate of binomial distribution)
		for(int k=0;k<G;k++) {
			pi[k]=0;
			for(int i=0;i<N;i++) pi[k]+=t[i][k];
			pi[k]/=N;
		}
	}
	return pi;
}

J'ai testé avec 3 distributions normales:

Code java :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
 
class NormalDistribution implements Law {
	double mu,sigma,normalizer;
	public NormalDistribution(double mu, double sigma) {
		this.mu=mu; this.sigma=sigma;
		this.normalizer = sigma*Math.sqrt(2*Math.PI);
	}
	public double proba(double x) {
		double t = (x-mu);
		return Math.exp(-(t*t)/(2*sigma*sigma)) / normalizer;
	}	
}
 
static Random random = new Random(0);
 
static double gaussianSample(double mu, double sigma) {
	return (random.nextGaussian())*sigma + mu;
}
 
public static void main(String[] args) {
	int N=1000;
	double[] x = new double[N];
 
	// generate random values according to the laws
	for(int i=0;i<N;i++) {
		double r = (double)i/N;
		// 50% of laws[0], 35% of laws[1], 15% of laws[2]
		int l =  (r<0.50)?0:(r<0.85)?1:2;
		switch(l) {
		case 0 : x[i]=gaussianSample(3.0, 1.0);  break;
		case 1 : x[i]=gaussianSample(5.0, 2.0);  break;
		case 2 : x[i]=gaussianSample(10.0, 3.0); break;
		} 
	}
 
	// Compute mixture coefs for the 3 following distributions
	Law[] laws = new Law[3]; 
	laws[0] = new NormalDistribution(3.0,1.0);
	laws[1] = new NormalDistribution(5.0,2.0);
	laws[2] = new NormalDistribution(10.0,3.0);
 
	double[] coefs = mixtureCoefficients(x,laws);
 
	System.out.println(Arrays.toString(coefs));
	// result is [0.5034928298008519, 0.33363317002051407, 0.16287400017863374]
 
}