kmeans.m

最新的模式识别分类工具箱,希望对朋友们有用!

% ---------------------------------------------------------------- 
% FUNCTION   kmeans.m    kmeans algorithm.
% ---------------------------------------------------------------- 
% 		builds a block training set. 
% Usage: 
%   cluster_centers = kmeans(orig_cluster_centers, training set,
%                       precision);
% 
% orig_cluster: matrix, each ROW is a cluster center.
% training set: matrix, each ROW is a sample.
% Precision = stopping criterion: If the difference between two following
%		mse's is less than epsilon, or if the decay in 
%		the dispersion is less than the threshold, then
%		inquire if better results come from splitting or merging,
%		then, if no splitting or merging occurs, terminate.
% __________________________
%  Author: Vittorio Castelli
%  Copyright IBM T. J. Watson Research Center
%      February 6, 1995;	Last modified: February 6, 1995


function clCen =  kmeans(origCen, trSet,prec)


n_dist = 8;
maxcount = 40;
count = 0;
contatore = 0;
dis = 1;
firsttime = 1;
perc = 4;	% percentage of the training set shown in the plot
its1 = 1;
its2 = 2;	% indices of the coordinates of the pattern shown in
its3 = 3;	% plot3


[m,n] = size(trSet);
[n_cent,n1] = size(origCen);
if(its1 > n1)
   its1 = 1;
end

if(its2 > n1)
   its2 = n1/3;
end

if(its3 > n1)
   its3 = n1*2/3;
end

if n ~= n1,
   errormes(4);
   return
end

aaa = ones(1,n);
aa = ones(1,m);
clCen = origCen;

%centroid = mean(trSet);
%origdisp = mean(sum((trSet'-centroid'*aa).^2))
%splitdisp = origdisp*threshsplit;

% --------------------
% Initialize the plots

figure(1);
clg;
figure(2);
clg;
figure(3);
clg;
asdfg = max(size(trSet));
disp([its1, its2, its3]); 
plot3(trSet(1:perc:asdfg,its1),trSet(1:perc:asdfg,its2), ...
         trSet(1:perc:asdfg,its3),'.');
grid;
hold;
oldplCen = clCen;  
plot3(clCen(:,its1),clCen(:,its2),clCen(:,its3),'g*');

% --------------------
% Initialize the variables used in the loop

contin =1;
disper = ones(1,n_cent)*1024;
oldisper = disper;


while contin==1,

  oldCen=clCen;			% Initialize variables used
  contatore = contatore +1;	% for termination of loop;
  olddisper = disper;		%

  
  for ii = 1:n_cent,  % Calculate distances from centers
      distan(ii,:) = sum((trSet'-((clCen(ii,:))'*aa)).^2);
  end;
  

  [distan,ord] = sort(distan);	% Sort distances
  ord = ord';			% ord contains now the cluster to which 
  ord = ord(:,1);		% each element belongs
  chancent=0;			%
  
  for i=1:n_cent,			% Calculates the number of pts. 
        nelcloud(i) = sum(ord==i);      % in the cloud. If it is > 0
        if nelcloud(i) ==0,		% it moves the center 
             fprintf(' %g',i)		% to the centroid
          else
           clCen(i,:) =  sum((((ord==i)*aaa).*trSet))/nelcloud(i);
           disper(i) =  sum((ord == i).* ...
                   (sum(((clCen(i,:)'*aa)-trSet').^2))')/nelcloud(i);
        end
        d_i = sum((clCen(i,:)-oldCen(i,:)).^2); % Calculate how much the 
        chancent=max(chancent, d_i);		% centers moved in this pass
        disi(i) = d_i;         			% and record the max
   end
 
   figure(1);
   subplot(211)
   plot(1:length(disper), sqrt(disper),'o');
   title('Dispersion');
   subplot(212)
   plot(disi,'x')
   titlestr = sprintf('Changes, iter. %g',contatore);
   title(titlestr);
   
   count = count+1;
   dispersion(count) = sum(nelcloud.*disper)/m;


   if rem(count,n_dist)==0, % Plot the mse
      figure(2);
      plot(dispersion,'+');
      grid
      figure(3)
      if firsttime == 0,
		plot3(oldplCen(:,its1),oldplCen(:,its2),oldplCen(:,its3),'b+');
      end

      plot3(clCen(:,its1),clCen(:,its2),clCen(:,its3),'rx');
      oldplCen = clCen;
   end
   if count==maxcount,      % Plots at most maxcount changes
       dispersion(1:count-n_dist+1) = dispersion(n_dist:count);
       count = count-n_dist+1;
   end


   % termination conditions: small change in the centroids, 
   % or small change in MSE 

   if firsttime == 0,
      if chancent < prec,
         contin = 0;
         fprintf('\n Should stop, centroids');        
      end
      if max(abs(disper-oldisper))/max(abs(disper))<prec/100,
         contin =0;
         fprintf('\n Should stop: dispersion<!');        
      end
   else
      firsttime = 0;
      oldisper = disper;
   end

end % while contin

fprintf('\n Stop! \n');
figure(3)
hold;

figure(1)
clg;
asdfg = max(size(trSet));
plot3(trSet(1:perc:asdfg,its1),trSet(1:perc:asdfg,its2), ...
               trSet(1:perc:asdfg,its3),'.');
grid;
hold;
plot3(clCen(:,its1),clCen(:,its2),clCen(:,its3),'mx');
hold;

figure(3)