[java] algorithme Expectation-maximization (EM)

Version imprimable

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

http://xphilipp.developpez.com/contribuez/algoEM.gif
exemple de convergence avec 3 distributions normales (voir post suivant)

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

Code:

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:

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;
}

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

http://xphilipp.developpez.com/contribuez/algoEM.png

Définition d'une loi pour une distribution normale

Code:

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:

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:

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 )

19/05/2010, 11h37
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
20/05/2010, 00h10
pseudocode

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
20/05/2010, 11h44
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.
20/05/2010, 12h21
pseudocode

Citation:

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. :D
Il faut faire des tests, je ne peux pas répondre a cette question/

Citation:

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.
21/05/2010, 15h26
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

Citation:

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:

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

25/05/2010, 11h59
ellgafsi

salut tres belle contribution!

Malgres ca me dit rien, layout, l implementation et les explications sont super! :) :ccool:
ca sert a quoi cet algo?
25/05/2010, 12h26
pseudocode

Comme le dit très bien wikipedia:

Citation:

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.
25/05/2010, 14h09
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 :)
25/05/2010, 16h16
pseudocode

Non, le code n'est pas compilable tel quel. Il s'agit juste du de la définition 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'exécuter la classe EM.Test

Salut,
j'ai fait ce code pour prendre le pixels d'image comme entrée , mais j'ai l'erreur :

ImageTest.java:30: '.class' expected
return d[];
          ^

Comme je suis nouveau à java j'ai une difficulté de corriger cette erreur, est-ce que vous peuvent m'aider de trouver le cause de cette erreur? Sinon, est-ce que vous aves un programme qui lit les niveaux de gris des pixels de l'image commme entree pour l'algorithme EM. Merci beaucoup en avance et je vous remercie pour tous l'aide que vous avez me donne.

Code:

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
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
 
import java.io.*;
import java.awt.*;
import java.util.*;
import java.awt.image.PixelGrabber;
 
 
 
public class ImageTest {
public static byte[] readImage(String file){
Image image = Toolkit.getDefaultToolkit().getImage(file);
MediaTracker tracker = new MediaTracker (new Component() {});
tracker.addImage(image,0);
try {
tracker.waitForID(0);
PixelGrabber grabber = new PixelGrabber(image, 0, 0, -1, -1, false);
if(grabber.grabPixels()){
int nc = grabber.getWidth();
int nr = grabber.getHeight();
if (bytesAvailable(grabber)){
byte[] data = (byte[])grabber.getPixels();
double[] d = getData(data);
}
 
}
 
}
 
//MediaTracker tracker = new MediaTracker (new Component() {});
//tracker.addImage(image,0);
//try{tracker.waitForID(0);}
catch (InterruptedException e) {e.printStackTrace();}
return d[];
}
 
public static final boolean bytesAvailable(PixelGrabber pg){
return pg.getPixels() instanceof byte[];
}
//convert image byte data to double
public static double[] getData (byte[] data) {
  ByteArrayInputStream byte_in = new ByteArrayInputStream (data);
  DataInputStream data_in = new DataInputStream (byte_in);
  double[] d = data_in.readDouble ();
return d;
}
 
public static void main(String[] argv) {
if (argv.length>0) { double[] I = readImage(argv[0]);
 
        // initial guess for the laws parameters (uniform)
        int G=2;
 
        Law[] laws = new Law[G];
 
        laws[0] = new NormalDistribution(30.0,1.0);
 
        laws[1] = new NormalDistribution(130.0,2.0);            
 
        // perform EM algorithm
 
        double pi[] = Algorithm.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]);
    }
 
}
 
}

1ere remarque sur le principe : convertir la couleur d'un pixel en "double" ca n'est pas toujours évident. Un pixel etant constitué de 3 composantes (R,G,B), on calcule généralement sa valeur d'intensité lumineuse : I=0.299*R + 0.587*G + 0.114*B

2eme remarque sur le code : pour lire une image en Java, il y a plu simple que le MediaTracker : la classe ImageIO.

3eme remarque sur l'algo : le paramètre initial "sigma" pour les 2 lois est trop faible : une variance de 1 ou 2, sur une plage de [0-255] ca ne fait pas beaucoup. La plupart des pixels auront au départ une proba quasi nulle d'appartenir a l'une des lois.

