cluster2b_ser.m

動態聚類k-means演算 將輸入在程式中的數據資料 給予適當的分群

function [dbelong, neach, clusters, niter, ST, DB] = ...
    cluster2_ser(data, nc, np, incnp, eta, deta, maxEpoch, info, psec);
%
% The size of the dataset
indata=data(:,2:97);
target=data(:,1);
[ndata,nvar] = size(indata);
class=ones(1,ndata);
col=class;
leaf=0;
if info >= 1
  disp(sprintf('Number of datapoints: %d', ndata));
  disp(sprintf('Number of variables:  %d', nvar));
end

% Normalize data to unit standard deviation and zero mean
% using prestd.
[indata, datamean, datastd] = kmeans_norm(indata); 

%%%%%%%%%% Parameters for k-means %%%%%%%%%%

    npoints = ndata;

    clust(1,1)=0;
    clust(1,2)=0;

CP=zeros(1,200);
clust=zeros(1000,nvar);
dotter=1000*ones(1000,2);
mor=dotter;
N=777*ones(1,500);
N(1)=200;
ant=1;
cl=0;
go2=1;
while leaf<63 & cl<900
    cl=cl+1;
       
    % The starting clusters
    u=find(class==cl);
    
    if length(u)>=3
        rand('state',sum(100*clock));
        clusters(1,:) = 0.0001 * (rand(1,nvar) - 0.5)+clust(cl,:);   % Slumpar startpunkter.
        clusters(2,:) = 0.0001 * (rand(1,nvar) - 0.5)+clust(cl,:);   % Slumpar startpunkter.
        bclusters=clusters;
        niter = 0;
        go = 1;
        qerr_old = 0.0;
        cold = zeros(nc,nvar);
      
    
        while go == 1

            niter = niter + 1;

            data = indata;
  
            % Compute the distance between each datapoint to the clusters
            dists = kmeans_dist(indata, clusters);
    
            % Assign datapoints to clusters and compute updates
            [minval, minidx] = min(dists');
            for i=1:nc

                idxtmp = find(minidx == i & class==cl); 
                if length(idxtmp > 0)
                    
                    cpoints = indata(idxtmp,:);
	
	                % Compute the mean of the indatapoints belonging to cluster i
	                if length(idxtmp)==1
                        centroid=cpoints;
                    else
                        centroid = mean(cpoints);
                    end
                    update(i,:) = centroid - clusters(i,:);
                    
                
                else
	                update(i,:) = zeros(1,2);
                    
                    
                    
                    
                end
            end
            
            % Update the cluster centers
            clusters = clusters + eta * update;    
            update_old=update;
            
            
            % Compute changes
            dcluster = sum( sqrt( sum( ( (clusters-cold).*(clusters-cold) )') ) ) / ...
            nc;
            cold = clusters;
            %dcluster = sum( sqrt( sum( (update.*update)') ) ) / nc;

            % Compute some clusters statistics
            [dbelong, neach, avedist, quanterr] = ...
            cluster_stat(indata, clusters, nc,cl,class);
  
            % Some plotting if 2 dimenstions
            if nvar == 2 & info > 0 & cl<8  
                plot(clusters(:,1), clusters(:,2),'+r');
                %cluster_plotclusters(fig, clusters,cl,clust);

            end
  
            % Stopping criteria
            if abs(quanterr - qerr_old) < 0.00001 | niter >= maxEpoch 
                go = 0;
            end
            qerr_old = quanterr;
  
            % Update the learning rate
            %eta = eta * deta;
  
            % This for visualization purposes
            %pause(psec);
        end

        
        % Final plotting if in 2 dimensions
        
        for i=1:2
            idxtmp = find(minidx == i & class==cl);
            ant=ant+1;
            class(idxtmp)=ant;
            N(ant)=length(idxtmp);  
            dotter(cl,i)=ant;    
            mor(ant)=cl;
        end

            
        
        clust(ant-1:ant,:)=clusters;

        
    elseif length(u)==2
        
        ant=ant+1;
        clust(ant,:)=indata(u(1),:);
        dotter(cl,1)=ant;
        mor(ant)=cl;
        class(u(1))=ant;
        N(ant)=1;
        ant=ant+1;
        clust(ant,:)=indata(u(2),:);
        dotter(cl,2)=ant;
        mor(ant)=cl;
        class(u(2))=ant;
        N(ant)=1;
        
    elseif length(u)==1
        CP(u)=2;
        leaf =leaf+1;
        
        floor(100*leaf/63)
        
    end

end




for i=1:4
    tmp=find(target==i);
    j=class(tmp);
    ccol(j)=i;
end
                

g=ones(1,1000);
y=zeros(1,1000);
x=y;
for i=1:cl
    d=dotter(i,:);
    if d(1)<1000
        
        g(d(1))=g(i)+1;
        x(d(1))=x(i)+.5^g(i);
        g(d(2))=g(i)+1;
        x(d(2))=x(i)-.5^g(i);

    end
end




u=find(N(1:cl)==1);

[tmp, q]=sort(x(u));
z=u(q);

y(z)=1:length(z);



[tmp, j]=sort(N(1:cl));




figure
hold on
for i=1:length(j)
    s=j(i);
    d=dotter(s,:);
    if d(1)<1000
        y(s)=(y(d(1))+y(d(2)))/2;
        plot ([g(d(1)), g(s), g(s), g(d(2))], [y(d(1)), y(d(1)), y(d(2)), y(d(2))],'k')

        if ccol(d(1)) == ccol(d(2))
            ccol(s)=ccol(d(1));
        end
    end
end
    plot ([g(1)-1,g(1)],[y(1),y(1)],'k')
    

for i=1:cl
    if ccol(i)==1
        plot ([g(i)-1,g(i)],[y(i),y(i)],'m')
    end
    if ccol(i)==2
        plot ([g(i)-1,g(i)],[y(i),y(i)],'g')
    end
    if ccol(i)==3
        plot ([g(i)-1,g(i)],[y(i),y(i)],'r')
    end
    if ccol(i)==4
        plot ([g(i)-1,g(i)],[y(i),y(i)],'b')
    end
end
axis ([0 max(g) -5 70])