k-means聚类算法(matlab实现)

k-means聚类算法


前言

k-means聚类算法是硬聚类算法的一种,即在n纬欧几里得空间把n个样本数据分为k类。首先根据用户要确定聚类的数目k,随机性的选取k个样本,把每一个对象成为一个种子,每一个种子代表一个类的中心,对其余的每个对象,采用近邻原则,将它们赋给最近的类。重新计算在每个类中对象的均值新形成的聚类中心,重复进行这个过程直到函数式收敛为止。

k-means聚类算法的改进:
k-means++算法链接: https://blog.csdn.net/qq_40276082/article/details/130237784?spm=1001.2014.3001.5501


一、k-means聚类算法简介

K-means算法是一种无监督学习算法,用于将n个数据点划分为k个簇(cluster)。该算法的主要思想是将数据点划分为k个簇,使得同一簇内的数据点彼此相似,而簇间的数据点相差较大。
K-means算法是一种迭代算法,其核心是利用距离度量来度量不同数据点之间的相似度,并利用质心来表示簇的中心。在每次迭代中,算法通过更新簇质心和重新分配数据点来优化簇的划分。具体来说,算法的步骤如下:

初始化:随机选择k个数据点作为初始质心,将所有数据点分配到距离最近的质心所在的簇中。

簇分配:对于每个数据点,计算其到k个质心的距离,将其分配到距离最近的质心所在的簇中。

更新质心:对于每个簇,计算其所有数据点的平均值,将其作为新的质心。

重复步骤2和步骤3,直到簇的分配不再改变或达到最大迭代次数。

二、k-means聚类算法步骤

在这里插入图片描述

三、代码实现

代码如下(示例):

clc
clear
a=[22.5661	113.9808	66
22.6862	113.9405	65.5
22.5648	114.2446	75
22.5589	113.9507	65.5
22.559	114.2413	75
22.5628	113.9566	65.5
22.5	113.8957	66
22.5249	113.9309	65.5
22.5191	113.9358	65.5
22.4982	113.8985	66
23.0057	112.9708	75
23.0391	113.0726	66
22.9408	113.0627	72.5
23.0408	113.1116	65
23.0387	113.1074	65.5
23.0042	113.0704	66.5
23.0489	113.0719	66
23.0104	113.1463	66
23.0058	113.1356	65.5
22.8821	113.2179	68.5
23.0401	113.1075	65.5
22.9946	113.1409	72
22.9789	113.0167	72
22.9456	113.09	72
23.0441	113.1258	65.5
];
N=length(a(:,1));
X=[];
Y=[];

for j=1:N
    f=a(j,1);
w=a(j,2);
X=vertcat(X,f);
Y=vertcat(Y,w);

end
N=length(X(:,1));
plot(X, Y, 'k*'); % 绘出原始的数据点
xlabel('X');
ylabel('Y');
title('聚类之前的数据点');
n = 2; %将所有的数据点分为两类
m = 1; %迭代次数
eps = 1e-7; % 迭代结束的阈值
u1 = [X(1),Y(1)]; %初始化第一个聚类中心
u2 = [X(2),Y(2)]; %初始化第二个聚类中心
U1 = zeros(2,100);
U2 = zeros(2,100); %U1,U2 用于存放各次迭代两个聚类中心的横纵坐标
U1(:,2) = u1;
U2(:,2) = u2;
D = zeros(2,N); %初始化数据点与聚类中心的距离
while(abs(U1(1,m) - U1(1,m+1)) > eps || abs(U1(2,m) - U1(2,m+1) > eps || abs(U2(1,m) - U2(1,m+1)) > eps || abs(U2(2,m) - U2(2,m+1)) > eps))
    m = m +1;
    % 计算所有点到两个聚类中心的距离
for i = 1 : N
    D(1,i) = sqrt((X(i) - U1(1,m))^2 + (Y(i) - U1(2,m))^2);
end
for i = 1 : N
    D(2,i) = sqrt((X(i) - U2(1,m))^2 + (Y(i) - U2(2,m))^2);
end
A = zeros(2,N); % A用于存放第一类的数据点
B = zeros(2,N); % B用于存放第二类的数据点
for k = 1: N
    [MIN,index] = min(D(:,k)); 
    if index == 1  % 点属于第一个聚类中心
        A(1,k) = X(k);
        A(2,k) = Y(k);
    else           % 点属于第二个聚类中心
        B(1,k) = X(k);
        B(2,k) = Y(k);
    end
