精华内容
下载资源
问答
  • kmeans聚类算法matlab代码
    千次阅读
    2019-09-06 13:26:27

    kmeans聚类算法是一种简单实用的聚类算法,matlab自带函数kmeans可直接对数据进行kmeans聚类。为了方便更好地掌握kmeans聚类算法,今天我们自己来实现一个弱化的版本mykmeans。

    mykmeans输入包含三项,分别为聚类所使用的数据data,data每一行代表一个样本,每一列代表一个特征;聚类中心数量numclass;第三项为所使用的距离的定义,默认情况下为欧式距离。

    function [cluster,clusterhistory]=mykmeans(data,numclass,varargin)
    % kmeans聚类。
    % 聚类过程动画显示
    % INPUTS:
    % data:每一行代表一个样本,每一列代表一个特征
    % numclass:聚类中心的数量
    % varargin{1}: 距离定义。支持euclidean(默认)、hamming
    % OUTPUT:
    % cluster: numsample行的列向量,代表每个样本的分类归属
    % clusterhistory{i}.cluster第i次迭代的聚类
    % clusterhistory{i}.core第i次迭代的聚类中心
    %
    % 公众号【数学建模公会】,HCLO4,20190902
    
    if nargin==2
        method='euclidean';
    else
        method=varargin{1};
    end
    % 初始化聚类中心坐标,在数据点中随机选择numclass个点作为初始聚类中心。
    numsample=size(data,1);
    temp=randperm(numsample);
    core=data(temp(1:numclass),:);
    dis=caldis(data,core,method); %存储每个样本到当前聚类中心的距离
    
    
    % 执行迭代过程
    maxIter=20; % 最大迭代次数
    numiter=1;
    clusterhistory=cell(1,maxIter);
    while 1
        
        newcore=zeros(size(core));
        % 迭代聚类中心
        [~,ind]=min(dis,[],2); %计算每个sample归属于哪个聚类中心,如果某个聚类中心没有一个点?
        for i=1:numclass
            newcore(i,:)=mean(dis(ind==i,:));
        end
        clusterhistory{numiter}.cluster=ind;
        clusterhistory{numiter}.core=core;
        
        if all(newcore(:)==core(:))||numiter>=maxIter  % 迭代终止条件,聚类中心不再改变
            cluster=ind;
            break
        end
        
        core=newcore;
        dis=caldis(data,core,method); 
        
        numiter=numiter+1;
        
    end
        
    
    clusterhistory=clusterhistory(1:numiter);
    
    
    end % mykmeans
    
    
    
    function dis=caldis(data,core,method)
    %计算每个样本到当前聚类中心的距离
    numsample=size(data,1);
    numclass=size(core,1);
    dis=zeros(numsample,numclass);
    switch method
        case 'euclidean'
            for i=1:numclass
                dis(:,i)=sqrt(sum((data-repmat(core(i,:),numsample,1)).^2,2));
            end
        case 'hamming'
            for i=1:numclass
                dis(:,i)=mean(data~=repmat(core(i,:),numsample,1),2);
            end
    end
    

    下面一段代码用于测试mykmeans和kmeans的运行结果。生成两组服从二维高斯分布的数据,作为两个数据类别,分别使用mykmeans和matlab自带的kmeans函数进行聚类分析。

    % 测试mykmeans函数
    % 公众号【数学建模公会】,HCLO4,20190902
    
    % 生成模拟数据,平面上的两组点,均服从二维高斯分布
    mu1=[2,8];
    sigma1=[2,2];
    mu2=[8,2];
    sigma2=[3,3];
    
    numsample1=100;
    numsample2=200;
    data1=zeros(numsample1,2);
    data1(:,1)=normrnd(mu1(1),sigma1(1),numsample1,1);
    data1(:,2)=normrnd(mu1(2),sigma1(2),numsample1,1);
    
    data2=zeros(numsample2,2);
    data2(:,1)=normrnd(mu2(1),sigma2(1),numsample2,1);
    data2(:,2)=normrnd(mu2(2),sigma2(2),numsample2,1);
    
    data=[data1;data2];
    
    tic,
    [cluster1,clusterhistory]=mykmeans(data,2);
    toc
    M=getFrame_Kmeans(data,clusterhistory); % 生成聚类动图
    
    tic,
    cluster2=kmeans(data,2);
    toc
    

    下面的函数实现mykmeans聚类过程的可视化,仅针对二维数据和三维数据类型。

    function M=getFrame_Kmeans(data,clusterhistory)
    % kmeans聚类过程可视化程序。只能对二维和三维数据可视化!
    %
    % INPUTS: 
    % data: 聚类用的数据,每行代表一个样本,每列代表一个特征
    % clusterhistory: mykmeans函数输出的聚类过程数据
    %
    % OUTPUT:
    % kmeans聚类过程的动画。
    
    dimension=size(data,2);
    if dimension>3 % 若三维以上,只能通过pca降维处理
        error('无法实现高于三维的数据的聚类可视化')
    end
    
    colorset=cell(1,1000);
    colorset(1:8)={'r','g','b','k','c','m','y','k'};
    for i=9:1000 %如果聚类中心数量大于8,就使用随机的颜色。
        colorset{i}=rand(1,3);
    end
    
    numcore=length(unique(clusterhistory{1}.cluster));
    numiter=length(clusterhistory);
    
    for k=1:numiter
        figure
        hold on
        switch dimension
            case 2
                for i=1:numcore
                    ind=clusterhistory{k}.cluster==i;
                    scatter(data(ind,1),data(ind,2),6,colorset{i})
                    core=clusterhistory{k}.core;
                    plot(core(i,1),core(i,2),[colorset{i},'.'],'MarkerSize',20)
                end
            case 3
                for i=1:numcore
                    ind=clusterhistory{k}.cluster==i;
                    scatter3(data(ind,1),data(ind,2),data(ind,3),6,colorset{i})
                    core=clusterhistory{k}.core;
                    plot3(core(i,1),core(i,2),core(i,3),colorset{i},'MarkerSize',6)
                end
        end
        
        M(k)=getframe(gcf);
        frame = getframe(gcf); 
        im = frame2im(frame);     %将影片动画转换为编址图像,因为图像必须是index索引图像
        imshow(im);
        [I,map] = rgb2ind(im,20); %将真彩色图像转化为索引图像
        if k==1
            imwrite(I,map,'kmeans.gif','gif','Loopcount',inf,'DelayTime',0.3);     %Loopcount只是在i==1的时候才有用
        else
            imwrite(I,map,'kmeans.gif','gif','WriteMode','append','DelayTime',1);%DelayTime:帧与帧之间的时间间隔
        end
        
        
        close(gcf);
        
    end
    
    更多相关内容
  • 该程序获取图像和所需的分区数,并找到不同类别的均值并提供分类图像(面具)。
  • k-means聚类算法matlab代码 数据挖掘实验 实验一:相似度、距离、最近邻分类器 1、实验目的 (1)理解相似度、距离的度量方式。 (2)理解最近邻分类器的工作原理。 2、实验内容 (1)、实现任意给定两个相同维度...
  • 该课题为基于kmeans聚类分割,输入一张彩色图像,可以选择需要分割成多少类,就会以不同颜色区分不同的块。
  • Kmeans聚类,kmeans聚类算法,matlab源码.rar
  • kmeans聚类分析,无监督学习实现Matlab代码
  • KMEANS聚类算法MATLAB代码 Algorithm 用Python,Matlab写的一些算法 \ (主目录) 文件名:算法名_功能 DeepLearning 来自吴恩达的深度学习课程 IntelligentAlgorithm 智能算法的代码 粒子群 模拟退火 鱼群算法 ...
  • matlab实现Kmeans聚类算法,非常详细。
  • K-means算法是很典型的基于距离的聚类算法,k-medoids聚类算法同样有效
  • K-means聚类算法利用matlab实现,可以查看每次迭代的效果
  • 对数据使用kmeans
  • 分析matlab代码K均值聚类 这是K-means算法在MATLAB和Python中的简单实现 K-means 聚类是一种矢量量化方法,最初来自信号处理,在数据挖掘中流行用于聚类分析。 k-means聚类旨在将n个观测值划分为k个簇,其中每个观测...
  • matlab均值聚类的基本代码代码参考自周志华《机器学习》9.4.1节而写的,代码使用matlab矩阵序列化操作,速度会快一点,本代码仅供参考,请尊重原创
  • kmeans聚类分析,无监督学习实现Matlab代码
  • k-means聚类算法matlab代码 数据挖掘-实验 . 目录 实验内容 实验说明文档 :link: 第三方库 数据处理模块: | 数据可视化模块: | 仓库文件内容说明 文件/目录 说明 实验一 代码 实验二 代码 运行结果 实验三 代码 ...
  • 经过简化的kmeans算法matlab程序,可满足基本需求。
  • 对数据使用kmeans
  • kmeans聚类算法,kmeans聚类算法优缺点,matlab源码.rar
  • kmeans聚类算法,kmeans聚类算法优缺点,matlab源码.zip
  • Kmeans聚类,kmeans聚类算法,matlab源码.zip
  • Kmeans聚类matlab代码

    2015-06-10 20:13:55
    Kmeans聚类算法matlab代码,可用于图像分割等数字图像处理领域。
  • 怎样用matlab实现多维K-means聚类算法小编觉得一个好的周末应该是这样的:睡到中午醒来,在床上躺着玩两个小时...怎样用matlab实现多维k-means聚类算法function [ labels ] = kmeans_clustering( data, k ) [num,~]...

    怎样用matlab实现多维K-means聚类算法小编觉得一个好的周末应该是这样的:睡到中午醒来,在床上躺着玩两个小时手机,起床随便吃点东西,下午去超市买一大堆零食,五六点的时候去约小伙伴们吃火锅烧烤,如果是一个人的话就去吃炸鸡,看电影。

    怎样用matlab实现多维k-means聚类算法

    function [ labels ] = kmeans_clustering( data, k ) [num,~]=size(data); ind = randperm(num); ind = ind(1:k); centers = data(ind,:); d=inf; labels = nan(num,1); while d>0 labels0 = labels; dist = pdist2(data, centers); [~,labels]小编不敢说自己一生都会喜欢你,但至少在能看见你的岁月里只想对你一个人好。

    1e9eada21267e1de1766ebec5772fe0e.png

    如何编写分享K-均值聚类算法的Matlab程序?喜欢,努力却无玩耍则心生疲惫,身累。喜欢,玩耍却不努力则无长进,心累。努力的玩耍而不喜爱则属应酬,身心俱累。

    在聚类分析中,K-均值聚类算法(k-means algorithm)是无监督分类中的一种基本方法,其也称为C-均值算法,其基本思想是:通过迭代的方法,逐次更新各聚类中心的值,直至得到最好的聚类结果。假设要把样本集分为c个类别。

    如何对点进行k均值聚类算法 matlab

    在聚类分析中,K-均值聚类算法(k-means algorithm)是无监督分类中的一种基本方法,其也称为C-均值算法,其基本思想是:通过迭代的方法,逐次更新各聚类中心的值,直至得到最好的聚类结果.\x0d假设要把样本集分为c个类别,算法如下:\x0d(1)适当选当一个人不再对一个人畅所欲言时,很多时候是心远了,当一个人学会默默承担时,很多时候是心凉了……

    急分享。。。matlab的K-均值聚类算法程序,采用下面采用下面的数据进行聚类分析: x1 x2 -0.5200 1.8539 2.5849 2.2481 0.9调用kmeans函数自甘堕落的人根本不值得别人可怜,越是没有人爱,就越要爱自己.

    matlab实现怎么kmeans算法总有一天你恍然大悟,父母是你花心思,花时间最少,却最爱你的人。

    利用kmeans算法对矩阵进行聚类,各个矩阵维数一样,使用1范数(两矩阵差小编开始喜欢最初的自己,那时候没有伤,不会哭泣。

    分享自适应k均值聚类算法 matlab

    function group=Kmeans(k,mid) %K均值聚类算法 person=load('sample.txt','height','weight');%从文本文件读入数据放入person结构体中 % person=person(1:10,:); num=size(person,1);%获得person结构体大小 for i=1:k%赋初始值,划分为k类。

    如何将k-means算法与de算法结合用matlab实现

    data=input('请输入样本数据矩阵:'); m=size(data,1); n=size(data,2); counter=0; k=input('请输入聚类数目:'); whilek>m disp('您输入的聚类数目过大,请输入正确的 k 值'); k=input('请输入聚类数目:'); end if k==1 disp('聚类数目不能为有一种默默的爱,叫放弃,轻轻地叹息,为了爱而离开。

    matlab 聚类算法silhouette小编不是那种几个飞吻就能打发的人,有本事你真亲过来呀。

    展开全文
  • K-means聚类算法 简介 聚类是一个将数据集中在某些方面相似的数据成员进行分类组织的过程,聚类就是一种发现这种内在结构的技术,聚类技术经常被称为无监督学习。 K均值聚类是最著名的划分聚类算法,由于简洁和效率...
  • 实验报告——Kmeans聚类方法.docx
  • Kmeans聚类算法及其matlab源码

    万次阅读 多人点赞 2016-10-24 14:57:34
    本文介绍了K-means聚类算法,并注释了部分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 hadNaNs
        warning(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 distance
        case 'cosine'
            Xnorm = sqrt(sum(X.^2, 2));%模长
            if any(min(Xnorm) <= eps(max(Xnorm)))
                error(message('stats:kmeans:ZeroDataForCos'));
            end
            X =  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'));
            end
            X =  bsxfun(@rdivide,X,Xnorm);
         case 'hamming'
           if  ~all( X(:) ==0 | X(:)==1)
                error(message('stats:kmeans:NonbinaryDataForHamm'));
          end
    end
    
    Xmins = [];
    Xmaxs = [];
    CC = [];
    if ischar(start)
        startNames = {'uniform','sample','cluster','plus','kmeans++'};
        j = find(strncmpi(start,startNames,length(start)));
        if length(j) > 1
            error(message('stats:kmeans:AmbiguousStart', start));
        elseif isempty(j)
            error(message('stats:kmeans:UnknownStart', start));
        elseif isempty(k)
            error(message('stats:kmeans:MissingK'));
        end
        start = startNames{j};
        if strcmp(start, 'uniform')
            if strcmp(distance, 'hamming')
                error(message('stats:kmeans:UniformStartForHamm'));
            end
            Xmins = 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'));
        end
        if 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
    else
        error(message('stats:kmeans:InvalidStart'));
    end
    
    emptyactNames = {'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);
    end
    
    options = 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 < k
        error(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' only
                warning(message('stats:kmeans:displayParallel2'));
                fprintf('    worker\t  iter\t phase\t     num\t         sum\n' );
            end
        else
            if useParallel
                warning(message('stats:kmeans:displayParallel'));
            end
            if display == 3 % 'iter' only
                fprintf('  iter\t phase\t     num\t         sum\n');
            end
        end
    end
    
    if 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;
    end
    
    emptyErrCnt = 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 > 3
        varargout{4} = ClusterBest{7}; %最优解的点到任意聚类中心的距离
    end
    
    if display > 1 % 'final' or 'iter'
        fprintf('%s\n',getString(message('stats:kmeans:FinalSumOfDistances',sprintf('%g',totsumDbest))));
    end
     
    if hadNaNs
        varargout{1} = statinsertnan(wasnan, varargout{1});% idxbest 
        if nargout > 3
            varargout{4} = statinsertnan(wasnan, varargout{4}); %Dbest
        end
    end
    
        function cellout = loopBody(rep,S)%循环体函数
            
            if isempty(S)
                S = RandStream.getGlobalStream;
            end
            
            if display > 1 % 'iter'
                if usePool
                    dispfmt = '%8d\t%6d\t%6d\t%8d\t%12g\n';
                    labindx = internal.stats.parallel.workerGetValue('workerID');
                else
                    dispfmt = '%6d\t%6d\t%8d\t%12g\n';
                end
            end
    
            %定义元胞数组
            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 start
                case '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));
                    end
                    if isa(X,'single')
                        C = single(C);
                    end
                case 'sample'
                    C = X(:,randsample(S,n,k));
                case 'cluster'
                    Xsubset = X(:,randsample(S,n,floor(.1*n)));
                    % Turn display off for the initialization
                    optIndex = find(strcmpi('options',varargin));
                    if isempty(optIndex)
                        opts = statset('Display','off');
                        varargin = [varargin,'options',opts];
                    else
                        varargin{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);
                    end
                case {'plus','kmeans++'}
                    % Select the first seed by sampling uniformly at random
                    index = 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 model
                   for 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;
                        end
                        sampleProbability = minDist/denominator;
                        [C(:,ii), index(ii)] = datasample(S,X,1,2,'Replace',false,...
                            'Weights',sampleProbability);        
                    end
            end
            if ~isfloat(C)      % X may be logical
                C = 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 online
                    converged = onlineUpdate();
                end
                
                
                if display == 2 % 'final'
                    fprintf('%s\n',getString(message('stats:kmeans:IterationsSumOfDistances',rep,iter,sprintf('%g',totsumD) )));
                end
                
                if ~converged
                    if reps==1
                        warning(message('stats:kmeans:FailedToConverge', maxit));
                    else
                        warning(message('stats:kmeans:FailedToConvergeRep', maxit, rep));
                    end
                end
                
                % 计算类内距离和
                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 ME
                if reps == 1 || (~isequal(ME.identifier,'stats:kmeans:EmptyCluster')  && ...
                             ~isequal(ME.identifier,'stats:kmeans:EmptyClusterRep'))
                    rethrow(ME);
                else
                    emptyErrCnt = emptyErrCnt + 1;
                    warning(message('stats:kmeans:EmptyClusterInBatchUpdate', rep, iter));
                    if emptyErrCnt == reps
                        error(message('stats:kmeans:EmptyClusterAllReps'));
                    end
                end
            end % catch
            
            %------------------------------------------------------------------
            
            function converged = batchUpdate()
                
                % 遍历每个点,更新每个类
                moved = 1:n;
                changed = 1:k;
                previdx = zeros(n,1);
                prevtotsumD = Inf;
                
                %
                % 开始第一阶段
                %
                
                iter = 0;
                converged = false;
                while true
                    iter = 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==1
                                error(message('stats:kmeans:EmptyCluster', iter));
                            else
                                error(message('stats:kmeans:EmptyClusterRep', iter, rep));
                            end
                        end
                        switch emptyact
                            case 'drop'
                                if reps==1
                                    warning(message('stats:kmeans:EmptyCluster', iter));
                                else
                                    warning(message('stats:kmeans:EmptyClusterRep', iter, rep));
                                end
                                % Remove the empty cluster from any further processing
                                D(:,empties) = NaN;
                                changed = changed(m(changed) > 0);
                            case 'singleton'
                                for i = empties
                                    d = D((idx-1)*n + (1:n)'); % use newly updated distances
                                    
                                    % 选取一个距离当前类最远的样本作为一个新的类
                                    [~, lonely] = max(d);
                                    from = idx(lonely); % taking from this cluster
                                    if 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');
                                    end
                                    C(:,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]);
                                end
                        end
                    end
                    
                    % 在当前配置下计算总距离
                    totsumD = sum(D((idx-1)*n + (1:n)'));
                    % 循环测试: 如果目标为减少,返回出去
                    % 最后一步,之后进行单个更新阶段
                    if prevtotsumD <= totsumD
                        idx = previdx;
                        [C(:,changed), m(changed)] = gcentroids(X, idx, changed, distance);
                        iter = iter - 1;
                        break;
                    end
                    if display > 2 % 'iter'
                        if usePool
                            fprintf(dispfmt,labindx,iter,1,length(moved),totsumD);
                        else
                            fprintf(dispfmt,iter,1,length(moved),totsumD);
                        end
                    end
                    if iter >= maxit
                        break;
                    end
                    
                    %对每个点根据就近原则归属到各自的类 
                    previdx = idx;
                    prevtotsumD = totsumD;
                    [d, nidx] = min(D, [], 2);
                    
                    % 决定哪个样本点移动
                    moved = find(nidx ~= previdx);
                    if ~isempty(moved)
                        % Resolve ties in favor of not moving
                        moved = moved(D((previdx(moved)-1)*n + moved) > d(moved));
                    end
                    if isempty(moved)
                        converged = true;
                        break;
                    end
                    idx(moved) = nidx(moved);
                    
                    % 寻找得到或者失去样本点的类
                    changed = unique([idx(moved); previdx(moved)])';
                    
                end % phase one
                
            end % 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 distance
                        case 'sqeuclidean'
                            for i = changed
                                mbrs = (idx == i)';
                                sgn = 1 - 2*mbrs; % -1 for members, 1 for nonmembers
                                if m(i) == 1
                                    sgn(mbrs) = 0; % prevent divide-by-zero for singleton mbrs
                                end
                              Del(:,i) = (m(i) ./ (m(i) + sgn)) .* sum((bsxfun(@minus, X, C(:,i))).^2, 1);
                            end
                          case {'cosine','correlation'}
                            % The points are normalized, centroids are not, so normalize them
                            normC = sqrt(sum(C.^2, 1));
                            if any(normC < eps(class(normC))) % small relative to unit-length data points
                                if reps==1
                                    error(message('stats:kmeans:ZeroCentroid', iter));
                                else
                                    error(message('stats:kmeans:ZeroCentroidRep', iter, rep));
                                end
                                
                            end
                            % This can be done without a loop, but the loop saves memory allocations
                            for i = changed
                                XCi =  C(:,i)'*X;
                                mbrs = (idx == i)';
                                sgn = 1 - 2*mbrs; % -1 for members, 1 for nonmembers
                                Del(:,i) = 1 + sgn .*...
                                    (m(i).*normC(i) - sqrt((m(i).*normC(i)).^2 + 2.*sgn.*m(i).*XCi + 1));
                            end
                    end
                    
                    % 对于任意一个样本点,确定可能是最好的移动方式。然后选择其中的一个进行移动
                    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 moving
                        moved = moved(Del((previdx(moved)-1)*n + moved) > minDel(moved));
                    end
                    if 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 points
                        if (iter == iter1) || nummoved > 0
                            iter = iter + 1;
                            if display > 2 % 'iter'
                                if usePool
                                    fprintf(dispfmt,labindx,iter,2,length(moved),totsumD);
                                else
                                    fprintf(dispfmt,iter,2,length(moved),totsumD);
                                end
                            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 moved <= lastmoved
                        iter = iter + 1;
                        if display > 2 % 'iter'
                            if usePool
                                fprintf(dispfmt,labindx,iter,2,length(moved),totsumD);
                            else
                                fprintf(dispfmt,iter,2,length(moved),totsumD);
                            end
                        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);
                    
                    %更新类的索引向量、新旧类别各自的样本数目和中心
                    idx(moved) = nidx;
                    m(nidx) = m(nidx) + 1;
                    m(oidx) = m(oidx) - 1;
                    switch distance
                        case {'sqeuclidean','cosine','correlation'}
                            C(:,nidx) = C(:,nidx) + (X(:,moved) - C(:,nidx)) / m(nidx);
                            C(:,oidx) = C(:,oidx) - (X(:,moved) - C(:,oidx)) / m(oidx);
                    end
                    changed = sort([oidx nidx]);
                end % phase two
                
            end % nested function
            
        end
    
    end % main function
    
    %------------------------------------------------------------------
    
    function D = distfun(X, C, dist, iter,rep, reps)
    %DISTFUN计算样本点到中心的距离
    
    switch dist
        case '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==1
                    error(message('stats:kmeans:ZeroCentroid', iter));
                else
                    error(message('stats:kmeans:ZeroCentroidRep', iter, rep));
                end
                
            end
            C = 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:num
        members = (index == clusts(i));
        if any(members)
           counts(i) = sum(members);
           switch dist
               case {'sqeuclidean','cosine','correlation'}
                  centroids(:,i) = sum(X(:,members),2) / counts(i);
          end
        end
    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和dbscan的聚类算法matlab实现
  • 这是 kmeans 聚类算法的超级快速实现。 代码完全矢量化并且非常简洁。 它比 Matlab 内置的 kmeans 函数快得多。 还包括 kmeans++ 播种算法 (kseeds.m) 以实现良好的初始化。 所以,这个包不仅是为了炫酷,它确实是...
  • 该课题为基于kmeans聚类分割,输入一张彩色图像,可以选择需要分割成多少类,就会以不同颜色区分不同的块,带有GUI界面,操作丰富。

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 2,656
精华内容 1,062
关键字:

kmeans聚类算法matlab代码

matlab 订阅