D'ou ma version de ton code :

Code:

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
55
56
57
58
import java.awt.image.BufferedImage;
import java.io.File;
import java.io.IOException;
 
import javax.imageio.ImageIO;
 
public class ImageTest {
 
    public static double[] readImage(String pathname) throws IOException {
        // read image
        File f = new File(pathname);
        BufferedImage bi = ImageIO.read( f );
 
        // convert pixels to double[]
        int W = bi.getWidth(), H=bi.getHeight();
        double[] d = new double[W*H];
        int i=0;
        for(int y=0;y<H;y++)
            for(int x=0;x<W;x++) {
                // get R,G,B components and convert to intensity
                int rgb = bi.getRGB(x, y);
                int r = (rgb>>16)&0xFF, g=(rgb>>8)&0xFF, b=(rgb>>0)&0xFF;
                d[i++] = 0.299*r + 0.587*g + 0.114*b;
            }
        return d;
    }
 
    public static void main(String[] argv) {
        if (argv.length > 0) {
 
            double[] x;
            try {
                x = readImage(argv[0]);
            } catch (IOException e) {
                e.printStackTrace();
                return;
            }
 
            // initial guess for the laws parameters (uniform)
            int G = 2;
 
            Law[] laws = new Law[G];
 
            laws[0] = new NormalDistribution(30.0, 10.0);
            laws[1] = new NormalDistribution(130.0, 10.0);
 
            // perform EM algorithm
 
            double pi[] = Algorithm.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]);
        }
    }
 
}

27/05/2010, 11h50
hsalam

Merci beaucoup, vous êtes très gentil. Mais mon image est déja en niveau de gris! je ne besoin pas de convertir les components en niveaux de gris. alors j'enlever le partie de conversion? ou quoi? et pour executer le programme je tape sur la terminale : javac ImageTest 46.jpeg?(46.jpeg est le nom de mon image placé dans la même repertoire)
Merci
27/05/2010, 12h13
pseudocode
Citation:

Envoyé par hsalam

Mais mon image est déja en niveau de gris! je ne besoin pas de convertir les components en niveaux de gris. alors j'enlever le partie de conversion? ou quoi?

Non, ce code fonctionne aussi avec les images en niveaux de gris.

Citation:

et pour executer le programme je tape sur la terminale : javac ImageTest 46.jpeg?(46.jpeg est le nom de mon image placé dans la même repertoire)

Il faut d'abord compiler le code (avec javac), puis l'executer (avec java):
```
javac  ImageTest.java
java  ImageTest  46.jpeg
```
01/06/2011, 04h33
Qt forever

j'ignorais qu'il ya une implémentation de EM sur developpez ^^" , un inconvénient de compter sur google ts les jours

je n code pas par java mais j'ai analysé votre programme , j'ai une simple question est-ce que votre programme traite les dimension , car j'ai vu que vous avez donné au X(i) qu'une seul valeur , et s'il admet n'importe quel nombre de clusters ( ya pas de limites ! )

( dans mon cas chaque X(i) est représenté par D dimensions
X(i) = ( v1, v2 ,...,v2) , actuellement je travail avec D = 2 )

merci pour ce code qui ma bcp aider a la conception :ccool:
01/06/2011, 10h11
pseudocode

Dans le code que j'ai posté, un échantillon est un réel (double). Pour utiliser d'autre type d'échantillons, il faut changer les signatures des méthodes (et le code associé naturellement).

Le principe de l'algo EM ne change pas. Les seules modifications du code de l'algo concernent le changement de signature.

Le gros changement est dans l'implémentation des Lois, qui sont bien sur différentes pour une loi normale 1D et une loi normale 2D
03/06/2011, 19h25
Qt forever

merci bcp pour votre explication

maintenant après avoir " les means et les poids " de clusters je pense d'implémenter votre image que vous avez met dans "" ce poste ""

