function EM;

%--------------------------------------------------------------------------
% Function: EM;
% Scope   : Illustra il funzionamento di EM per un problema di stima di due
%           rette passanti per l'origine con rumore gaussiano sui dati.
% Input   : None.
% Output  : None.
% Notes   : None.
% Author  : Iuri Frosio, frosio@dsi.unimi.it.
% Date    : 2007.12.04.
%--------------------------------------------------------------------------

% Genera i dati
P1=0.15;  m1=-2;
P2=1-P1; m2=1;
Sigma=25;
Nyis=300;
xis=zeros(1,Nyis);
yis=zeros(1,Nyis);
for nyi=1:Nyis
    xis(nyi)=randn(1)*100;
    if (rand(1)<P1)
        yis(nyi)=m1*xis(nyi)+Sigma*randn(1);
    else
        yis(nyi)=m2*xis(nyi)+Sigma*randn(1);
    end
end

% Inizializzazione casuale
Estm1=m1+randn(1)*10*m1;
Estm2=m2+randn(1)*10*m2;

% Plot time
Minxis=min(xis);
Maxxis=max(xis);
Minyis=min(yis);
Maxyis=max(yis);
figure(1); subplot(1,2,1);
           plot(xis,yis,'k.',...
                [Minxis Maxxis],Estm1*[Minxis Maxxis],'r',...
                [Minxis Maxxis],Estm2*[Minxis Maxxis],'b');
           xlabel('x'); ylabel('y'); grid on; title('Dati'); axis equal; axis([Minxis Maxyis Minyis Maxyis]);
           
% Numero di iterazioni
NIters=50;

% Iterazioni algoritmo EM
for nIter=1:NIters
    
    % Vecchia stima dei parametri
    m1_old=Estm1;
    m2_old=Estm2;
    
    % Passo expectation: calcolo p(Z|X,m_old)
    for nyi=1:Nyis
        p1_xi_old(nyi)=(Gaussian(yis(nyi),m1_old*xis(nyi),Sigma)*P1)/...
                       max(eps,(Gaussian(yis(nyi),m1_old*xis(nyi),Sigma)*P1+Gaussian(yis(nyi),m2_old*xis(nyi),Sigma)*P2));
        p2_xi_old(nyi)=(Gaussian(yis(nyi),m2_old*xis(nyi),Sigma)*P2)/...
                       max(eps,(Gaussian(yis(nyi),m1_old*xis(nyi),Sigma)*P1+Gaussian(yis(nyi),m2_old*xis(nyi),Sigma)*P2));
    end
    
    % Passo maximization
    Estm1=sum(p1_xi_old.*yis.*xis)/sum(p1_xi_old.*(xis.^2));
    Estm2=sum(p2_xi_old.*yis.*xis)/sum(p2_xi_old.*(xis.^2));
    
    % Plot time
    subplot (1,2,2);
    plot(xis,yis,'k.',...
         [Minxis Maxxis],Estm1*[Minxis Maxxis],'r',...
         [Minxis Maxxis],Estm2*[Minxis Maxxis],'b');
    xlabel('x'); ylabel('y'); grid on; title(['nIter = ' num2str(nIter)]); axis equal; axis([Minxis Maxyis Minyis Maxyis]);
    pause;   
    
end

%--------------------------------------------------------------------------
% Function: p=Gaussian(x,Mu,Sigma);
% Scope   : Calcola il valore in x della gaussiana con media Mu e deviazione
%           standard Sigma.
% Input   : x, Mu, Sigma.
% Output  : p(x|Mu,Sigma).
% Notes   : None.
% Author  : Iuri Frosio, frosio@dsi.unimi.it.
% Date    : 2007.12.04.
%--------------------------------------------------------------------------
function p=Gaussian(x,Mu,Sigma);

p=(1/(sqrt(2*pi)*Sigma))*exp(-0.5*(((x-Mu)/Sigma).^2));