模糊C均值聚类代码参考(Matlab) - leivo

时间：2024-03-08 20:43:18
  1. FCM算法的两种迭代形式的MATLAB代码写于下,也许有的同学会用得着:
  2. m文件1/7:
  3. function [U,P,Dist,Cluster_Res,Obj_Fcn,iter]=fuzzycm(Data,C,plotflag,M,epsm)
  4. % 模糊 C 均值聚类 FCM: 从随机初始化划分矩阵开始迭代
  5. % [U,P,Dist,Cluster_Res,Obj_Fcn,iter] = fuzzycm(Data,C,plotflag,M,epsm)
  6. % 输入:
  7. %     Data: N×S 型矩阵,聚类的原始数据,即一组有限的观测样本集,
  8. %           Data 的每一行为一个观测样本的特征矢量,S 为特征矢量
  9. %           的维数,N 为样本点的个数
  10. %     C:    聚类数,1<C<N
  11. %     plotflag: 聚类结果 2D/3D 绘图标记,0 表示不绘图,为缺省值        
  12. %     M:    加权指数,缺省值为 2
  13. %     epsm: FCM 算法的迭代停止阈值,缺省值为 1.0e-6
  14. % 输出:
  15. %     U:    C×N 型矩阵,FCM 的划分矩阵
  16. %     P:    C×S 型矩阵,FCM 的聚类中心,每一行对应一个聚类原型
  17. %     Dist: C×N 型矩阵,FCM 各聚类中心到各样本点的距离,聚类中
  18. %           心 i 到样本点 j 的距离为 Dist(i,j)
  19. %     Cluster_Res: 聚类结果,共 C 行,每一行对应一类
  20. %     Obj_Fcn: 目标函数值
  21. %     iter: FCM 算法迭代次数
  22. % See also: fuzzydist maxrowf fcmplot
  23. if nargin<5
  24.     epsm=1.0e-6; 
  25. end
  26. if nargin<4
  27.     M=2;
  28. end
  29. if nargin<3
  30.     plotflag=0;
  31. end
  32. [N,S]=size(Data);m=2/(M-1);iter=0;
  33. Dist(C,N)=0; U(C,N)=0; P(C,S)=0;
  34. % 随机初始化划分矩阵
  35. U0 = rand(C,N); 
  36. U0=U0./(ones(C,1)*sum(U0));
  37. % FCM 的迭代算法
  38. while true 
  39.     % 迭代计数器
  40.     iter=iter+1; 
  41.     % 计算或更新聚类中心 P
  42.     Um=U0.^M;
  43.     P=Um*Data./(ones(S,1)*sum(Um\'))\';   
  44.     % 更新划分矩阵 U
  45.     for i=1:C
  46.         for j=1:N
  47.             Dist(i,j)=fuzzydist(P(i,:),Data(j,:));
  48.         end
  49.     end         
  50.     U=1./(Dist.^m.*(ones(C,1)*sum(Dist.^(-m))));          
  51.     % 目标函数值: 类内加权平方误差和
  52.     if nargout>4 | plotflag
  53.         Obj_Fcn(iter)=sum(sum(Um.*Dist.^2));
  54.     end
  55.     % FCM 算法迭代停止条件
  56.     if norm(U-U0,Inf)<epsm
  57.         break
  58.     end
  59.     U0=U;   
  60. end
  61. % 聚类结果
  62. if nargout > 3
  63.     res = maxrowf(U);
  64.     for c = 1:C
  65.         v = find(res==c);
  66.         Cluster_Res(c,1:length(v))=v;
  67.     end
  68. end
  69. % 绘图
  70. if plotflag
  71.     fcmplot(Data,U,P,Obj_Fcn);
  72. end
  73. m文件2/7:
  74. function [U,P,Dist,Cluster_Res,Obj_Fcn,iter]=fuzzycm2(Data,P0,plotflag,M,epsm)
  75. % 模糊 C 均值聚类 FCM: 从指定初始聚类中心开始迭代
  76. % [U,P,Dist,Cluster_Res,Obj_Fcn,iter] = fuzzycm2(Data,P0,plotflag,M,epsm)
  77. % 输入: Data,plotflag,M,epsm: 见 fuzzycm.m
  78. %       P0: 初始聚类中心
  79. % 输出: U,P,Dist,Cluster_Res,Obj_Fcn,iter: 见 fuzzycm.m    
  80. % See also: fuzzycm
  81. if nargin<5
  82.     epsm=1.0e-6; 
  83. end
  84. if nargin<4
  85.     M=2;
  86. end
  87. if nargin<3
  88.     plotflag=0;
  89. end
  90. [N,S] = size(Data); m = 2/(M-1); iter = 0;
  91. C=size(P0,1);Dist(C,N)=0;U(C,N)=0;P(C,S)=0;
  92. % FCM 的迭代算法
  93. while true 
  94.     % 迭代计数器
  95.     iter=iter+1; 
  96.     % 计算或更新划分矩阵 U
  97.     for i=1:C
  98.         for j=1:N
  99.             Dist(i,j)=fuzzydist(P0(i,:),Data(j,:));
  100.         end
  101.     end         
  102.     U=1./(Dist.^m.*(ones(C,1)*sum(Dist.^(-m))));      
  103.     % 更新聚类中心 P
  104.     Um=U.^M;
  105.     P=Um*Data./(ones(S,1)*sum(Um\'))\';   
  106.     % 目标函数值: 类内加权平方误差和
  107.     if nargout>4 | plotflag
  108.         Obj_Fcn(iter)=sum(sum(Um.*Dist.^2));
  109.     end
  110.     % FCM 算法迭代停止条件
  111.     if norm(P-P0,Inf)<epsm
  112.         break
  113.     end
  114.     P0=P;
  115. end
  116. % 聚类结果
  117. if nargout > 3
  118.     res = maxrowf(U);
  119.     for c = 1:C
  120.         v = find(res==c);
  121.         Cluster_Res(c,1:length(v))=v;
  122.     end
  123. end
  124. % 绘图
  125. if plotflag
  126.     fcmplot(Data,U,P,Obj_Fcn);
  127. end
  128. m文件3/7:
  129. function fcmplot(Data,U,P,Obj_Fcn)
  130. % FCM 结果绘图函数
  131. % See also: fuzzycm maxrowf ellipse
  132. [C,S] = size(P); res = maxrowf(U);
  133. str = \'po*x+d^v><.h\'; 
  134. % 目标函数绘图
  135. figure(1),plot(Obj_Fcn)
  136. title(\'目标函数值变化曲线\',\'fontsize\',8)
  137. % 2D 绘图
  138. if S==2 
  139.     figure(2),plot(P(:,1),P(:,2),\'rs\'),hold on
  140.     for i=1:C
  141.         v=Data(find(res==i),:); 
  142.         plot(v(:,1),v(:,2),str(rem(i,12)+1))      
  143.         ellipse(max(v(:,1))-min(v(:,1)), ...
  144.                 max(v(:,2))-min(v(:,2)), ...
  145.                 [max(v(:,1))+min(v(:,1)), ...
  146.                 max(v(:,2))+min(v(:,2))]/2,\'r:\')    
  147.     end
  148.     grid on,title(\'2D 聚类结果图\',\'fontsize\',8),hold off
  149. end
  150. % 3D 绘图
  151. if S>2 
  152.     figure(2),plot3(P(:,1),P(:,2),P(:,3),\'rs\'),hold on
  153.     for i=1:C
  154.         v=Data(find(res==i),:);
  155.         plot3(v(:,1),v(:,2),v(:,3),str(rem(i,12)+1))      
  156.         ellipse(max(v(:,1))-min(v(:,1)), ...
  157.                 max(v(:,2))-min(v(:,2)), ...
  158.                 [max(v(:,1))+min(v(:,1)), ...
  159.                 max(v(:,2))+min(v(:,2))]/2, ...
  160.                 \'r:\',(max(v(:,3))+min(v(:,3)))/2)   
  161.     end
  162.     grid on,title(\'3D 聚类结果图\',\'fontsize\',8),hold off
  163. end
  164. m文件4/7:
  165. function D=fuzzydist(A,B)
  166. % 模糊聚类分析: 样本间的距离
  167. % D = fuzzydist(A,B)
  168. D=norm(A-B);
  169. m文件5/7:
  170. function mr=maxrowf(U,c)
  171. % 求矩阵 U 每列第 c 大元素所在行,c 的缺省值为 1
  172. % 调用格式: mr = maxrowf(U,c)
  173. % See also: addr
  174. if nargin<2
  175.     c=1;
  176. end
  177. N=size(U,2);mr(1,N)=0;
  178. for j=1:N
  179.     aj=addr(U(:,j),\'descend\');
  180.     mr(j)=aj(c);
  181. end
  182. m文件6/7:
  183. function ellipse(a,b,center,style,c_3d)
  184. % 绘制一个椭圆
  185. % 调用: ellipse(a,b,center,style,c_3d)
  186. % 输入:
  187. %     a: 椭圆的轴长(平行于 x 轴)
  188. %     b: 椭圆的轴长(平行于 y 轴)
  189. %     center: 椭圆的中心 [x0,y0],缺省值为 [0,0]
  190. %     style: 绘制的线型和颜色,缺省值为实线蓝色
  191. %     c_3d:   椭圆的中心在 3D 空间中的 z 轴坐标,可缺省
  192. if nargin<4
  193.     style=\'b\';
  194. end
  195. if nargin<3 | isempty(center)
  196.     center=[0,0];
  197. end
  198. t=1:360;
  199. x=a/2*cosd(t)+center(1);
  200. y=b/2*sind(t)+center(2);
  201. if nargin>4
  202.     plot3(x,y,ones(1,360)*c_3d,style)
  203. else
  204.     plot(x,y,style)
  205. end
  206. m文件7/7:
  207. function f = addr(a,strsort)
  208. % 返回向量升序或降序排列后各分量在原始向量中的索引
  209. % 函数调用:f = addr(a,strsort)
  210. % strsort: \'ascend\' or \'descend\'
  211. %          default is \'ascend\'
  212. % -------- example --------
  213. % addr([ 4 5 1 2 ]) returns ans:
  214. %       [ 3   4   1   2 ]
  215. if nargin==1
  216.     strsort=\'ascend\';
  217. end
  218. sa=sort(a); ca=a;
  219. la=length(a);f(la)=0;
  220. for i=1:la
  221.     f(i)=find(ca==sa(i),1);
  222.     ca(f(i))=NaN;
  223. end
  224. if strcmp(strsort,\'descend\')
  225.     f=fliplr(f);
  226. end

相关文章