Discussion :

[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 :

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
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 :

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
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 :

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
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 :

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
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 :

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
// 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 )

[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]

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]

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