function [means, covs, priors] = mogEM(x,k,ctype,reg,nriters)% MOGEM Optimize Mixture of Gaussian by EM%%       [means, covs, priors] = mogEM(x,k,ctype)%% Optimize a mixture of Gaussians with k centers on data x. The type% of covariance matrix can be given in variable ctype:%      ctype = 'sphr'  : diagonal cov. matrix with equal values%      ctype = 'diag'  : diagonal cov. matrix%      ctype = 'full'  : full cov. matrix%% The means, the covariance matrices and the prior probabilities for% each of the clusters is returned.% %     [means, covs, priors] = mogEM(x,k,ctype,reg,nriters)%% The user can add a regularization weight reg, to avoid singular% covariance matrix estimations. Finally the number of iterations can% be given (default: nriters=100)% Copyright: D. Tax, R.P.W. Duin, davidt@ph.tn.tudelft.nl% Faculty of Applied Physics, Delft University of Technology% P.O. Box 5046, 2600 GA Delft, The Netherlandsif nargin<5  nriters = 100;endif nargin<4  reg = 0.01;end% initialize the stuff:[n,d] = size(x);COVWIDTH = 1;%random start meansI = randperm(n);means = x(I(1:k),:);%priorsD = distm(x,means);[md,I] = min(D,[],2);nr_x_i = sum(full(ind2vec(I))',1);priors = nr_x_i./sum(nr_x_i);%covariance matrix:.switch ctype  case 'sphr'    if (k==1)      covs = COVWIDTH;    else      sD = sort(distm(means));      covs = sD(2,:)';      covs(covs<eps) = COVWIDTH;    end  case 'diag'    covs = zeros(k,d);    for i=1:k      x_i = x(find(I==i),:);      dif = x_i - ones(size(x_i,1),1)*means(i,:);      covs(i,:) = sum(dif.*dif)/nr_x_i(i);    end    covs(covs<eps) = COVWIDTH;  case 'full'    covs = zeros(k,d,d);    for i=1:k      x_i = x(find(I==i),:);      dif = x_i - ones(size(x_i,1),1)*means(i,:);      covs(i,:,:) = (dif'*dif)/nr_x_i(i);      if rank(squeeze(covs(i,:,:)))<d        covs(i,:,:) = squeeze(covs(i,:,:)) + COVWIDTH*eye(d);      end    endend% now we start to optimize:for iter=1:nriters  % calculate the old P:  P = mogP(x,means,covs,priors);  % normalize it:  sumP = sum(P,2); sumP(sumP==0) = 1;  normP = P./(sumP*ones(1,k));  % update the params:  new_priors = sum(normP);  new_means = normP' * x;  % thus:  priors = new_priors/n;  means = new_means ./ (new_priors' * ones(1,d));  switch ctype    case 'sphr'      D = distm(x,means);      news = zeros(k,1);      for i=1:k        news(i) = (normP(:,i)'*D(:,i));      end      covs = (news./new_priors')/d + reg*ones(k,1);    case 'diag'      for i=1:k        dif = x - ones(n,1)*means(i,:);        covs(i,:) = sum((dif.*dif).*(normP(:,i)*ones(1,d)))./new_priors(i) +...                    reg*ones(1,d);      end    case 'full'      for i=1:k        dif = (x - ones(n,1)*means(i,:)) .* sqrt(normP(:,i)*ones(1,d));        covs(i,:,:) = (dif'*dif)/new_priors(i) + reg*eye(d);      end     otherwise      error('oef2');  endendreturn