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)
%
% See also: mogP, mog_dd, gauss_dd

% Copyright: D.M.J. 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 = 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.');
  end
end

return


function 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