%----------------
randn('seed',1) %set the seeds to have constancy of results
rand ('seed',1)
%--------------generate the data --------------------------
y=[];
n=150;
for i=1:n
    add = 2 * randn - 5; %N(-5,4) 50% comp 1
    if rand < 0.2
       add =  randn + 2; %N(2,1) 20% comp 2
    elseif rand > 0.7
       add = 0.5*randn; %N(0,0.25) 30% comp 3
    end
   y = [y add];
end
% we got 150 obsercations from
% 0.5  N(-5, 2^2) + 0.3 N(0, 0.5^2) + 0.2 N(2,1).
%------- forget the weights now ------------------------
g = 3; %number of distributions
pi_current = 1/3 * ones(1, 3); %------ignorance proposal---
%---------
mus    = [ -5   2   0 ];   %---- we know the distributions
sig2s  = [4   1   0.25 ];  %---- exactly, so means, vars are known
%---------------EM starts!-------------------------------------
for  k = 1 : 10            %---- # of E-M cycles. 10 is enough!
szs=[]; z=zeros(g,n);
%E-Step, conditional expectations of z's
for i=1:g
    for j=1:n
        z(i,j) = pi_current(i) * 1./(sqrt(2*pi*sig2s(i))) * ...
            exp( -(y(j)-mus(i))^2/(2*sig2s(i)) );
    end
end
%  z(i,j)  are probabilities that observation j is coming
%  from the component i.    But the probs are not normalized.
%  The following cycle will compute the normalizing constant.
for j=1:n
    sz = 0;
    for i=1:g
        sz = sz + z(i,j);
    end
    szs =[szs, sz];
end
norm = [szs' szs' szs']'; % norm needed to normalize z_{ij}'s
zk = z ./ norm;           % now the z's are normalized
%-------------------------------------
% M-Step, just plug in zk's
pi_new = sum(zk')/n; % and repeat...
pi_current = pi_new;
end
pi_current

%%> pi_current =  0.5599    0.1765    0.2636