end
indexA = find(A(1,:) ~= 0); % 找出第一类中的点
indexB = find(B(1,:) ~= 0); % 找出第二类中的点
U1(1,m+1) = mean(A(1,indexA));
U1(2,m+1) = mean(A(2,indexA));
U2(1,m+1) = mean(B(1,indexB));
U2(2,m+1) = mean(B(2,indexB)); % 更新两个聚类中心
end
figure;
plot(A(1,indexA) , A(2,indexA), '*b'); % 作出第一类点的图形
hold on
plot(B(1,indexB) , B(2,indexB), '*g'); %作出第二类点的图形
hold on
centerx = [U1(1,m) U2(1,m)]
centery = [U1(2,m) U2(2,m)]
plot(centerx , centery, 'or'); % 画出两个聚类中心点
xlabel('X');
ylabel('Y');
title('聚类之后的数据点');
disp(['迭代的次数为:',num2str(m)]);

运行结果:
在这里插入图片描述


总结

K-means算法的优点在于简单、易于实现,适用于大型数据集。但是,它也有一些缺点。首先,K-means算法对于数据点分布的先验知识敏感,且容易收敛到局部最优解。其次,K-means算法要求每个簇的大小相等,这在某些情况下可能会导致不合理的结果。最后,K-means算法对于噪声和异常值敏感,因为它只考虑数据点之间的欧几里德距离,而没有考虑其他距离度量或数据结构。

