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本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。
