Kmeans聚类算法及其matlab源码

本文介绍了K-means聚类算法，并注释了部分matlab实现的源码。

K-means算法

K-means算法是一种硬聚类算法，根据数据到聚类中心的某种距离来作为判别该数据所属类别。K-means算法以距离作为相似度测度。

假设将对象数据集分为个不同的类，k均值聚类算法步骤如下：

Step1：随机从对象集中抽取个对象作为初始聚类中心；

Step2：对于所有的对象，分别计算其到各个聚类中的欧氏距离，相互比较后将其归属于距离最小的那一类；

Step3：根据step2得到的初始分类，对每个类别计算均值用来更新聚类中心；

Step4：根据新的聚类中心，重复进行step2和step3,直至满足算法终止条件。

K-means算法是基于划分的思想，因此算法易于理解且实现方法简单易行，但需要人工选择初始的聚类数目即算法是带参数的。类的数目确定往往非常复杂和具有不确定性，因此需要专业的知识和行业经验才能较好的确定。而且因为初始聚类中心的选择是随机的，因此会造成部分初始聚类中心相似或者处于数据边缘，造成算法的迭代次数明显增加，甚至会因为个别数据而造成聚类失败的现象。

其流程图大致如下：

matlab源码

function varargout = kmeans(X, k, varargin)
%K均值聚类.
%   IDX = KMEANS(X, K) 分割X[N P]的数据矩阵中的样本为K个类，是一种最小化类内点到中心距离和的总和的分割。
%   矩阵X中的行对应的是数据样本，列对应的是变量。
%   提示: 当X是一个向量，本函数会忽略它的方向，将其当作一个[N 1]的数据矩阵。
%   KMEANS 函数返回一个代表各个数据样本所属类别索引的[N 1]维向量，函数默认使用平方的欧氏距离。
%   KMEANS 将NaNs当作丢失的数据并且忽略X中任何包含NaNs的行
%
%
%   [IDX, C] = KMEANS(X, K) 返回一个包含K个聚类中心的[K P]维的矩阵C.
%
%   [IDX, C, SUMD] = KMEANS(X, K) 返回一个类间点到聚类中心距离和的[K 1]维向量SUMD。
%
%   [IDX, C, SUMD, D] = KMEANS(X, K) 返回一个每个点到任一聚类中心距离的[N K]维矩阵D。
%
%   [ ... ] = KMEANS(..., 'PARAM1',val1, 'PARAM2',val2, ...) 指定了可选参数对（参数名/参数值）来控制算法的迭代。
%   参数如下：
%
%   'Distance' - 距离测度,  P维空间,  KMEANS算法需要最小化的值
%         可以选择:
%          'sqeuclidean'  - 平方的欧氏距离 (默认)
%          'cityblock'    - 曼哈顿距离，各维度差异的绝对值之和。
%          'cosine'       - 1减去两个样本（当作向量）夹角的余弦值
%          'correlation'  - 1减去两个样本(当作值的序列)的相关系数
%
%          'hamming'      - 汉明距离，二进制数据相匹配位置的不同比特百分比。
%
%   'Start' - 选择初始聚类中心的方法,有时候也称作种子。
%         可以选择:
%          'plus'    - 默认值。 利用k-means++算法从X中选择K个观测值：从X中随机的选取第一个聚类中心；之后的
%                      聚类中心以一定的概率从剩余的样本中根据其到最近的聚类中心的比例来随机的选取。
%          'sample'  - 随机的从X中选取K个观测值。
%          'uniform' - 根据X的取值范围均匀的随机选取K个样本，对汉明距离不适用。
%          'cluster' - 随机的利用X中10%的样本进行一个预聚类的阶段，预聚类阶段的初始聚类中心选取采用‘sample’。
%           matrix   - 一个初始聚类中心的[K P]维矩阵。此时,你可以用[]代替K，算法会自动的根据矩阵的第一个维度推算K值。
%                      你也可以使用3D数组，暗含着第三维为参数'Replicates'的值。
%
%   'Replicates' - 重复聚类的次数,默认为1。 每次都会有一个新的初始聚类中心。
%
%   'EmptyAction' - 发生空类时的处理措施。
%         可以选择:
%          'singleton' - 默认方法。利用据该中心最远的一个观测值建立一个新的类。
%          'error'     - 将产生空类作为一个错误(error)。
%          'drop'      - 移除空类并将对应的C和D中的值设置为NaN。
%
%
%   'Options' - 迭代算法最小化拟合准则(？)的选项，通过STATSET创建。 Choices of STATSET
%          STATSET参数可以选择:
%
%          'Display'  - 显示输出的哪一阶段的值，可以为 'off'(默认），‘iter’和‘final’;
%          'MaxIter'  - 最大的迭代次数，默认值为100。
%
%          'UseParallel'  - 在满足条件下，如果为真则开启并行计算否则使用串行模式。默认使用串行模式。
%          'UseSubstreams'  - 默认不使用。
%          'Streams'  - 这些区域指明是否执行并行的多个‘Start’值和当产生初始聚类中心时如何使用随机数值，
%                       更详细的参考 PARALLELSTATS。
%                       提示: 如果 'UseParallel'为TRUE且 'UseSubstreams'为FALSE,
%                       那么'Streams'的长度必须等于KMEANS使用的workers的数目。
%                       如果打开了并行池，那么它的大小和并行池一样。如果没有打开并行池，
%                       那么MATLAB可能会自动的打开(这取决于你的安装设置)。为了得到更好的结果，
%                       建议运用PARPOOL命令创建并行池的优先级以便当'UseParallel'为TRUE时执行算法。
%
%         'OnlinePhase' - 标志位，表示KMEANS是否除了运行一个"batch update"阶段还需一个"on-line
%                         update"阶段 。on-line阶段在大数据量时耗时很多。默认为‘off’。
%
%   示例:
%
%       X = [randn(20,2)+ones(20,2); randn(20,2)-ones(20,2)];
%       opts = statset('Display','final');
%       [cidx, ctrs] = kmeans(X, 2, 'Distance','city', ...
%                             'Replicates',5, 'Options',opts);
%       plot(X(cidx==1,1),X(cidx==1,2),'r.', ...
%            X(cidx==2,1),X(cidx==2,2),'b.', ctrs(:,1),ctrs(:,2),'kx');
%
%   也可以参考LINKAGE, CLUSTERDATA, SILHOUETTE。%   KMEANS 运用两阶段迭代算法来最小化K个类中样本到中心的距离和。
%   第一阶段利用文献中经常描述的"batch" 更新, 其中每次迭代中都一
%   次性地将样本分配到最近的聚类中心，然后更新聚类中心。这一阶段
%   偶尔(特别实在小样本的时候)会陷入局部最优。因此，"batch"阶段可
%   以考虑为第二阶段提供一个快速且可能为近似解的初始聚类中心。第二
%   阶段利用文献中常提及的"on-line"更新， 其中。如果能够减小距离
%   的总和那么其中的样本点都是单独地重新分配且每次分配后都重新计算
%   聚类中心。第二阶段中的每次迭代都会遍历所有的点，但是on-line阶段会收
%   敛到一个局部最小值。寻找全局最优的问题一般只能通过详细(幸运)地选择初始
%   聚类中心，但是使用重复多次的使用随机初始聚类中心中的典型结果是一个全局最小。
%
%  参考文献:
%
%   [1] Seber, G.A.F. (1984) Multivariate Observations, Wiley, New York.
%   [2] Spath, H. (1985) Cluster Dissection and Analysis: Theory, FORTRAN
%       Programs, Examples, translated by J. Goldschmidt, Halsted Press,
%       New York.%判断输入变量是否少于两个
if nargin < 2 error(message('stats:kmeans:TooFewInputs'));
end
%判断X是否是实数矩阵；
if ~isreal(X) error(message('stats:kmeans:ComplexData'));
end
%查找是否有NaN数据，有的话就删除，更新X矩阵；
wasnan = any(isnan(X),2);
hadNaNs = any(wasnan);
if hadNaNswarning(message('stats:kmeans:MissingDataRemoved'));X = X(~wasnan,:);
end% 获取X矩阵的维数
[n, p] = size(X);
%参数名与默认参数值设置
pnames = {   'distance'  'start' 'replicates' 'emptyaction' 'onlinephase' 'options' 'maxiter' 'display'};
dflts =  {'sqeuclidean' 'plus'          []  'singleton'         'off'        []        []        []};
[distance,start,reps,emptyact,online,options,maxit,display] ...= internal.stats.parseArgs(pnames, dflts, varargin{:});distNames = {'sqeuclidean','cityblock','cosine','correlation','hamming'};
distance = internal.stats.getParamVal(distance,distNames,'''Distance''');switch distancecase 'cosine'Xnorm = sqrt(sum(X.^2, 2));%模长if any(min(Xnorm) <= eps(max(Xnorm)))error(message('stats:kmeans:ZeroDataForCos'));endX =  bsxfun(@rdivide,X,Xnorm);%标准化case 'correlation'X = bsxfun(@minus, X, mean(X,2));Xnorm = sqrt(sum(X.^2, 2));if any(min(Xnorm) <= eps(max(Xnorm)))error(message('stats:kmeans:ConstantDataForCorr'));endX =  bsxfun(@rdivide,X,Xnorm);case 'hamming'if  ~all( X(:) ==0 | X(:)==1)error(message('stats:kmeans:NonbinaryDataForHamm'));end
endXmins = [];
Xmaxs = [];
CC = [];
if ischar(start)startNames = {'uniform','sample','cluster','plus','kmeans++'};j = find(strncmpi(start,startNames,length(start)));if length(j) > 1error(message('stats:kmeans:AmbiguousStart', start));elseif isempty(j)error(message('stats:kmeans:UnknownStart', start));elseif isempty(k)error(message('stats:kmeans:MissingK'));endstart = startNames{j};if strcmp(start, 'uniform')if strcmp(distance, 'hamming')error(message('stats:kmeans:UniformStartForHamm'));endXmins = min(X,[],1);%求每一列的最小值Xmaxs = max(X,[],1);%求每一列的最大值end
elseif isnumeric(start) %如果初始中心是数值类型(numeric)CC = start;start = 'numeric';if isempty(k)k = size(CC,1);%如果K为空通过数值的初始聚类中心获取K值elseif k ~= size(CC,1);%检测初始聚类中心行是否合法error(message('stats:kmeans:StartBadRowSize'));elseif size(CC,2) ~= p %检测初始聚类中心列是否合法error(message('stats:kmeans:StartBadColumnSize'));endif isempty(reps) reps = size(CC,3);%如果重复次数参数为空，检测初始聚类中心的第三维获取elseif reps ~= size(CC,3);error(message('stats:kmeans:StartBadThirdDimSize'));end% Need to center explicit starting points for 'correlation'. (Re)normalization% for 'cosine'/'correlation' is done at each iteration.if isequal(distance, 'correlation')CC = bsxfun(@minus, CC, mean(CC,2));%如果距离测度为相关性需要中心化数据end
elseerror(message('stats:kmeans:InvalidStart'));
endemptyactNames = {'error','drop','singleton'};
emptyact = internal.stats.getParamVal(emptyact,emptyactNames,'''EmptyAction''');[~,online] = internal.stats.getParamVal(online,{'on','off'},'''OnlinePhase''');
online = (online==1);% 'maxiter' and 'display' are grandfathered as separate param name/value pairs
if ~isempty(display)options = statset(options,'Display',display);
end
if ~isempty(maxit)options = statset(options,'MaxIter',maxit);
endoptions = statset(statset('kmeans'), options);
display = find(strncmpi(options.Display, {'off','notify','final','iter'},...length(options.Display))) - 1;
maxit = options.MaxIter;%确定最大迭代次数if ~(isscalar(k) && isnumeric(k) && isreal(k) && k > 0 && (round(k)==k))error(message('stats:kmeans:InvalidK'));% elseif k == 1% this special case works automatically
elseif n < kerror(message('stats:kmeans:TooManyClusters'));
end% Assume one replicate 检测重复次数的值
if isempty(reps)reps = 1;
elseif ~internal.stats.isScalarInt(reps,1)error(message('stats:kmeans:BadReps'));
end[useParallel, RNGscheme, poolsz] = ...internal.stats.parallel.processParallelAndStreamOptions(options,true);usePool = useParallel && poolsz>0;%检测是否使用并行池% Prepare for in-progress
if display > 1 % 'iter' or 'final'if usePool% If we are running on a parallel pool, each worker will generate% a separate periodic report.  Before starting the loop, we% seed the parallel pool so that each worker will have an% identifying label (eg, index) for its report.internal.stats.parallel.distributeToPool( ...'workerID', num2cell(1:poolsz) );% Periodic reports behave differently in parallel than they do% in serial computation (which is the baseline).% We advise the user of the difference.if display == 3 % 'iter' onlywarning(message('stats:kmeans:displayParallel2'));fprintf('    worker\t  iter\t phase\t     num\t         sum\n' );endelseif useParallelwarning(message('stats:kmeans:displayParallel'));endif display == 3 % 'iter' onlyfprintf('  iter\t phase\t     num\t         sum\n');endend
endif issparse(X) || ~isfloat(X) || strcmp(distance,'cityblock') || ...strcmp(distance,'hamming')[varargout{1:nargout}] = kmeans2(X,k, distance, emptyact,reps,start,...Xmins,Xmaxs,CC,online,display, maxit,useParallel, RNGscheme,usePool,...wasnan,hadNaNs,varargin{:});return;
endemptyErrCnt = 0;% Define the function that will perform one iteration of the
% loop inside smartFor
loopbody = @loopBody;%定义循环体函数% Initialize nested variables so they will not appear to be functions here
%初始化循环嵌套变量
totsumD = 0;
iter = 0;%将数据转置
X = X';
Xmins = Xmins';
Xmaxs = Xmaxs';% 执行KMEANS多次(reps)在各自的工作区上.
ClusterBest = internal.stats.parallel.smartForReduce(...reps, loopbody, useParallel, RNGscheme, 'argmin');% 选出最优解
varargout{1} = ClusterBest{5};%最优解的索引idx
varargout{2} = ClusterBest{6}';%最优解的聚类中心C
varargout{3} = ClusterBest{3}; %最优解的类内距离和sumD
totsumDbest = ClusterBest{1};%最优解的所有类内距离和的总和if nargout > 3varargout{4} = ClusterBest{7}; %最优解的点到任意聚类中心的距离
endif display > 1 % 'final' or 'iter'fprintf('%s\n',getString(message('stats:kmeans:FinalSumOfDistances',sprintf('%g',totsumDbest))));
endif hadNaNsvarargout{1} = statinsertnan(wasnan, varargout{1});% idxbest if nargout > 3varargout{4} = statinsertnan(wasnan, varargout{4}); %Dbestend
endfunction cellout = loopBody(rep,S)%循环体函数if isempty(S)S = RandStream.getGlobalStream;endif display > 1 % 'iter'if usePooldispfmt = '%8d\t%6d\t%6d\t%8d\t%12g\n';labindx = internal.stats.parallel.workerGetValue('workerID');elsedispfmt = '%6d\t%6d\t%8d\t%12g\n';endend%定义元胞数组cellout = cell(7,1);  % cellout{1}类间距离总和% cellout{2}重复次数% cellout{3}类内距离总和% cellout{4}迭代次数% cellout{5}索引% cellout{6}聚类中心% cellout{7}距离% Populating total sum of distances to Inf. This is used in the% reduce operation if update fails due to empty cluster.cellout{1} = Inf;%赋值cellout{2} = rep;%初始化聚类中心switch startcase 'uniform'%C = Xmins(:,ones(1,k)) + rand(S,[p,k]).*(Xmaxs(:,ones(1,k))-Xmins(:,ones(1,k)));C = Xmins(:,ones(1,k)) + rand(S,[k,p])'.*(Xmaxs(:,ones(1,k))-Xmins(:,ones(1,k)));% For 'cosine' and 'correlation', these are uniform inside a subset% of the unit hypersphere.仍需要为'correlation'进行中心化.  %  'cosine'/'correlation'的正交化在每次迭代中完成if isequal(distance, 'correlation')C = bsxfun(@minus, C, mean(C,1));endif isa(X,'single')C = single(C);endcase 'sample'C = X(:,randsample(S,n,k));case 'cluster'Xsubset = X(:,randsample(S,n,floor(.1*n)));% Turn display off for the initializationoptIndex = find(strcmpi('options',varargin));if isempty(optIndex)opts = statset('Display','off');varargin = [varargin,'options',opts];elsevarargin{optIndex+1}.Display = 'off';end[~, C] = kmeans(Xsubset', k, varargin{:}, 'start','sample', 'replicates',1);C = C';case 'numeric'C = CC(:,:,rep)';if isa(X,'single')C = single(C);endcase {'plus','kmeans++'}% Select the first seed by sampling uniformly at randomindex = zeros(1,k);[C(:,1), index(1)] = datasample(S,X,1,2);minDist = inf(n,1);% Select the rest of the seeds by a probabilistic modelfor ii = 2:k                    minDist = min(minDist,distfun(X,C(:,ii-1),distance));denominator = sum(minDist);if denominator==0 || isinf(denominator) || isnan(denominator)C(:,ii:k) = datasample(S,X,k-ii+1,2,'Replace',false);break;endsampleProbability = minDist/denominator;[C(:,ii), index(ii)] = datasample(S,X,1,2,'Replace',false,...'Weights',sampleProbability);        endendif ~isfloat(C)      % X may be logicalC = double(C);end% 计算点到聚类中心的距离和归属到各个类别D = distfun(X, C, distance, 0, rep, reps);%计算点到个中心的距离[d, idx] = min(D, [], 2);%根据最短距离归属到各个类m = accumarray(idx,1,[k,1])';%计算各个类中样本的个数try % catch空类错误并转移到下一个重复次%开始第一阶段:批分配converged = batchUpdate();% 开始第二阶段：单个分配if onlineconverged = onlineUpdate();endif display == 2 % 'final'fprintf('%s\n',getString(message('stats:kmeans:IterationsSumOfDistances',rep,iter,sprintf('%g',totsumD) )));endif ~convergedif reps==1warning(message('stats:kmeans:FailedToConverge', maxit));elsewarning(message('stats:kmeans:FailedToConvergeRep', maxit, rep));endend% 计算类内距离和nonempties = find(m>0);%判断没有空类，生成非空类的线性目录D(:,nonempties) = distfun(X, C(:,nonempties), distance, iter, rep, reps);d = D((idx-1)*n + (1:n)');sumD = accumarray(idx,d,[k,1]);% 计算类内距离和totsumD = sum(sumD(nonempties));% 计算所有类内距离和的总和% 保存目前最好的解cellout = {totsumD,rep,sumD,iter,idx,C,D}';% 如果在重复运行中发生空类现象,进行捕获并警告，然后继续下一次重复运行，%  只有在所有的重复运行失败才会ERROR，再次引发另一种ERROR。catch MEif reps == 1 || (~isequal(ME.identifier,'stats:kmeans:EmptyCluster')  && ...~isequal(ME.identifier,'stats:kmeans:EmptyClusterRep'))rethrow(ME);elseemptyErrCnt = emptyErrCnt + 1;warning(message('stats:kmeans:EmptyClusterInBatchUpdate', rep, iter));if emptyErrCnt == repserror(message('stats:kmeans:EmptyClusterAllReps'));endendend % catch%------------------------------------------------------------------function converged = batchUpdate()% 遍历每个点，更新每个类moved = 1:n;changed = 1:k;previdx = zeros(n,1);prevtotsumD = Inf;%% 开始第一阶段%iter = 0;converged = false;while trueiter = iter + 1;% 更新新的聚类中心和数目以及每个样本到新聚类中心的距离 [C(:,changed), m(changed)] = gcentroids(X, idx, changed, distance);D(:,changed) = distfun(X, C(:,changed), distance, iter, rep, reps);%处理空类empties = changed(m(changed) == 0);if ~isempty(empties)if strcmp(emptyact,'error')if reps==1error(message('stats:kmeans:EmptyCluster', iter));elseerror(message('stats:kmeans:EmptyClusterRep', iter, rep));endendswitch emptyactcase 'drop'if reps==1warning(message('stats:kmeans:EmptyCluster', iter));elsewarning(message('stats:kmeans:EmptyClusterRep', iter, rep));end% Remove the empty cluster from any further processingD(:,empties) = NaN;changed = changed(m(changed) > 0);case 'singleton'for i = emptiesd = D((idx-1)*n + (1:n)'); % use newly updated distances% 选取一个距离当前类最远的样本作为一个新的类[~, lonely] = max(d);from = idx(lonely); % taking from this clusterif m(from) < 2% In the very unusual event that the cluster had only% one member, pick any other non-singleton point.from = find(m>1,1,'first');lonely = find(idx==from,1,'first');endC(:,i) = X(:,lonely);m(i) = 1;idx(lonely) = i;D(:,i) = distfun(X, C(:,i), distance, iter, rep, reps);% Update clusters from which points are taken[C(:,from), m(from)] = gcentroids(X, idx, from, distance);D(:,from) = distfun(X, C(:,from), distance, iter, rep, reps);changed = unique([changed from]);endendend% 在当前配置下计算总距离totsumD = sum(D((idx-1)*n + (1:n)'));% 循环测试: 如果目标为减少,返回出去% 最后一步，之后进行单个更新阶段if prevtotsumD <= totsumDidx = previdx;[C(:,changed), m(changed)] = gcentroids(X, idx, changed, distance);iter = iter - 1;break;endif display > 2 % 'iter'if usePoolfprintf(dispfmt,labindx,iter,1,length(moved),totsumD);elsefprintf(dispfmt,iter,1,length(moved),totsumD);endendif iter >= maxitbreak;end%对每个点根据就近原则归属到各自的类 previdx = idx;prevtotsumD = totsumD;[d, nidx] = min(D, [], 2);% 决定哪个样本点移动moved = find(nidx ~= previdx);if ~isempty(moved)% Resolve ties in favor of not movingmoved = moved(D((previdx(moved)-1)*n + moved) > d(moved));endif isempty(moved)converged = true;break;endidx(moved) = nidx(moved);% 寻找得到或者失去样本点的类changed = unique([idx(moved); previdx(moved)])';end % phase oneend % nested function%------------------------------------------------------------------function converged = onlineUpdate()%% 第二阶段开始: 单个分配%changed = find(m > 0);lastmoved = 0;nummoved = 0;iter1 = iter;converged = false;Del = NaN(n,k); % 重新分配的准则while iter < maxit%计算每个样本点到各个类的距离以及因每个类中添加或者移除样本点引起的误差和的变化%没有发生变化的类并不用更新。仅含单个样本点的类是总距离计算中的特殊情况。%移除它们仅有的样本点并不是最好的选择，根据分配准则最好保证它们能够得到保留， %令人高兴地是，对于这种情况我们自动的使用Del(i,idx(i)) == 0。 switch distancecase 'sqeuclidean'for i = changedmbrs = (idx == i)';sgn = 1 - 2*mbrs; % -1 for members, 1 for nonmembersif m(i) == 1sgn(mbrs) = 0; % prevent divide-by-zero for singleton mbrsendDel(:,i) = (m(i) ./ (m(i) + sgn)) .* sum((bsxfun(@minus, X, C(:,i))).^2, 1);endcase {'cosine','correlation'}% The points are normalized, centroids are not, so normalize themnormC = sqrt(sum(C.^2, 1));if any(normC < eps(class(normC))) % small relative to unit-length data pointsif reps==1error(message('stats:kmeans:ZeroCentroid', iter));elseerror(message('stats:kmeans:ZeroCentroidRep', iter, rep));endend% This can be done without a loop, but the loop saves memory allocationsfor i = changedXCi =  C(:,i)'*X;mbrs = (idx == i)';sgn = 1 - 2*mbrs; % -1 for members, 1 for nonmembersDel(:,i) = 1 + sgn .*...(m(i).*normC(i) - sqrt((m(i).*normC(i)).^2 + 2.*sgn.*m(i).*XCi + 1));endend% 对于任意一个样本点，确定可能是最好的移动方式。然后选择其中的一个进行移动previdx = idx;prevtotsumD = totsumD;[minDel, nidx] = min(Del, [], 2);moved = find(previdx ~= nidx);moved(m(previdx(moved))==1)=[];if ~isempty(moved)% Resolve ties in favor of not movingmoved = moved(Del((previdx(moved)-1)*n + moved) > minDel(moved));endif isempty(moved)% Count an iteration if phase 2 did nothing at all, or if we're% in the middle of a pass through all the pointsif (iter == iter1) || nummoved > 0iter = iter + 1;if display > 2 % 'iter'if usePoolfprintf(dispfmt,labindx,iter,2,length(moved),totsumD);elsefprintf(dispfmt,iter,2,length(moved),totsumD);endendendconverged = true;break;end% Pick the next move in cyclic order%循环地选择下一次的移动moved = mod(min(mod(moved - lastmoved - 1, n) + lastmoved), n) + 1;% 遍历完所有的样本点，则完成一次迭代if moved <= lastmovediter = iter + 1;if display > 2 % 'iter'if usePoolfprintf(dispfmt,labindx,iter,2,length(moved),totsumD);elsefprintf(dispfmt,iter,2,length(moved),totsumD);endendif iter >= maxit, break; endnummoved = 0;endnummoved = nummoved + 1;lastmoved = moved;oidx = idx(moved);nidx = nidx(moved);totsumD = totsumD + Del(moved,nidx) - Del(moved,oidx);%更新类的索引向量、新旧类别各自的样本数目和中心idx(moved) = nidx;m(nidx) = m(nidx) + 1;m(oidx) = m(oidx) - 1;switch distancecase {'sqeuclidean','cosine','correlation'}C(:,nidx) = C(:,nidx) + (X(:,moved) - C(:,nidx)) / m(nidx);C(:,oidx) = C(:,oidx) - (X(:,moved) - C(:,oidx)) / m(oidx);endchanged = sort([oidx nidx]);end % phase twoend % nested functionendend % main function%------------------------------------------------------------------function D = distfun(X, C, dist, iter,rep, reps)
%DISTFUN计算样本点到中心的距离switch distcase 'sqeuclidean'D = pdist2mex(X,C,'sqe',[],[],[]);  case {'cosine','correlation'}% 样本点已被标准化而中心点没有,因此对它们进行标准化 normC = sqrt(sum(C.^2, 1));if any(normC < eps(class(normC))) % small relative to unit-length data points(?)if reps==1error(message('stats:kmeans:ZeroCentroid', iter));elseerror(message('stats:kmeans:ZeroCentroidRep', iter, rep));endendC = bsxfun(@rdivide,C,normC);D = pdist2mex(X,C,'cos',[],[],[]);
end
end % function%------------------------------------------------------------------
function [centroids, counts] = gcentroids(X, index, clusts, dist)
%GCENTROIDS Centroids and counts stratified by group.
%计算各类的样本数目和中心点
p = size(X,1);
num = length(clusts);centroids = NaN(p,num,'like',X);
counts = zeros(1,num,'like',X);for i = 1:nummembers = (index == clusts(i));if any(members)counts(i) = sum(members);switch distcase {'sqeuclidean','cosine','correlation'}centroids(:,i) = sum(X(:,members),2) / counts(i);endend
end
end % function

Python 中的Kmeans

from sklearn.cluster import KMeans
import numpy as np
X = np.array([[1, 2], [1, 4], [1, 0], [4, 2], [4, 4], [4,0]])
kmeans=KMeans(n_clusters=2,random_state=0).fit(X)

Kmeans聚类算法及其matlab源码相关推荐

C语言实现聚类K-means cluster算法(附完整源码)
聚类K-means cluster算法实现聚类K-means cluster算法的完整源码(定义,实现,main函数测试) 实现聚类K-means cluster算法的完整源码(定义,实现,main ...
matlab的数值计算方法,数值计算方法中的一些常用算法的Matlab源码
数值计算方法中的一些常用算法的Matlab源码,这些程序都是原创,传上来仅供大家参考,不足之处请大家指正,切勿做其它用途-- 说明:这些程序都是脚本函数,不可直接运行,需要创建函数m文件,保存时文件名 ...
多目标人工秃鹫优化算法（MATLAB源码分享，智能优化算法）提出了一种多目标版本的人工秃鹫优化算法（AVOA）
多目标人工秃鹫优化算法(MATLAB源码分享,智能优化算法) 提出了一种多目标版本的人工秃鹫优化算法(AVOA),用于多目标优化问题. AVOA的灵感来源于非洲秃鹫的生活方式. 档案.网格和领导者选择 ...
FA（萤火虫算法）MATLAB源码详细中文注解
以优化SVM算法的参数c和g为例,对FA(萤火虫算法)MATLAB源码进行了逐行中文注解. 完整程序和示例文件地址:http://download.csdn.net/detail/u013337691 ...
MATLAB实战系列（三十八）-基于K-means聚类算法的MATLAB图像分割
前言以下是我为大家准备的几个精品专栏,喜欢的小伙伴可自行订阅,你的支持就是我不断更新的动力哟! MATLAB-30天带你从入门到精通 MATLAB深入理解高级教程(附源码) tableau可视化数据 ...
k-means聚类算法及matlab实现(简单实现)
k-means简介 k-means算法也称k均值算法,是一种常用的聚类算法.聚类算法是研究最多.应用最广的一种无监督学习算法. 聚类试图将数据集中的样本划分为若干个通常是不相交的子集,每个子集 ...
K-Means聚类算法（matlab）
定义 k均值聚类算法(k-means clustering algorithm)是一种迭代求解的聚类分析算法,其步骤是,预将数据分为K组,则随机选取K个对象作为初始的聚类中心,然后计算每个对象与各个种 ...
python：实现DBSCAN聚类算法(附完整源码)
python:实现DBSCAN聚类算法 print(__doc__)# 引入相关包import numpy as npfrom sklearn.cluster import DBSCANfrom sk ...
【滤波跟踪】Singer-Kalman模型下的机动目标跟踪算法含Matlab源码
1 简介现实跟踪场景中,运动方式多样性.随机性运动和运动规则不确定性是机动目标的典型运动特征.机动目标跟踪的难点主要有: 建模一个准确的.通用的先验数学模型来表示机动目标的运动特性; 设计基于完善的 ...

Kmeans聚类算法及其matlab源码

K-means算法

matlab源码

Kmeans聚类算法及其matlab源码相关推荐

最新文章

热门文章