Matlab学习系列23模糊聚类分析原理及实现

23.模糊聚类分析原理及实现聚类分析，就是用数学方法研究和处理所给定对象，按照事物间的相似性进行区分和分类的过程。传统的聚类分析是一种硬划分，它把每个待识别的对象严格地划分到某个类中，具有非此即彼的性质，这种分类的类别界限是分明的。随着模糊理论的建立，人们开始用模糊的方法来处理聚类问题，称为模糊聚类分析。由于模糊聚类得到了样本数与各个类别的不确定性程度，表达了样本类属的中介性，即建立起了样本对于类别的不确定性的描述，能更客观地反映现实世界。本篇先介绍传统的两种（适合数据量较小情形，及理解模糊聚类原理）：基于择近原则、模糊等价关系的模糊聚类方法。（一）预备知识一、模糊等价矩阵定义1设R=(rij)n×n为模糊矩阵，I为n阶单位矩阵，若R满足i)自反性：I≤R（等价于rii=1）；ii)对称性：RT=R;则称R为模糊相似矩阵，若再满足iii)传递性：R2≤R（等价于1()nikkjijkrrr）则称R为模糊等价矩阵。定理1设R为n阶模糊相似矩阵，则存在一个最小的自然数k（kn）,使得Rk为模糊等价矩阵，且对一切大于k的自然数l，恒有Rl=Rk.Rk称为R的传递闭包矩阵，记为t(R).二、模糊矩阵的λ-截矩阵定义2设A=(aij)n×m为模糊矩阵，对任意的λ∈[0,1],作矩阵()ijnmAa其中，()1,0,ijijijaaa称为模糊矩阵A的λ-截矩阵。显然，Aλ为布尔矩阵，且其等价性与与A一致。意义：将模糊等价矩阵转化为等价的布尔矩阵，可以得到有限论域上的普通等价关系，而等价关系是可以分类的。因此，当λ在[0,1]上变动时，由Aλ得到不同的分类。若λ1＜λ2,则Aλ1≥Aλ2,从而由Aλ2确定的分类是由Aλ1确定的分类的加细。当λ从1递减变化到0时，Aλ的分类由细变粗，逐渐归并，形成一个分级聚类树。例1设U={u1,u2,u3,u4,u5},对给定的U上的模糊等价关系让λ从1到0变化，观察分类过程。(1)当λ=1时，11000001000001000001000001R分类结果为5类：（每行代表一类，1代表对应元素在该类）{u1},{u2},{u3},{u4},{u5}(2)当λ=0.8时，0.81010001000101000001000001R分类结果为4类：{u1,u3},{u2},{u4},{u5}(3)当λ=0.6时，0.61010001000101000001100011R分类结果为3类：{u1,u3},{u2},{u4,u5}(4)当λ=0.5时，0.51011101000101111011110111R分类结果为2类：{u1,u3,u4,u5},{u2}(4)当λ=0.4（R中的最小值）时，0.41111111111111111111111111R分类结果为1类：{u1,u2,u3,u4,u5}整个动态分类过程如下：（二）基于择近原则的模糊聚类择近原则就是利用贴近度来实现分类操作，贴近度用来衡量两个模糊集A和B的接近程度，用N(A,B)表示。贴近度越大，表明二者越接近。设论域有限或者在一定区间，即U={u1,u2,…,un}或U=[a,b],常用的贴近度有以下三种：(1)海明贴近度11(,)1|()()|niiiNABAuBun1(,)1|()()|dbiiaNABAuBuuba(2)欧氏贴近度12211(,)1[()()]niiiNABAuBun1221(,)1[()()]dbiiaNABAuBuuba(3)格贴近度(,)()()ccNABABAB其中，1()()niiiABAuBu.Matlab实现：格贴近度的实现函数fuz_closing.mfunctiony=fuz_closing(A,B,type)%要求A与B列数相同的行向量[m,n]=size(A);switchtypecase1%海明贴近度y=1-sum(abs(A-B))/n;case2%欧氏贴近度y=1-(sum(A-B).^2)^(1/2)/sqrt(n);case3%格贴近度y1=max(min(ones(m,n)-A,ones(m,n)-B));%ones(m,n)-A等于A^cy2=max(min(A,B));y=min(y1,y2);end例2设某产品的质量等级分为5级，其中一级有5种评判因素u1,u2,u3,u4,u5.每一等级的模糊集为B1={0.50.50.60.40.3}B2={0.30.30.40.20.2}B3={0.20.20.30.10.1}B4={0.10.10.20.10}B5={0.10.10.10.10}假设某产品各评判因素的值为A={0.40.30.20.10.2},问该产品属于哪个等级？代码：A=[0.40.30.20.10.2];B=[0.50.50.60.40.3;0.30.30.40.20.2;0.20.20.30.10.1;0.10.10.20.10;0.10.10.10.10];fori=1:5haiming(i)=fuz_closing(A,B(i,:),1);oushi(i)=fuz_closing(A,B(i,:),2);ge(i)=fuz_closing(A,B(i,:),3);endhaimingoushige运行结果：haiming=0.78000.92000.90000.86000.8400oushi=0.50810.91060.86580.68700.6422ge=0.40000.30000.20000.20000.1000可见样本A与各等级的格贴近度分别为0.4,0.3,0.2,0.2,0.1,故可认为该产品属于B1等级。若按令两种贴近度判断，该产品属于B2等级。（三）基于模糊等价关系的模糊聚类一、算法步骤1.样本数据归一化设X={x1,x2,…,xn}为要分类的n个样本，每个样本有m个指标，即xi={xi1,xi2,…,xim},i=1,2,..,n得到原始数据矩阵X=(xij)n×m.由于不同指标的数据量纲不同，为了使数据能够比较，要先对X做归一化处理。2.建立模糊相似矩阵R先建立样本xi与xj相似程度rij,进而构造模糊相似矩阵R=(rij)n×n建立rij常用的方法有：(1)相似系数法①夹角余弦法：12211mikjkkijmmikjkkkxxrxx②相关系数法：12211||||()()mikijkjkijmmikijkjkkxxxxrxxxx(2)距离法一般取rij=1-c(d(xi,xj))α,其中c和α为适当选取的参数，使得0≤rij≤1.常用的距离有：①海明距离：1(,)||mijikjkkdxxxx②欧氏距离：21(,)()mijikjkkdxxxx③切比雪夫距离：1(,)max||ijikjkkmdxxxx(3)贴近度法①最大最小法：11()()mikjkkijmikjkkxxrxx②算术平均最小法：11()1()2mikjkkijmikjkkxxrxx③几何平均最小法：11()mikjkkijmikjkkxxrxx3.求出R的传递闭包t(R)即改造相似关系为等价关系：令2RRR,再令422RRR,…,直到满足2lllRRR与Rl相等，即为t(R),仍记为R.4.选取合适的λ,利用λ-截矩阵Rλ进行分类（参考例1）。二、Matlab实现求模糊相似矩阵R的函数：fuz_distance.mfunctionR=fuz_distance(x,type)%x为归一化的数据矩阵,type选择计算相似程度的方法%返回模糊相似矩阵R[n,m]=size(x);%距离法的选择参数c和a,需要根据具体情况修改以保证R(i,j)属于[0,1]c=0.1;a=1;fori=1:nforj=1:nswitchtypecase1%夹角余弦法R(i,j)=(x(i,:)*x(j,:)')/(norm(x(i,:),2)*norm(x(j,:),2));case2%相关系数法Dxi=abs(x(i,:)-mean(x(i,:)));Dxj=abs(x(j,:)-mean(x(j,:)));R(i,j)=(Dxi*Dxj')/(norm(Dxi,2)*norm(Dxj,2));case3%海明距离法d=sum(abs(x(i,:)-x(j,:)));R(i,j)=1-c*d^a;case4%欧氏距离法d=norm(x(i,:)-x(j,:),2);R(i,j)=1-c*d^a;case5%切比雪夫距离法d=max(abs(x(i,:)-x(j,:)));R(i,j)=1-c*d^a;case6最大最小(贴近度)法R(i,j)=sum(min([x(i,:);x(j,:)]))/sum(max([x(i,:);x(j,:)]));case7算术平均最小(贴近度)法R(i,j)=2*sum(min([x(i,:);x(j,:)]))/sum(x(i,:)+x(j,:));case8%几何平均最小(贴近度)法R(i,j)=sum(min([x(i,:);x(j,:)]))/sum(sqrt(x(i,:).*x(j,:)));endendend求R的传递闭包t(R)的函数：tran_R.mfunction[B,k]=tran_R(R)%R为模糊相似矩阵,循环构造满足传递性的t(R)%k为满足R^2k=R^k的最小的自然数kn=length(R);B=zeros(n,n);flag=0;k=1/2;whileflag==0B=fco(R,R);%做模糊合成运算k=2*k;ifB==Rflag=1;elseR=B;%循环计算R传递闭包endend上面的函数tran_R.m调用函数矩阵模糊合成算子函数：fco.mfunctionB=fco(Q,R)%实现模糊合成算子的计算,要求Q的列数等于R的行数[n,m]=size(Q);[m,l]=size(R);B=zeros(n,l);fori=1:nfork=1:lB(i,k)=max(min([Q(i,:);R(:,k)']));endend求t(R)的λ-截矩阵的函数：fuz_lamda.mfunctiony=fuz_lamda(X,m)%用λ-截矩阵将样本分成m类,m≤总样本数lamda=unique(X)';%根据R中的值取λ值%unique函数取矩阵不重复元素组成向量并从小到大排好序X(find(Xlamda(m)))=0;X(find(X=lamda(m)))=1;y=X;例3某地区设有11个雨量站，其分布如图所示：10年来各雨量站测得的年降雨量表如下：现因经费问题，希望撤销几个雨量站，问撤销哪些雨量站而不会太多地减少降雨信息？分析：对11个雨量站进行模糊聚类，同一类的只需保留一个即可。比如，已知该市决定撤销6个只保留5个雨量站，则模糊聚类为5类。代码：loaddata;%数据归一化[X,ps]=mapminmax(data',0,1);X=X';%选择计算相似程度的方法type=3;%c=0.1,a=1,此时也称绝对值减数法%求模糊相似矩阵R0R0=fuz_distance(X,type)%将模糊相似矩阵R0改造成模糊等价矩阵R[R,k]=tran_R(R0)%求将样本分成8类的λ-截矩阵R_lamda=fuz_lamda(R,8)运行结果及说明：归一化后的数据矩阵X:模糊相似矩阵R0:由R0改造成的模糊等价矩阵R:k=8说明R16=R8.将样本分为5类的λ-截矩阵R_lamda:可以判断5类分别是：{x1,x7}{x2,x4,x5,x6}{x3,x9}{x8,x11}{x10}注：对于这类C均值模糊聚类问题，也可以直接调用Matla