function [centr,clst,sse] = kmeans_st(X,k,restarts)
% KMEANS: K-means clustering of n points into k clusters so that the
% within-cluster sum of squares is minimized. Based on Algorithm
% AS136, which seeks a local optimum such that no movement of a
% point from one cluster to another will reduced the within-cluster
% sum of squares. Tends to find spherical clusters.
% Not globally optimal against all partitions, which is an NP-complete
% problem; thus allows for multiple restarts.
%
% Usage: [centr,clst,sse] = kmeans(X,k,restarts)
%
% X = [n x p] data matrix.
% k = number of desired clusters.
% restarts = optional number of restarts, after finding a new minimum
% sse, needed to end search [default=0].
% --------------------------------------------------------------------
% centr = [k x p] matrix of cluster centroids.
% clst = [n x 1] vector of cluster memberships.
% sse = total within-cluster sum of squared deviations.
%
% Hartigan,JA & MA Wong. 1979. Algorithm AS 136: A k-means clustering
% algorithm. Appl. Stat. 200:100-108.
% Milligan,G.W. & L. Sokol. 1980. A two-stage clustering algorithm with
% robustness recovery characteristics. Educational and Psychological
% Measurement 40:755-759.
% RE Strauss, 8/26/98
% 8/21/99 - changed misc statements for Matlab v5.
% 10/14/00 - use MEANS rather than MEAN to update cluster centers;
% remove references to TRUE and FALSE.
if (nargin<3) restarts = []; end
if (isempty(restarts))
restarts = 0;
end
[n,p] = size(X);
if (k<1 || k>n)
error('KMEANS: k out of range 1-N');
end
max_iter = 50;
highval = 10e8;
upgma_flag = 1; % Do UPGMA first time thru
% Repeat optimization until have 'restart' attempts with no better solution
restart_iter = restarts + 1;
best_sse = highval;
while (restart_iter > 0)
restart_iter = restart_iter - 1;
c1 = zeros(n,1);
c2 = c1;
clst = c1;
nclst = zeros(k,1);
% Estimate initial cluster centers via UPGMA the first time, and randomly
% thereafter.
if (upgma_flag) % If UPGMA,
upgma_flag = 0;
dist = eucl(X); % All interpoint distances
topo = upgma(dist,[],0); % UPGMA dendrogram topology
grp = [1:n]'; % Identify obs in each of k grps
for step = 1:(n-k)
t1 = topo(step,1);
t2 = topo(step,2);
t3 = topo(step,3);
indx = find(grp==t1);
grp(indx) = t3 * ones(length(indx),1);
indx = find(grp==t2);
grp(indx) = t3 * ones(length(indx),1);
end
grpid = uniquef(grp); % Unique group identifiers
centr = zeros(k,p);
for kk = 1:length(grpid); % Calculate group centroids
g = grpid(kk);
indx = find(grp == g);
Xk = X(indx,:);
if (size(Xk,1)>1)
centr(kk,:) = mean(Xk);
else
centr(kk,:) = Xk;
end
end
else % Else if randomized,
rndprm = randperm(n); % Pull out random set of points
centr = X(rndprm(1:k),:);
end
% For each point i, find its two closest centers, c1(i) and c2(i), and
% assign it to c1(i)
dist = zeros(n,k);
for kk = 1:k % Dists to centers of all clusters
dist(:,kk) = eucl(X,centr(kk,:));
end
[d,indx] = sort(dist'); % Sort points separately
c1 = indx(1,:)';
c2 = indx(2,:)';
% Update cluster centers to be the centroids of points contained within them
for kk = 1:k
indx = find(c1==kk);
len_indx = length(indx);
% if (len_indx>1)
% centr(kk,:) = mean(X(indx,:));
% else
% centr(kk,:) = X(indx,:);
% end
centr(kk,:) = means(X(indx,:));
nclst(kk) = len_indx;
end
% Initialize working matrices
live = ones(k,1);
R1 = zeros(n,1);
R2 = zeros(k,1);
adj1 = nclst ./ (nclst+1);
adj2 = nclst ./ (nclst-1+eps);
update = n * ones(k,1);
for i = 1:n % Adjusted dist from each pt to own cluster center
R1(i) = adj2(c1(i)) * eucl(X(i,:),centr(c1(i),:))^2;
end
% Iterate
iter = 0;
while ((iter < max_iter) & (any(live)))
iter = iter+1;
% Optimal-transfer stage
i = 0;
while ((i<n) & any(live)) % Cycle thru points
i = i+1;
update = update-1; % Decrement cluster update counters
L1 = c1(i); % Pt i is in cluster L1
if (nclst(L1)>1) % If point i is not sole member of L1,
R2 = adj1(L1) .* eucl(centr,X(i,:)).^2; % Dists from pt to cluster centers
if (live(L1)) % If L1 is in live set,
R2(L1) = highval; % find min R2 over all clusters
[r2min,L2] = min(R2);
else % If L1 is not in live set,
indx = find(~live); % find min R2 over live clusters only
R2([L1;indx]) = highval * ones(length(indx)+1,1);
[r2min,L2] = min(R2);
end
if (r2min >= R1(i)) % No reallocation necessary
c2(i) = L2;
else % Else reallocate pt to new cluster
c1(i) = L2;
c2(i) = L1;
for kk = [L1,L2] % Recalc cluster centers and sizes
indx = find(c1==kk);
len_indx = length(indx);
if (len_indx>1)
centr(kk,:) = mean(X(indx,:));
else
centr(kk,:) = X(indx,:);
end
nclst(kk) = len_indx;
end
adj1([L1,L2]) = nclst([L1,L2]) ./ (nclst([L1,L2])+1);
adj2([L1,L2]) = nclst([L1,L2]) ./ (nclst([L1,L2])-1+eps);
R1(i) = adj2(c1(i)) * eucl(X(i,:),centr(c1(i),:))^2;
live([L1,L2]) = [1,1]; % Put into live set
update([L1,L2]) = [n,n]; % Reinitialize update counters
end
end; % if nclst(L1)>1
indx = find(update<1); % Remove clusters from live set
if (~isempty(indx)) % that haven't been recently updated
live(indx) = zeros(length(indx),1);
end
end; % while (i<n & any(live))
if (any(live))
% Quick transfer stage
for i = 1:n % Cycle thru points
L1 = c1(i);
L2 = c2(i);
if (any(live([L1,L2])) & nclst(L1)>1)
R2 = adj1(L2) .* eucl(centr(L2,:),X(i,:)).^2; % Dists from pt to L2 center
if (R1>=R2)
temp = c1(i); % Switch L1 & L2
c1(i) = c2(i);
c2(i) = temp;
for kk = [L1,L2] % Recalc cluster centers and sizes
indx = find(c1==kk);
centr(kk,:) = mean(X(indx,:));
nclst(kk) = length(indx);
end
adj1([L1,L2]) = nclst([L1,L2]) ./ (nclst([L1,L2])+1);
adj2([L1,L2]) = nclst([L1,L2]) ./ (nclst([L1,L2])-1+eps);
R1(i) = adj2(c1(i)) * eucl(X(i,:),centr(c1(i),:))^2;
live([L1,L2]) = [1,1]; % Put into live set
update([L1,L2]) = [n,n]; % Reinitialize update counters
end
end
end
end; % if any(live)
end; % while iter<max_iter & any(live)
sse = 0; % Calc final total sum-of-squares
for i=1:n
sse = sse + eucl(X(i,:),centr(c1(i),:))^2;
end
if (sse < best_sse) % Update best solution so far
best_sse = sse;
best_c1 = c1;
best_centr = centr;
restart_iter = restarts;
end
end; % while (restart_iter > 0)
clst = best_c1;
centr = best_centr;
sse = best_sse;
return;
Datasets
Models
