//// 1d example EM algorithm
 clear; 
 x0=-10:0.1:10;
 var1=1; var2=1; mu1=-2; mu2=2;
 
 function y=normal(x,mu,var)
   y=exp(-(x-mu).^2/(2*var))/sqrt(2*%pi*var);
 endfunction
 
 for trial=1:10
////plot distribution
    clf; 
    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');
//// data
    x=[2*rand(50,1,'normal')-1;0.5*rand(50,1,'normal')+4;];
//// recogintion
    c=1*(normal(x,mu1,var1)>normal(x,mu2,var2));
//// maximization
    mu1=sum(x(c>0.5))/sum(c);
    var1=sum((x(c>0.5)-mu1).^2)/sum(c);
    mu2=sum(x(c<0.5))/(100-sum(c));
    var2=sum((x(c<0.5)-mu2).^2)/(100-sum(c));
 end