%% 1d example EM algorithm
clear;	hold on; x0=-10:0.1:10;
var1=1; var2=1; mu1=-2; mu2=2; phi=0.5
normal= @(x,mu,var) exp(-(x-mu).^2/(2*var))/sqrt(2*pi*var);

%% data
x=[2*randn(50,1)-1;0.5*randn(100 ,1)+4];
while 1
    %%plot distribution
    clf; hold on; ylim([0 1]);
    plot(x0, normal(x0,-1,4),'k:');
    plot(x0, normal(x0,4,.25),'k:');
    plot(x0, normal(x0,mu1,var1),'r');
    plot(x0, normal(x0,mu2,var2),'b');
    waitforbuttonpress ;
    %% expectation
    p1=normal(x,mu1,var1)*phi;
    p2=normal(x,mu2,var2)*(1-phi);
    nrm=p1+p2; 
    p1=p1./nrm; p2=p2./nrm;
    %% maximization
    phi=1/size(x,1)*sum(p1);
    mu1=sum(x.*p1)/sum(p1); var1=sum(p1.*(x-mu1).^2)./sum(p1);
    mu2=sum(x.*p2)/sum(p2); var2=sum(p2.*(x-mu2).^2) ./sum(p2);
end