本期是聚类方法的第二讲——系统聚类法,第一讲中实现了简化版的kmeans聚类算法:
(注:原kmeans算法第37行存在bug,经群友“没有昵称”指出,已改正!)
matlab自带的系统聚类函数linkage功能比较复杂,定义了各种样本距离和类间距离,对于初学者而言不容易掌握方法的精髓。今天实现的简化版的系统聚类仅实现了欧几里得距离和汉明距离两种点距离,以及最小距离作为类与类之间距离,更容易理解。
系统聚类法(自下而上)先将聚类的样本或变量各自看成一群,然后确定类与类间的相似统计量,并选择最接近的两类或若干个类合并成一个新类,计算新类与其他各类间的相似性统计量,再选择最接近的两群或若干群合并成一个新类,直到所有的样本或变量都合并成一类为止。以n个样本的聚类分析为例,系统聚类法的步骤如下:
(1)定义k维空间中一种距离,k为样本维度;
(2)计算n个样本两两之间的距离;
(3)将每个样本归为一类,根据计算出的样本间的距离合并距离最近的两类为一个新类;
(4)再计算新类与其他各类的距离,同样再根据计算出的距离合并距离最近的两类为一个新类;
(5)循环以上过程直至类别个数为1。
hierarchicalCluster函数是我们实现的系统聚类函数,该函数调用caldis函数计算样本之间的距离,详细注释版代码如下:
function [clusterpath,clusterind]=hierarchicalCluster(data,distype)% 实现自下而上的系统聚类,样本和样本之间的距离采用欧几里得距离或者汉明距离,类与类之间的距离用两类点之间最小距离代替。% INPUTS:% data是聚类数据,每一行代表一个样本,每一列代表一个维度。% dis代表使用的度量距离的方法,这里实现了欧几里得距离和汉明距离%% OUTPUTS% clusterpath:聚类的途径,例如在十个样本的聚类例子当中,% clusterpath =% % [4 2 11% 11 1 12% 6 3 13% 12 7 14% 13 5 15% 15 8 16% 16 14 17% 17 9 18% 18 10 19]% 表明样本4和2首先聚在一起,形成类11,然后类11和点1聚在一起,形成类12% clusterind:聚类过程中每个中间节点包含的节点,例如上例中clusterind{11}包含4和2两个点。% 公众号【数学建模公会】,HCLO4原创,20190923% 计算所有样本两两之间的距离dis=caldis(data,distype);numsample=size(data,1); % 样本总数量clusterind=cell(1,2*numsample); % 存储所有中间聚类集合包含的元素for i=1:numsample % 最初始的叶子节点集合下标设置为相应的在data矩阵中的序号即可。 clusterind{i}=i;enddis_cluster=dis; % dis_cluster表示现有的cluster之间距离矩阵。tempclusterind=1:numsample;% 当前的dis_cluster距离矩阵对应的cluster编号clusterpath=zeros(2*numsample,3); % 记录聚类历史p=1;newind=numsample+1;while length(dis_cluster)>1 % 聚类结束的标志:所有点聚为一类 [ind1,ind2]=find(dis_cluster==min(dis_cluster(:)),1); % 寻找最小的距离的两个类 clusterpath(p,:)=[tempclusterind(ind1),tempclusterind(ind2),newind]; % 更新聚类路径 clusterind{newind}=[clusterind{tempclusterind(ind1)},clusterind{tempclusterind(ind2)}]; %更新类包含的点 tempclusterind=[setdiff(tempclusterind,tempclusterind([ind1,ind2])),newind]; numcluster=length(tempclusterind); dis_cluster=zeros(numcluster); for i=1:numcluster-1 for j=i+1:numcluster subdis=dis(clusterind{tempclusterind(i)},clusterind{tempclusterind(j)}); % 两个类中包含的点的相互距离矩阵 dis_cluster(i,j)=min(subdis(:)); % cluster之间点的距离的最小值代表两个cluster的距离 dis_cluster(j,i)=dis_cluster(i,j); end end dis_cluster(logical(eye(size(dis_cluster))))=inf; % 将对角线设置为无穷大 newind=newind+1; p=p+1; endclusterpath=clusterpath(1:p-1,:);clusterind=clusterind(1:newind-1);endfunction dis=caldis(data,distype)%计算每个样本到当前聚类中心的距离,distype='euclidean' 或者'hamming'numsample=size(data,1);dis=zeros(numsample);switch distype case 'euclidean' for i=1:numsample-1 for j=i+1:numsample dis(i,j)=sqrt(sum((data(i,:)-data(j,:)).^2)); dis(j,i)=dis(i,j); end end case 'hamming' for i=1:numsample-1 for j=i+1:numsample dis(i,j)=mean(data(i,:)~=data(j,:)); dis(j,i)=dis(i,j); end endenddis(logical(eye(size(dis))))=inf; % 给对角线元素赋值无穷,防止在聚类时将点和自身聚为一类end % caldis
以下代码是测试hierarchicalCluster的代码。
data=rand(10,2); % 随意生成一组待聚类样本,包含10个二维数据tic,distype='euclidean';[clusterpath,clusterind]=hierarchicalCluster(data,distype);tocdisp(clusterpath)