%%EM example

% from a mixture model
% p*N(mu_1,v=1)+(1-p)N(mu_2,v=1)
% estimate p, mu_1, mu_2
clear all;

muTrue_1 = 4;  %% true means
muTrue_2 = 2;

pTrue=.7; %% true mixing proportion

v=1; %% known variance

N=1000; % number of data points


%simulate data (2 component mixture)
for(i=1:N)
    if(rand<pTrue)
        y = randn*sqrt(v) + muTrue_1;
    else
        y = randn*sqrt(v) + muTrue_2
    end
    Y(i) = y;
end


%%plot data
figure(1)
clf
hist(Y,30);
title('Simulated Mixture Data')
pause(1)

%EM for parameter estimation

%initial values
mu(1) = Y(ceil(rand*N));
mu(2) = Y(ceil(rand*N));;
p(1)=rand
p(2) = 1-p(1);  % this could be ommitted

iters=100;
for(i=1:iters)
    %update mixing proportion
    for(j=1:N)
        for(k=1:2)
            pp(j,k) = exp(-normlike([mu(k),v],Y(j)))*p(k);   
        end
        pp(j,:) = pp(j,:)./sum(pp(j,:));  % simple averaging
    end
    
    for(k=1:2)
        p(k) = mean(pp(:,k));
    end
    
    %update the means
    for(k=1:2)
        mu(k)=0;
        for(j=1:N)
            mu(k) = mu(k)+pp(j,k)*Y(j); % weighted averaging
        end
        mu(k) = mu(k)/sum(pp(:,k));
    end
    fprintf('------------------------\n')
    mu
    p 
end