pouvez vous svp m'orienter pour pouvoir l'implémenter ? les outils nécessaires ( si vous avez utilisé un outil(logiciel , biblio ... ) ou les paramètres d'entrées que vous donnez à l'algo pour la dessiner !

1 pièce(s) jointe(s)

Bonjour,

Je me permet de "déterrer" ce sujet car j'ai actuellement besoin d'un tel algorithme pour une approximation. Je me suis donc permis de réimplémenter votre travail en Python, le code étant tout à fait similaire mais dans les standards de Python. Je rencontre toutefois quelques problèmes.

J'ai l'histogramme suivant, que je souhaiterais approcher avec une somme de deux lois normales :
Pièce jointe 259105

On peut très clairement voir 2 modes. Le code avant lancement de l'algorithme est le suivant :

Code:

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
probabilities = [0.6236059479553904, 0.637646554189305, 0.3711181012296254, 0.8209465255933658, 0.05173005433228524, 0.3573062625107237, 0.13344295110094412, 0.14206462682299148, 0.04645410351730055,
                 0.6488132685158708, 0.17620817843866216, 0.2253645982270521, 0.04941378324277953, 0.2106519874177867, 0.17899628252788152, 0.17989705461824468, 0.22883900486131004, 0.2938232770946526,
                 0.19019159279382364, 0.28713182728052633, 0.3198598798970549, 0.505161567057478, 0.41765799256505587, 0.5947240491850156, 0.5154132113239918, 0.9011152416356876, 1.1022447812410636,
                 1.1028452959679726, 1.2561910208750358, 1.786231055190163, 1.9953817557906774, 1.9962396339719757, 1.7160566199599654, 1.9174006291106662, 1.6134114955676293, 1.7424363740348872,
                 1.1333857592221905, 0.5862453531598513, 0.58684586788676, 0.153846153846153, 0.4595224478124107, 0.4082499285101515, 0.4287675150128684, 0.4460680583357164, 0.4667143265656279,
                 0.4875464684014872, 0.4679582499285101, 0.4169002001715756, 0.4576923076923076, 0.4662710895052904, 0.4512582213325709, 0.4437517872462113, 0.4724192164712611, 0.4731913068344295,
                 0.889791249642551, 0.8183299971404061, 0.8190591935945095, 1.0605519016299687, 1.0772376322562196, 1.1741778667429226, 1.207034601086646, 1.1755075779239348, 1.6416356877323421,
                 1.3533886188161282, 1.7631112382041754, 1.9885616242493562, 2.0775664855590503, 2.102316271089505, 2.1190877895338858, 2.424578209894195, 2.5295824992851013, 2.498269945667715,
                 2.531040892193308, 2.2991135258793247, 2.227695167286245, 2.2204460966542747, 1.9724764083500141, 1.8287818129825564, 1.6210037174721188, 1.3330140120102945, 1.3976980268801826,
                 1.4227051758650273, 1.4634114955676292, 1.223591649985702, 1.0959822705175866, 0.9683728910494712, 0.9289533886188162, 1.0179582499285103, 0.6656705747783817, 0.5460823563054048,
                 0.3863454389476695, 0.9726622819559625, 0.7164998570203032, 0.8112810980840718, 0.6298541607091793, 0.6760651987417787]
 
# Value repartition based on above probabilities.
ev_list=[]
ev_number = 3300
for i in range(len(probabilities)):
    number = int(ev_number * (probabilities[i] / 100))
    ev_list.append(number)
print ev_list
 
# Construct histogram.
histo = []
for i in range(len(ev_list)):
    value = ev_list[i]
    for j in range(value):
        histo.append(i)
 
# Compute the algorithm.
law_1 = NormalLaw(50.0, 15.0)
law_2 = NormalLaw(10.0, 25.0)
laws = [law_1, law_2]
results = em_algorithm(histo, laws, 5000)
print"Law 1 var: {}, mean: {}\nLaw 2 var: {}, mean: {}".format(law_1.variance, law_1.mean, law_2.variance, law_2.mean)
 
# Plot the histogram
x = np.linspace(0, 96, 200)
plt.hist(histo, x)
plt.show()

Je passe mon historique de valeurs en paramètre, est-ce correct ? Car après l'exécution de l'algorithme, j'ai les valeurs suivantes :

Code:

1
2
Law 1 var: 61.7435834864, mean: 57.1325894322
Law 2 var: 61.7435834864, mean: 57.1325894322

Ce qui me semble très surprenant.

Merci de votre aide.

30/03/2017, 11h10
Ecterion

J'ai réussi à régler le problème et j'ai également utiliser votre algorithme en Java. Toutefois, comment faites-vous pour tracer les courbes en Java ?