mix.m

EM算法处理高斯混和模型

function [Z] = mix(x,k,iterations);
% Simple function to do EM for a mixture of Gaussians.
% by Carl Edward Rasmussen, 2006-01-18.

D = 2;

% Initialise parameters

n = length(x);        % number of observations
p = ones(1,k)/k;      % mixing proportions
mu = randn(D,k);      % means
for i=1:k
  s2(1:D,1:D,i) = -10*diag(log(rand(D,1))); % variances
end
Z = [];

for t=1:iterations

  % Do the E-step:

  figure(1)
  hold off
  plot(x(:,1),x(:,2),'+');
  hold on
  for i=1:k
    draw_ellipse(mu(:,i),s2(:,:,i),'r');
  end

  for i=1:k
    di = x-repmat(mu(:,i)',n,1);
    Q(:,i)=p(i)*exp(-sum(di'.*(inv(s2(:,:,i))*di'),1)/2)'/sqrt(det(s2(:,:,i)));
   end
   Z(t) = sum(log(sum(Q/(2*pi),2)));

   figure(2)
   plot(Z)

% waitforbuttonpress

   fprintf('Iteration: %i  log likelihood: %4.3e\r', t, Z(t));
   Q = Q ./ (sum(Q,2)*ones(1,k));    % Normalise

   % Do the M-step:

   for i=1:k
     di = x-repmat(mu(:,i)',n,1);
     mu(:,i) = sum(x.*repmat(Q(:,i),1,D),1)'/sum(Q(:,i))';
     di = x-repmat(mu(:,i)',n,1);
     s2(:,:,i)=(di'*(di.*repmat(Q(:,i),1,D)))/sum(Q(:,i));
    end
    p = sum(Q,1)/sum(sum(Q));

end