mogem.m

模式识别工具箱,希望对大家有用!

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, [e-mail]
% Faculty of Applied Physics, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands


if nargin<5
  nriters = 100;
end
if nargin<4
  reg = 0.01;
end

% initialize the stuff:
[n,d] = size(x);
COVWIDTH = 1;

%random start means
I = randperm(n);
means = x(I(1:k),:);

%priors
D = 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
    end
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);
  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');
  end
end

return