// Evaluation d'une densité function f=gauss(m,sq,x)
f=exp((-(x-ones(x)*m)^2)./(2*sq))/sqrt(2*%pi*sq)
endfunction;

// Génération de l'échantillon
function X=echantillon(n,p,m,v)
X=[];
k=length(p);
z=grand(1,"mul",n,p(1:$-1)')';
for i=1:k
  X=[X grand(1,z(i),"nor",m(i),sqrt(v(i)))]
end;
endfunction;

// Fonction principale 
//alv : probas de mélange vraies 
//muv : moyennes vraies
//vv : variances vraies
//n : taille de l'échantillon
//iter : nombre d'itérations de l'algorithme
function em(alv,muv,vv,n,iter)
k=length(alv)
X=echantillon(n,alv,muv,vv)
al=ones(1,k)/k
mu=[1:k]
sq=ones(1,k)
h=zeros(n,k)
//Boucle principale
for i=1:iter
  //etape Expectation
    // A faire 

  //etape Minimisation
    // A faire

end
// Comparaison des paramètres
disp('al vrai')
disp(alv)
disp('al estimé')
disp(al)
disp('mu vrai')
disp(muv)
disp('mu estimé')
disp(mu)
disp('variance vraie')
disp(vv)
disp('variance estimée')
disp(sq)
// Tracé des deux densités
xbasc()
x=[min(muv)-5:.01:max(muv)+5]
Y=zeros(x)
Z=zeros(x)
for l=1:k
  Y=Y+al(l)*gauss(mu(l),sq(l),x)
  Z=Z+alv(l)*gauss(muv(l),vv(l),x)
end
plot2d(x,Y,5)
plot2d(x,Z,2)
legends(['Densité vraie','Estimation'],[2,5],1)
endfunction;