function [idx, C, sumD, D] = kmeans(X, k, varargin) % varargin:实际输入参量 if nargin 1 % 大于1刚至少有一种距离 error(sprintf('Ambiguous ''distance'' parameter value: %s.', distance)); elseif isempty(i) % 如果是空的,则表明没有合适的距离 error(sprintf('Unknown ''distance'' parameter value: %s.', distance)); end % 针对不同的距离,处理不同 distance = distNames{i}; switch distance case 'cityblock' % sort 列元素按升序排列,Xord中存的是元素在原始矩阵中的列中对应的大小位置 [Xsort,Xord] = sort(X,1); case 'cosine' % 余弦 % 计算每一行的和的平方根 Xnorm = sqrt(sum(X.^2, 2)); if any(min(Xnorm) <= eps * max(Xnorm)) error(['Some points have small relative magnitudes, making them ', ... 'effectively zero.\nEither remove those points, or choose a ', ... 'distance other than ''cosine''.'], []); end % 标量化 Xnorm(:,ones(1,p))得到n*p的矩阵 X = X ./ Xnorm(:,ones(1,p)); case 'correlation' % 线性化 X = X - repmat(mean(X,2),1,p); % 计算每一行的和的平方根 Xnorm = sqrt(sum(X.^2, 2)); if any(min(Xnorm) 1 error(sprintf('Ambiguous ''start'' parameter value: %s.', start)); elseif isempty(i) error(sprintf('Unknown ''start'' parameter value: %s.', start)); elseif isempty(k) error('You must specify the number of clusters, K.'); end start = startNames{i}; % strcmp比较两个字符串 uniform是在X中随机选择K个点 if strcmp(start, 'uniform') if strcmp(distance, 'hamming') error('Hamming distance cannot be initialized with uniform random values.'); end % 标量化后的X Xmins = min(X,1); Xmaxs = max(X,1); end elseif isnumeric(start) % 判断输入是否是一个数 这里的start是一个K*P的矩阵,表示初始聚类中心 CC = start; % CC表示初始聚类中心 start = 'numeric'; if isempty(k) k = size(CC,1); elseif k ~= size(CC,1); error('The ''start'' matrix must have K rows.'); elseif size(CC,2) ~= p error('The ''start'' matrix must have the same number of columns as X.'); end if isempty(reps) reps = size(CC,3); elseif reps ~= size(CC,3); error('The third dimension of the ''start'' array must match the ''replicates'' parameter value.'); end % Need to center explicit starting points for 'correlation'. % 线性距离需要指定具体的开始点 % (Re)normalization for 'cosine'/'correlation' is done at each % iteration.每一次迭代进行“余弦和线性”距离正规化 % 判断是否相等 if isequal(distance, 'correlation') % repmat复制矩阵1*P*1 线性化 CC = CC - repmat(mean(CC,2),[1,p,1]); end else error('The ''start'' parameter value must be a string or a numeric matrix or array.'); end % ------------------------------------------------------------------ % 如果一个聚类丢失了所有成员,这个进程将被调用 if ischar(emptyact) emptyactNames = {'error','drop','singleton'}; i = strmatch(lower(emptyact), emptyactNames); if length(i) > 1 error(sprintf('Ambiguous ''emptyaction'' parameter value: %s.', emptyact)); elseif isempty(i) error(sprintf('Unknown ''emptyaction'' parameter value: %s.', emptyact)); end emptyact = emptyactNames{i}; else error('The ''emptyaction'' parameter value must be a string.'); end % ------------------------------------------------------------------ % 控制输出的显示示信息 if ischar(display) % strvcat 垂直连接字符串 i = strmatch(lower(display), strvcat('off','notify','final','iter')); if length(i) > 1 error(sprintf('Ambiguous ''display'' parameter value: %s.', display)); elseif isempty(i) error(sprintf('Unknown ''display'' parameter value: %s.', display)); end display = i-1; else error('The ''display'' parameter value must be a string.'); end % ------------------------------------------------------------------ % 输入参数K的控制 if k == 1 error('The number of clusters must be greater than 1.'); elseif n 2 % 'iter',\t 表示向后空出8个字符 disp(sprintf(' iter\t phase\t num\t sum')); end % ------------------------------------------------------------------ % Begin phase one: batch reassignments 第一队段:分批赋值 converged = false; iter = 0; while true % Compute the distance from every point to each cluster centroid % 计算每一个点到每一个聚类中心的距离,D中存放的是N*K个数值 D(:,changed) = distfun(X, C(changed,:), distance, iter); % Compute the total sum of distances for the current configuration. % Can't do it first time through, there's no configuration yet. % 计算当前配置的总距离,但第一次不能执行,因为还没有配置 if iter > 0 totsumD = sum(D((idx-1)*n + (1:n)')); % Test for a cycle: if objective is not decreased, back out % the last step and move on to the single update phase % 循环测试:如果目标没有减少,退出到最后一步,移动到第二阶段 % prevtotsumD表示前一次的总距离,如果前一次的距离比当前的小,则聚类为上一次的,即不发生变化 if prevtotsumD 2 % 'iter' disp(sprintf(dispfmt,iter,1,length(moved),totsumD)); end if iter >= maxit, % break(2) break; % 跳出while end end % Determine closest cluster for each point and reassign points to clusters % 决定每一个点的最近分簇,重新分配点到各个簇 previdx = idx; % 大小为n*1 % totsumD 被初始化为无穷大,这里表示总距离 prevtotsumD = totsumD; % 返回每一行中最小的元素,d的大小为n*1,nidx为最小元素在行中的位置,其大小为n*1,D为n*p [d, nidx] = min(D, [], 2); if iter == 0 % iter==0,表示第一次迭代 % Every point moved, every cluster will need an update % 每一个点需要移动,每一个簇更新 moved = 1:n; idx = nidx; changed = 1:k; else % Determine which points moved 决定哪一个点移动 % 找到上一次和当前最小元素不同的位置 moved = find(nidx ~= previdx); if length(moved) > 0 % Resolve ties in favor of not moving % 重新分配而不是移动 括号中是一个逻辑运算 确定需要移动点的位置 moved = moved(D((previdx(moved)-1)*n + moved) > d(moved)); end % 如果没有不同的,即当前的是最小元素,跳出循环,得到的元素已经是各行的最小值 if length(moved) == 0 % break(3) break; end idx(moved) = nidx(moved); % Find clusters that gained or lost members 找到获得的或者丢失的成员的分簇 % 得到idx(moved)和previdx(moved)中不重复出现的所有元素,并按升序排列 changed = unique([idx(moved); previdx(moved)])'; end % Calculate the new cluster centroids and counts. 计算新的分簇中心和计数 % C(changed,:)表示新的聚类中心,m(changed)表示聚类标号在idx中出现的次数 % sort 列元素按升序排列,Xsort存放对的元素,Xord中存的是元素在原始矩阵中的列中对应的大小位置 [C(changed,:), m(changed)] = gcentroids(X, idx, changed, distance, Xsort, Xord); iter = iter + 1; % Deal with clusters that have just lost all their members % 处理丢失所有成员的分簇,empties表示丢失元素的聚类标号 不用考虑 empties = changed(m(changed) == 0); if ~isempty(empties) switch emptyact case 'error' % 默认值,一般情况下不会出现 error(sprintf('Empty cluster created at iteration %d.',iter)); case 'drop' % Remove the empty cluster from any further processing % 移走空的聚类 D(:,empties) = NaN; changed = changed(m(changed) > 0); if display > 0 warning(sprintf('Empty cluster created at iteration %d.',iter)); end case 'singleton' if display > 0 warning(sprintf('Empty cluster created at iteration %d.',iter)); end for i = empties % Find the point furthest away from its current cluster. % Take that point out of its cluster and use it to create % a new singleton cluster to replace the empty one. % 找到离聚类中心最远距离的点,把该点从它的聚类中移走,用它来创建一个新的聚类,来代替空的聚类 % 得到列的最大元素(dlarge),以及该元素在列中的标号(lonely) [dlarge, lonely] = max(d); from = idx(lonely); % taking from this cluster 从当前聚类移走 C(i,:) = X(lonely,:); m(i) = 1; idx(lonely) = i; d(lonely) = 0; % Update clusters from which points are taken % 更新那些点被移走的聚类 [C(from,:), m(from)] = gcentroids(X, idx, from, distance, Xsort, Xord); changed = unique([changed from]); end end end end % phase one % ------------------------------------------------------------------ % Initialize some cluster information prior to phase two % 为第二阶段初始化一些先验聚类信息 针对特定的距离,默认的是欧氏距离 switch distance case 'cityblock' Xmid = zeros([k,p,2]); for i = 1:k if m(i) > 0 % Separate out sorted coords for points in i'th cluster, % and save values above and below median, component-wise % 分解出第i个聚类中挑选的点的坐标,保存它的上,下中位数 % reshape把矩阵分解为要求的行列数m*p % sort 列元素按升序排列,Xord中存的是元素在原始矩阵中的列中对应的大小位置 Xsorted = reshape(Xsort(idx(Xord)==i), m(i), p); % floor取比值小或者等于的最近的值 nn = floor(.5*m(i)); if mod(m(i),2) == 0 Xmid(i,:,1:2) = Xsorted([nn, nn+1],:)'; elseif m(i) > 1 Xmid(i,:,1:2) = Xsorted([nn, nn+2],:)'; else Xmid(i,:,1:2) = Xsorted([1, 1],:)'; end end end case 'hamming' Xsum = zeros(k,p); for i = 1:k if m(i) > 0 % Sum coords for points in i'th cluster, component-wise % 对基于分量对第i个聚类的坐标点求和 Xsum(i,:) = sum(X(idx==i,:), 1); end end end % ------------------------------------------------------------------ % Begin phase two: single reassignments 第二阶段:单独赋值 % m中保存的是每一个聚类的个数,元素和为n % find(m' > 0)得到m'中大于0的元素的位置(索引) % 实际情况(默认情况下)changed=1:k changed = find(m' > 0); lastmoved = 0; nummoved = 0; iter1 = iter; while iter < maxit % Calculate distances to each cluster from each point, and the % potential change in total sum of errors for adding or removing % each point from each cluster. Clusters that have not changed % membership need not be updated. % 计算从每一个点到每一个聚类的距离,潜在由于总距离发生变化移除或添加引起的错误。 % 那些成员没有改变的聚类不需要更新。 % % Singleton clusters are a special case for the sum of dists % calculation. Removing their only point is never best, so the % reassignment criterion had better guarantee that a singleton % point will stay in its own cluster. Happily, we get % Del(i,idx(i)) == 0 automatically for them. % 单独聚类在计算距离时是一个特殊情况,仅仅移除点不是最好的方法,因此重新赋值准则能够保证一个 % 单独的点能够留在它的聚类中,可喜的是,自动有Del(i,idx(i)) == 0。 switch distance case 'sqeuclidean' for i = changed % idx存放的距离矩阵D中每一行的最小元素所处的列,也即位置 mbrs = (idx == i); sgn = 1 - 2*mbrs; % -1 for members, 1 for nonmembers % 表示只有一个聚类 if m(i) == 1 % prevent divide-by-zero for singleton mbrs 防止单个成员做0处理 sgn(mbrs) = 0; end Del(:,i) = (m(i) ./ (m(i) + sgn)) .* sum((X - C(repmat(i,n,1),:)).^2, 2); end case 'cityblock' for i = changed if mod(m(i),2) == 0 % this will never catch singleton clusters ldist = Xmid(repmat(i,n,1),:,1) - X; rdist = X - Xmid(repmat(i,n,1),:,2); mbrs = (idx == i); sgn = repmat(1-2*mbrs, 1, p); % -1 for members, 1 for nonmembers Del(:,i) = sum(max(0, max(sgn.*rdist, sgn.*ldist)), 2); else Del(:,i) = sum(abs(X - C(repmat(i,n,1),:)), 2); end end case {'cosine','correlation'} % The points are normalized, centroids are not, so normalize them normC(changed) = sqrt(sum(C(changed,:).^2, 2)); if any(normC 0 % Resolve ties in favor of not moving % 重新分配而不是移动 确定移动的位置 moved = moved(Del((previdx(moved)-1)*n + moved) > minDel(moved)); end if length(moved) == 0 % Count an iteration if phase 2 did nothing at all, or if we're % in the middle of a pass through all the points if (iter - iter1) == 0 | nummoved > 0 iter = iter + 1; if display > 2 % 'iter' disp(sprintf(dispfmt,iter,2,nummoved,totsumD)); end end converged = true; break; end % Pick the next move in cyclic order moved = mod(min(mod(moved - lastmoved - 1, n) + lastmoved), n) + 1; % If we've gone once through all the points, that's an iteration if moved 2 % 'iter' disp(sprintf(dispfmt,iter,2,nummoved,totsumD)); end if iter >= maxit, break; end nummoved = 0; end nummoved = nummoved + 1; lastmoved = moved; oidx = idx(moved); nidx = nidx(moved); totsumD = totsumD + Del(moved,nidx) - Del(moved,oidx); % Update the cluster index vector, and rhe old and new cluster % counts and centroids idx(moved) = nidx; m(nidx) = m(nidx) + 1; m(oidx) = m(oidx) - 1; switch distance case 'sqeuclidean' C(nidx,:) = C(nidx,:) + (X(moved,:) - C(nidx,:)) / m(nidx); C(oidx,:) = C(oidx,:) - (X(moved,:) - C(oidx,:)) / m(oidx); case 'cityblock' for i = [oidx nidx] % Separate out sorted coords for points in each cluster. % New centroid is the coord median, save values above and % below median. All done component-wise. Xsorted = reshape(Xsort(idx(Xord)==i), m(i), p); nn = floor(.5*m(i)); if mod(m(i),2) == 0 C(i,:) = .5 * (Xsorted(nn,:) + Xsorted(nn+1,:)); Xmid(i,:,1:2) = Xsorted([nn, nn+1],:)'; else C(i,:) = Xsorted(nn+1,:); if m(i) > 1 Xmid(i,:,1:2) = Xsorted([nn, nn+2],:)'; else Xmid(i,:,1:2) = Xsorted([1, 1],:)'; end end end case {'cosine','correlation'} C(nidx,:) = C(nidx,:) + (X(moved,:) - C(nidx,:)) / m(nidx); C(oidx,:) = C(oidx,:) - (X(moved,:) - C(oidx,:)) / m(oidx); case 'hamming' % Update summed coords for points in each cluster. New % centroid is the coord median. All done component-wise. Xsum(nidx,:) = Xsum(nidx,:) + X(moved,:); Xsum(oidx,:) = Xsum(oidx,:) - X(moved,:); C(nidx,:) = .5*sign(2*Xsum(nidx,:) - m(nidx)) + .5; C(oidx,:) = .5*sign(2*Xsum(oidx,:) - m(oidx)) + .5; end changed = sort([oidx nidx]); end % phase two % ------------------------------------------------------------------ if (~converged) & (display > 0) warning(sprintf('Failed to converge in %d iterations.', maxit)); end % Calculate cluster-wise sums of distances nonempties = find(m(:)'>0); D(:,nonempties) = distfun(X, C(nonempties,:), distance, iter); d = D((idx-1)*n + (1:n)'); sumD = zeros(k,1); for i = 1:k sumD(i) = sum(d(idx == i)); end if display > 1 % 'final' or 'iter' disp(sprintf('%d iterations, total sum of distances = %g',iter,totsumD)); end % Save the best solution so far if totsumD 3 Dbest = D; end end end % Return the best solution
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

一寸光阴不可轻

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值