% Chad Myers
% EE 571
% File: em.m
% ---------------------------------
% The em function implements the pure RBF EM algorithm.  It takes as input
% the number of classes to assume when running the algorithm, and a vector of 
% the input.  This code is largely composed of Prof. Kung's em function.

function [h, mu, sigmas, P] = em(M, Y);
markers = ['o' 'x' '+' '*' 'd' 'v' '^' '<' '>' 'p' 'h'];
maxEpoch = 20; [T,D] = size(Y);
h = ones(T,M)/M;  P = ones(M,1)/M; 
for j=1:M,
	k = round(rand(1)*T);
    mu(j,1:D) = Y(k,1:D);
	sigmas(j,1:D) = abs(mean(mu(j,:)))/2*ones(1,D);   
end

for epoch=1:maxEpoch,
   % E-step
   for t=1:T,
      for k=1:M,
         g(k) = P(k)*Gaussian(Y(t,:),mu(k,:),sigmas(k,:));
      end
      s0 = sum(g);
      for j=1:M,
         if(s0 ~= 0) h(t,j) = g(j)/s0;         
         else h(t,j)=0;
         end
      end
   end
   
   % M-step
   for j=1:M,
      sigmasum=0;
      for i=1:D,
         s1=0; s2=0; s3=0;
         for t=1:T,
         s1 = s1 + h(t,j)*Y(t,i);
         s2 = s2 + h(t,j);      
         end
         
         if(s2 ~= 0) mu(j,i) = s1/s2;
         else mu(j,i)=0;
         end 
         P(j) = s2/T;
         s3=0;
         
      end    
      
      for t=1:T,
         s3 = s3 + h(t,j)*(norm((Y(t,:)-mu(j,:))))^2;         
      end
         
      for k=1:D,
          if(s2 ~= 0) sigmas(j,k)=s3/s2;
          else sigmas(j,k) = 0;
          end
      end
          
  end
   %mu

   
   %continued
   % Animation. Plot cluster by cluster 
% cla reset;
%    for j=1:M,
%       n=1; C=[]; Z=[];
%       for t=1:T,
%          [maxh,maxpos] = max(h(t,:));
%          if maxpos == j
%             C(n)=j+4;
%             Z(n,:)=Y(t,:);
%        	 	n = n+1;    
%       	end 
%       end
%       mindex = mod(j,length(markers));
%       if mindex == 0
%          mindex = length(markers);
%       end
%   		if(size(Z) ~= 0) scatter(Z(:,1),ones(size(Z,1),1),10,C,markers(mindex));
%       %if(size(Z) ~= 0) scatter(Z(:,1),Z(:,2),10,C,markers(mindex));
%         end
%         
%    	hold on;
%       Cmu=ones(length(mu),1); 
%       scatter(mu(:,1),ones(size(mu,1),1),10,Cmu,'s');
%       
%       %scatter(mu(:,1),mu(:,2),10,Cmu,'s');
% 	end
% 	pause;
%     %pause(.01);   
end
