[java] algorithme Expectation-maximization (EM)

[pseudocode]

Une implémentation Java de l'algorithme Expectation-maximization (EM).


exemple de convergence avec 3 distributions normales (voir post suivant)

L'interface qui définit une loi de probabilité

Code java :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
interface Law {
	/**
	 * @param x some value
	 * @return the probability of the value x
	 */
	double proba(double x);
 
	/**
	 * improve law parameters
	 * 
	 * @param N number of samples
	 * @param x samples
	 * @param tk probability of each sample
	 */
	void improveParameters(int N, double[] x, double[] tk);
}

Le code de l'algorithme EM

Code java :

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
48
49
50
51
52
53
54
/**
 * Compute the mixture coefficients using EM algorithm
 * 
 * @param x sample values
 * @param laws instances of the laws
 * @return mixture coefficients 
 */
public double[] algorithmEM(double[] x, Law[] laws) {
	int N=x.length;
	int G=laws.length;
 
	double[] pi = new double[G];
	double[][] t = new double[G][N];
 
	// initial guess for the mixture coefficients (uniform)
	for(int k=0;k<G;k++) pi[k]=1.0/G;
 
	// iterative loop (until convergence or 5000 iterations)
	double convergence;
	for(int loop=0;loop<5000;loop++) {
		convergence=0;
 
		// ---- 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[k][i]=numerator/denominator;
			}
		}
 
		// ---- M Step ----
 
		// mixture coefficients (maximum likelihood estimate of binomial distribution)
		for(int k=0;k<G;k++) {
			double savedpi=pi[k];
			pi[k]=0;
			for(int i=0;i<N;i++) pi[k]+=t[k][i];
			pi[k]/=N;
			convergence+=(savedpi-pi[k])*(savedpi-pi[k]);
		}
 
		// update the parameters of the laws
		for(int k=0;k<G;k++)
			laws[k].improveParameters(N, x, t[k]);
 
		if( convergence < 1E-10 ) break;
	}
 
	return pi;
}

[pseudocode]

Un exemple d'application avec des distributions normales (loi de probabilité gaussiennes).



Définition d'une loi pour une distribution normale

Code java :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
public class NormalDistribution implements Law {
	private double mean=0, sigma=0;
	public NormalDistribution(double mean, double sigma) {
		this.mean=mean; this.sigma=sigma;
	}
	public double proba(double x) {
		double t = (x-mean);
		return Math.exp(-(t*t)/(2*sigma*sigma)) / (sigma*Math.sqrt(2*Math.PI));
	}	
	public void improveParameters(int N, double[] x, double[] tk) {
		double sumTkX=0,sumTk=0;
		for(int i=0;i<N;i++) sumTkX+=tk[i]*x[i];
		for(int i=0;i<N;i++) sumTk+=tk[i];
		mean = sumTkX/sumTk;
		double sumTkXc2=0;
		for(int i=0;i<N;i++) sumTkXc2+=tk[i]*(x[i]-mean)*(x[i]-mean);
		sigma = Math.sqrt(sumTkXc2/sumTk);
	}
	@Override public String toString() {
		return "Gaussian( mean="+mean+" , sigma="+sigma+" )";
	}
}

Génération aléatoire de valeurs suivant 3 lois gaussiennes

Code java :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
static Random random = new Random(1);
 
public static double gaussianSample(double mu, double sigma) {
	return (random.nextGaussian())*sigma + mu;
}
 
public double[] generateSample(int N) {
	double[] x = new double[N];
 
	// generate random values according to some laws
	for(int i=0;i<N;i++) {
		double r = (100*i)/N;
		if (r<20) // 20% of law #1
			x[i]=gaussianSample(2.0,0.5); 
		else if (r<70) // 50% of law #2
			x[i]=gaussianSample(6.0,1.0); 
		else // 30% of law #3
			x[i]=gaussianSample(10.0,2.0); 
	}
	return x;
}

Et enfin l'appel de l'algo EM:

Code java :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
// initial guess for the laws parameters (uniform)
int G=3;
Law[] laws = new Law[G];
laws[0] = new NormalDistribution(4.0,1.0);
laws[1] = new NormalDistribution(8.0,1.0);
laws[2] = new NormalDistribution(12.0,1.0);
 
// performe EM algorithm
double pi[] = algorithmEM(x, laws);
 
// display mixture coefficients and law parameters
for(int k=0;k<G;k++)
	System.out.printf("%f * %s\n",pi[k],laws[k]);
 
// display:
// 0,200126 * Gaussian( mean=2.0079143733645988 , sigma=0.4937122623373714 )
// 0,498746 * Gaussian( mean=6.016170280615971 , sigma=0.9894557053479284 )
// 0,301128 * Gaussian( mean=9.982667768339546 , sigma=2.0338918908916916 )

[pseudocode]

Envoyé par daurey68

Bonjour,

je suis en stage pour mon projet de fin d'étude. Les données que je dois traiter sont des courbes dont l'équation est inconue et qui doivent être décomposer en une combinaison linéaire de gaussienne.

J'ai tout d'abord essayé de coder moi-même la méthode d'expectation-maximization mais elle ne fonctionne pas.

