function [a,m,cm, lp]=EM(x, N, iter, m, cm)
% Expectation Maximazation algorithn
% x: the samples
% N: how many Gaussian

[dimx,lenx] = size(x);

a = ones(1, N)/N;                  % initialization
%m = .03*randn(dimx, N);
if(nargin<5)
   cm = repmat(eye(dimx), 1, N);
end
if(nargin<4)
  m = x(:,1:N);
end

disp('ok');
an = zeros(1, N);                  % new values
mn = zeros(dimx, N);
cmn = repmat(eye(dimx), 1, N);

flagstop=0;
lp = zeros(1,iter);
disp('ok');
liter=1;
while(flagstop==0)
    
    nf = ppe(x, a, m, cm);
    lp(liter) = sum(nf);
    if(liter>1) 
        thres=(lp(liter)-lp(liter-1))/lp(liter);
        if(thres<2*10^(-5))
            flagstop=1;
        end
    end
    fprintf(1,'Iter: %d Likelihood : %f\n', liter , lp(liter));
    
    for i = 1:N
        vv = gau(x, m(:,i), cm(:, (i-1)*dimx+1 : i*dimx))./nf;
        an(i) = a(i)*sum(vv, 2) / lenx;
      
        mn(:,i) = a(i)*sum(x.*repmat(vv, dimx, 1), 2) / (lenx*an(i));
  
        for k = 1:dimx
            for kk = 1:dimx
                cmn(k, (i-1)*dimx+kk) = a(i)*sum(vv .* (x(k,:)-mn(k,i)).*(x(kk,:)-mn(kk,i))) / (lenx*an(i));
            end
        end

    end
    a = an;
    m = mn;
    cm = cmn;
    liter=liter+1;
end
lp = lp(1:liter-1);