function [mu, Sigma, mixcoef] = learn_mog(data, K, cov_type, max_iter, mu, Sigma, mixcoef)
% LEARN_MOG  Learn a mixture of K Gaussians using EM.
% [mu, Sigma] = LEARN_MOG(data, K, cov_type)
%
% Inputs:
% data(:,l) (col l) is the l'th observation vector
% K is the number of centers
% cov_type = 'full', 'diag', or 'spherical'
%
% Outputs:
% mu(:,k), Sigma(:,:,k) is the k'th mean and covariance
% mixcoef(k) is the prior for the k'th component

previous_loglik = -inf;
loglik = 0;
converged = 0;
num_iter = 1;
threshold = 1e-2;

[O numex] = size(data);

if nargin < 5
  mu = repmat(mean(data')', [1 K]);
  Sigma = repmat(cov(data'), [1 1 K]);
  mixcoef = normalise(rand(K,1)); % have to break symmetry somewhere!
end

while (num_iter <= max_iter) & ~converged

  % E step
  
   % lik(l,k) = Pr(y(l) | Q=k) where Q is the hidden mixture variable
  lik = zeros(numex, K);
  for k=1:K
    lik(:,k) = gaussian_prob(data, mu(:,k), Sigma(:,:,k));
  end
  % post(l,k) = Pr(Q=k|y(l)) = Pr(y(l) | Q=k) Pr(Q=k) / Z = responsibility of k'th component
  post = mk_stochastic(lik .* repmat(mixcoef(:)', [numex 1]));
  % loglik = sum_l log (sum_k lik(l,k) * mixcoef(k))
  loglik = sum(log(lik * mixcoef(:)));
  
  % M step
  w = sum(post, 1); % w(k) = sum_l post(l,k)
  mixcoef = normalise(w(:));
  mu = (post' * data')' ./ (repmat(w(:)', [O 1]));

  if cov_type(1) == 'f' | cov_type(1) == 'd'
    for k=1:K
      d = data - repmat(mu(:,k), [1 numex]); % d(o,l) = y(o,l) - mu(o)
      d = d .* repmat(sqrt(post(:,k)'), [O 1]); 
      Sigma(:,:,k) = d*d' / w(k); % S = sum_l post(l) d(:,l) d(:,l)'
      if cov_type(1) == 'd'
	Sigma(:,:,k) = diag(diag(Sigma(:,:,k)));
      end
    end
  elseif cov_type(1) == 's'
    d = data - repmat(mu(:,k), [1 numex]); % d(o,l) = y(o,l) - mu(o)
    d = diag(d'*d); % d(o) = sum_l d(o,l) d(o,l)
    for k = 1:K
      v(k) = (post(:,k)'*d);
      Sigma(:,:,k) = (1/O) * (v(k) / w(k));
    end
  end
  
  fprintf(1, 'iteration %d, loglik = %f\n', num_iter, loglik);
  num_iter =  num_iter + 1;

  converged = em_converged(loglik, previous_loglik);
  previous_loglik = loglik;
end


