% mog.m - trains a mixture of Gaussians model
%
% assumes X defined

[N L]=size(X);

% number of Gaussians
K=3;

% parameters
mu=2*randn(N,K);
gamma=randn(K,1);

clf
subplot(122)
plot(X(1,:),X(2,:),'.')
hold on
nphi=20;
phi=2*pi*[0:nphi-1]'/(nphi-1);
h=plot(cos(phi)*ones(1,K)+ones(nphi,1)*mu(1,:),...
    sin(phi)*ones(1,K)+ones(nphi,1)*mu(2,:),'EraseMode','xor');

num_trials=10000;
eta=.1/L;

for t=1:num_trials

    % compute priors
    Pa=exp(gamma);
    Pa=Pa/sum(Pa);

    % compute P(alpha|x)
    for alpha=1:K
        dX=X-mu(:,alpha)*ones(1,L);
        Pax(alpha,:)=exp(-0.5*sum(dX.^2))*Pa(alpha);
    end
    Pax=Pax*diag(1./sum(Pax));
    
    % update parameters
    for alpha=1:K
        dX=X-mu(:,alpha)*ones(1,L);
        mu(:,alpha)=mu(:,alpha)+eta*dX*Pax(alpha,:)';
        gamma(alpha)=gamma(alpha)+eta*(sum(Pax(alpha,:))-L*Pa(alpha));
        set(h(alpha),'XData',cos(phi)+mu(1,alpha),...
            'YData',sin(phi)+mu(2,alpha))
    end
    subplot(121)
    bar(Pa)
    drawnow
    
end
hold off