L'algo EM sert à déterminer les paramètres d'une loi de probabilité connaissant des valeurs prises par une variable aléatoire. L'histogramme gris représenté en fond d'image n'est pas la "valeur" des données d'entrée mais leur distribution.

L'algo EM ne sert pas à faire de l'ajustement de courbe (trouver les parametres d'une courbe en connaissant des "points de passages"). Pour cela, il faut plutot voir des algorithmes de http://en.wikipedia.org/wiki/Curve_fitting

[daurey68]

Merci pour cette précision sur la distribution des probabilités.
L'erreur dans le graphique doit donc provenir du fait que j'ai affiché les valeurs obtenues et non la probabilité d'obtenir une certain valeur.

Cependant, cette méthode est celle que je dois utiliser pour traiter ma courbe :
"To estimate the numbers and characterics of the echoes in the waveforms, Gaussians distribution functions are fitted to the siganl. The waveforms are modelled as a sum of Gaussian distributions." Extrait d'un article scientifique

Donc si j'ai bien compris, votre programme permet :
- de générer un certain nombre de points dont les valeurs sont aléatoires
- de créer l'histogramme de fréquence de ces valeurs
- d'appliquer à cet histogramme l'algo EM afin de lui faire correspondre des gaussiennes.

Dans quelle partie de votre code, l'histogramme est-il créé ?

Cela me permettrait de modifier mes données initiales afin de les normaliser pour avoir des valeurs entre 0 et 1 et de n'utiliser que la partie EM de votre programme.

[pseudocode]

Envoyé par daurey68

Donc si j'ai bien compris, votre programme permet :
- de générer un certain nombre de points dont les valeurs sont aléatoires
- de créer l'histogramme de fréquence de ces valeurs
- d'appliquer à cet histogramme l'algo EM afin de lui faire correspondre des gaussiennes.

Dans quelle partie de votre code, l'histogramme est-il créé ?

Non. L'algo EM ne travaille pas sur les histogrammes. Il travaille directement sur les valeurs. C'est là toute sa force !

Il n'y a donc pas besoin de calculer l'histogramme dans le code. Et heureusement, car une variable gaussienne peut prendre une infinité de valeurs (entre -infini et + infini), ce qui obligerait à créer un histogramme de taille infinie !

J'ai dessiné l''histogramme ainsi que les gaussiennes uniquement pour "montrer" l'évolution de l'algorithme. Mais ce n'est pas ainsi que l'algo travaille. Il ne cherche pas à "ajuster" les gaussiennes sur l'histogramme.

[daurey68]

Ok, l'histogramme n'est pas créé.
J'ai transformé des données de bases pour les faire correspondre à un ensemble de variable aléatoire (correspondant à votre x= randon(nbvaleurs) ).
Le programme fonctionne apparemment correctement.

Il me reste un seul problème : le nombre de gaussienne déterminé n'est pas celui que je souhaiterai.
Le fichier en PJ illustre ce problème. L'histogramme pourrait être assimilé à la courbe originelle en jaune. Le nombre de pics à détecter est de 2 et le nombre de gaussienne ajustées doit donc aussi être de 2. En pratique seule une seule gaussienne est véritablement déterminée (les autres ont un poids trop faible).

Savez-vous s'il est possible d'introduire un seuil dans l'algorithme EM afin de résoudre mon problème ?

[pseudocode]

Envoyé par daurey68

Il me reste un seul problème : le nombre de gaussienne déterminé n'est pas celui que je souhaiterai.
Le fichier en PJ illustre ce problème. L'histogramme pourrait être assimilé à la courbe originelle en jaune. Le nombre de pics à détecter est de 2 et le nombre de gaussienne ajustées doit donc aussi être de 2. En pratique seule une seule gaussienne est véritablement déterminée (les autres ont un poids trop faible).

Savez-vous s'il est possible d'introduire un seuil dans l'algorithme EM afin de résoudre mon problème ?

L'algo EM est assez sensible aux conditions initiales (valeurs initiales des lois).

Si vous savez au départ le nombre de gaussienne (2), le plus simple est de faire une première analyse rapide de l'histogramme pour déterminer les 2 plus hauts sommets, et de centrer les gaussiennes sur ces 2 sommets.

[daurey68]

j'avais déjà esasyé de centrer les gaussiennes sur les valeurs que je devais obtenir.
Mais cela ne change rien sur le fait qu'une seule gaussienne est calculée.

[daurey68]

Je pense que mes pics sont trop rapprochés pour pouvoir être détectés comme 2 maximums de gaussiennes différents.
Cela est-il possible ?

J'ai essayé d'étaler mon histogramme sur un plus grand intervalle pour espacer les pics. Cela ne donne pas des résultats concluants.

Y aurait-il un autre moyen de contourner ce problème ?

[pseudocode]

Envoyé par daurey68

Y aurait-il un autre moyen de contourner ce problème ?

Heu... arrêter de prendre l'algo EM pour faire ce genre de chose. Mais je l'ai déjà dit je crois.

[hsalam]

Bonjour,
j'ai un question concernant le data d'entrée. Si le data est une image alors l'entrée est des samples des valeurs des niveaux de gris de l'image??? Merci beaucoup en avance

[pseudocode]

Envoyé par hsalam

Bonjour,
j'ai un question concernant le data d'entrée. Si le data est une image alors l'entrée est des samples des valeurs des niveaux de gris de l'image??? Merci beaucoup en avance

oui, tu peux utiliser les niveaux de gris des pixels d'une image comme entrée. L'algo EM te donnera la "meilleure" modélisation de l'histogramme

[hsalam]

Est-que vous peuvez me donne un exemple de comment initialiser par des samples des niveaux de gris de l'image. Est-que je met tous les valeurs des niveaux de gris dans un vecteur et après faire un tirage de ce vecteur, ou je peut utiliser tous les valeurs? le taille de l'image est grand et si je veux met tous les valeurs ça va coût cher par rapport au temps de calcul! quest-que vous me consiel?
Est-que vous avez par chance un programme en language C qui fait la même comme le votre. Merci beaucoup en avance.

[pseudocode]

Envoyé par hsalam

Est-que je met tous les valeurs des niveaux de gris dans un vecteur et après faire un tirage de ce vecteur, ou je peut utiliser tous les valeurs? le taille de l'image est grand et si je veux met tous les valeurs ça va coût cher par rapport au temps de calcul!

Tout dépend de la taille de l'image et de la puissance du processeur.
Il faut faire des tests, je ne peux pas répondre a cette question/

Est-que vous avez par chance un programme en language C qui fait la même comme le votre. Merci beaucoup en avance.

Non, je ne l'ai pas porté en C. Le portage de la partie calcul ( méthode algorithmEM() ) ne pose pas de problème, c'est pratiquement du C. Pour le reste, il faut remplacer les interfaces/classes par des structures.

[hsalam]

qu'est que votre programme prend comme entrée? est-qu'il prend la probabilité des lois comme entrée?
normalement l'alg EM estime le variance et le mean des distributions gaussiennes à partir des valeurs ou aussi il faut met les probabilités comme input?
Merci

[pseudocode]

Envoyé par hsalam

qu'est que votre programme prend comme entrée? est-qu'il prend la probabilité des lois comme entrée?
normalement l'alg EM estime le variance et le mean des distributions gaussiennes à partir des valeurs ou aussi il faut met les probabilités comme input?
Merci

Voila la signature de la méthode de l'algorithme

Code java :

1
2
3
4
5
6
7
8
/**
 * Compute the mixture coefficients using EM algorithm
 * 
 * @param x sample values
 * @param laws instances of the laws
 * @return mixture coefficients 
 */
public double[] algorithmEM(double[] x, Law[] laws) {

En entrée :

double[] x = liste de "valeurs" (réalisations d'une variable aléatoire)

En entrée/sortie :

Law[] laws = liste de "Loi de probabilité" (par exemple, 3 gaussiennes)
--> en entrée, des valeurs de paramètres quelconques (= "initial guess").
--> en sortie, les valeurs de paramètres estimées par l'algo

En sortie :

double[] : liste des coefficients multiplicateurs de chaque loi

[ellgafsi]

salut tres belle contribution!

Malgres ca me dit rien, layout, l implementation et les explications sont super!
ca sert a quoi cet algo?

[pseudocode]

Envoyé par ellgafsi

salut tres belle contribution!

Malgres ca me dit rien, layout, l implementation et les explications sont super!
ca sert a quoi cet algo?

Comme le dit très bien wikipedia:

L'algorithme espérance-maximisation (en anglais Expectation-maximisation algorithm, souvent abrégé EM), proposé par Dempster et al. (1977), est une classe d'algorithmes qui permettent de trouver le maximum de vraisemblance des paramètres de modèles probabilistes lorsque le modèle dépend de variables latentes non observables.

Ca sert souvent à faire de la classification. Des données d'entrée sont utilisées comme exemple d'apprentissage, l'algo EM construit un modèle probabiliste basé sur cet exemple, et au final on peut classifier de nouvelles données en se basant sur la probabilité d'appartenance a chaque loi.

[hsalam]

Bonjour,
Je suis nouveau à java. j'ai compilé votre programme, mais je reçois cette erreur: class, interface, or enum expected dans 56 lignes. est-que vous pouvez m'aider ? je ne sais pas est-que ça veut dire cette erreur?
Est-que je besoin de faire un fonction main ou votre programme est executable comme il est?
Merci beacoup en avance

[pseudocode]

Envoyé par hsalam

Bonjour,
Je suis nouveau à java. j'ai compilé votre programme, mais je reçois cette erreur: class, interface, or enum expected dans 56 lignes. est-que vous pouvez m'aider ? je ne sais pas est-que ça veut dire cette erreur?
Est-que je besoin de faire un fonction main ou votre programme est executable comme il est?
Merci beacoup en avance

Non, le code n'est pas compilable tel quel. Il s'agit juste du de la definition de l'interface et du corps de la méthode implémentant l'algorithme.

Le code source complet de l'exemple est disponible dans cette archive. Il suffit de le compiler et d'executer la classe EM.Test