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)%% See also: mogP, mog_dd, gauss_dd% Copyright: D.M.J. 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 = sqeucldistm(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(sqeucldistm(means,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    end otherwise	 error('Covariance type is not defined.');end% 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);  I = find(new_priors==0); new_priors(I) = realmax/10;  new_means = normP' * x;  % thus:  priors = new_priors/n;  means = new_means ./ (new_priors' * ones(1,d));  switch ctype    case 'sphr'      D = sqeucldistm(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('Unknown covariance type requested.');  endendreturnfunction v=ind2vec(i)%IND2VEC Convert indices to vectors.%%	Syntax%%	  vec = ind2vec(ind)%%	Description%%	  IND2VEC and VEC2IND allow indices to be represented%	  either by themselves, or as vectors containing a 1 in the%	  row of the index they represent.%%	  IND2VEC(IND) takes one argument,%	    IND - Row vector of indices.%	  and returns sparse matrix of vectors, with one 1 in%	  each column, as indicated by IND.%%	Examples%%	  Here four indices are defined and converted to vector%	  representation.%%	    ind = [1 3 2 3]%	    vec = ind2vec(ind)%%	See also VEC2IND.% Mark Beale, 2-15-96.% Copyright 1992-2000 The MathWorks, Inc.% $Revision: 1.7 $vectors = length(i);v = sparse(i,1:vectors,ones(1,vectors));return