
/////////////////////////////////////////////////////////////
//ALGORITHME EXPECTATION MAXIMISATION 
//Auteur: F. Richard
//Date: 24/04/2006
/////////////////////////////////////////////////////////////


//Ces fonctions font l'estimation des parametres d'un modele de melange avec l'algorithme EM.

//Modele statistique:
// Y=(Yk)k variables aleatoires iid selon la loi de densite 
// f(y) =  (1-p) f0(y) + p f1(y)
// avec:
// - (1-p) et p sont les probabilites qu'un individu appartienne aux classes 0 et 1, respectivement.
// - f0 est la densite de loi sur la classe 0 (gaussienne de moyenne m0 et de variance sig0)
// - f1 est la densite de loi sur la classe 1 (gaussienne de moyenne m1 et de variance sig1)




//Les parametres estimes sont p, m0, m1, sig0, sig1
//ainsi que la probabilite d'appartenance a la classe 1 pour chaque individu (proba1).



//PARAMETRES D'ENTREE
///////////////////////

//Y: echantillon dont la loi de densite est f.
//nb_iter : nb maximum d'iterations

//EN SORTIE
////////////
//m0, sig0: estimation des parametres de la densite de loi de la classe 0
//m1, sig1: estimation des parametres de  la densite de loi de la classe 1
//p: estimation de la probabilite de presence d'un individu de la classe 1
//proba1:  estimation des probabilites d'appartenance de chq individu a la classe 1


function [m0,sig0,m1,sig1,p,proba1]=algo_em(Y,nb_iter),


//initialisation des parametres a estimer
N=size(Y,"c")*size(Y,"r");
m0=min(Y);
m1=max(Y);
p=0.5;
sig0=st_deviation(Y)^2;
sig1=sig0;


for iter=1:nb_iter,
	m1_ref=m1;
	s0=sqrt(sig0);
	s1=sqrt(sig1)	

	//nouvelle estimation des probabilites d'appartenance a la classe 0 pour chq individu
	proba1=p*1/sqrt(2*%pi)/s1*exp(-0.5*((Y-m1)/s1).^2);
	proba1=proba1./((1-p)*1/sqrt(2*%pi)/s0*exp(-0.5*((Y-m0)/s0).^2)+proba1);


	//nouvelle estimation des parametres de loi
	if (p==0),
	   m0=mean(Y);
	   m1=0;
	else,
	   norm1=sum(proba1);
	   norm0=N-norm1;
	   //estimation des moyennes sur chq classe	
           m0=sum((1-proba1).*(Y))/norm0;	
	   m1=sum(proba1.*(Y))/norm1;
	   //estimation des variances sur chq classe
	   sig0=sum((1-proba1).*(Y-m0).^2)/norm0;
           sig1=sum(proba1.*(Y-m1).^2)/norm1;
	end;

	//nouvelle estimation de la probabilite d'apparition d'un individu de la classe 1
	p=mean(proba1);

	//critere d'arret (stabilisation de l'estimation de m1)
	if (abs(m1_ref-m1)<1) iter=nb_iter; end;	
end;

endfunction;



//EM: version 2, on n'estime que m1, sig1, et p
function [m1,sig1,p,proba1]=algo_em_v2(Y,m0,sig0,nb_iter),


//initialisation des parametres a estimer
N=size(Y,"c")*size(Y,"r");
m1=max(Y);
p=0.5;
sig1=sig0;
s0=sqrt(sig0);

for iter=1:nb_iter,
	m1_ref=m1;
	s1=sqrt(sig1)	

	//nouvelle estimation des probabilites d'appartenance a la classe 0 pour chq individu
	proba1=p*1/sqrt(2*%pi)/s1*exp(-0.5*((Y-m1)/s1).^2);
	proba1=proba1./((1-p)*1/sqrt(2*%pi)/s0*exp(-0.5*((Y-m0)/s0).^2)+proba1);


	//nouvelle estimation des parametres de loi
	if (p==0),
	   m1=0;
	else,
	   norm1=sum(proba1);
	   //estimation des moyennes sur chq classe	
	   m1=sum(proba1.*(Y))/norm1;
	   //estimation des variances sur chq classe
           sig1=sum(proba1.*(Y-m1).^2)/norm1;
	end;

	//nouvelle estimation de la probabilite d'apparition d'un individu de la classe 1
	p=mean(proba1);

	//critere d'arret (stabilisation de l'estimation de m1)
	if (abs(m1_ref-m1)<1) iter=nb_iter; end;	
end;

endfunction;

//EM: version 3, on n'estime que m1
function [m1,proba1]=algo_em_v3(Y,m0,sig0,sig1,p,nb_iter),


//initialisation des parametres a estimer
N=size(Y,"c")*size(Y,"r");
m1=max(Y);
s0=sqrt(sig0);
s1=sqrt(sig1)	



for iter=1:nb_iter,
	m1_ref=m1;

	//nouvelle estimation des probabilites d'appartenance a la classe 0 pour chq individu
	proba1=p*1/sqrt(2*%pi)/s1*exp(-0.5*((Y-m1)/s1).^2);
	proba1=proba1./((1-p)*1/sqrt(2*%pi)/s0*exp(-0.5*((Y-m0)/s0).^2)+proba1);


	//nouvelle estimation des parametres de loi
	if (p==0),
	   m0=mean(Y);
	   m1=0;
	else,
	   norm1=sum(proba1);
	   //estimation des moyennes sur chq classe	
           m1=sum(proba1.*(Y))/norm1;
	end;

	//nouvelle estimation de la probabilite d'apparition d'un individu de la classe 0
	p=mean(proba1);

	//critere d'arret (stabilisation de l'estimation de m1)
	if (abs(m1_ref-m1)<1) iter=nb_iter; end;	
end;

endfunction;



