我的K均值算法的matlab实现

时间：2022-05-15 22:39:40

这是我的第一篇博客；


K-Means算法过程，略；

这是一次课程的任务2333，是利用所学K-means聚类分析方法，对iris数据集进行聚类分析，并利用已知的样本类别标
签进行聚类分析评价；

我的K均值算法以iris.data为例（附在文末）；

数据集：Iris数据集
 (http://archive.ics.uci.edu/ml/datasets/Iris)
数据描述：Iris数据集包含150个鸢尾花模式样本，其中 每个模式样本采用5维的特征描述
X = (x1,x2,x3,x4,w);
x1: 萼片长度(厘米）
x2: 萼片宽度(厘米）
x3: 花瓣长度(厘米）
x4：花瓣宽度(厘米）
w(类别属性 ): 山鸢尾 (Iris setosa)，变色鸢尾(Iris versicolor)和维吉尼亚鸢尾(Iris virginica)


先贴上我的函数结构：

函数结构

	    —— FindCluster(~)		 聚类算法主函数
            |
MyKmeans ——  MyPlot2(~),MyPlot3(~)     画图
	    |
	    —— Accuracy(~)             聚类精度评价

MyKmean.m，程序运行的主函数

% 作者 ： YG1501  LQY
% 日期 ： 2018.01.13  星期六
% 函数功能 ： 实现对iris.data数据集的分类，并根据分类结果进行精度评价

clear;clc;close all;
%手动选取文件，选用iris.data
[filename,fpath] = uigetfile(...
{...
	'*.m;*.txt;*.data',...
		'Data(*.m;*.txt;*.data)';...
			'*.*','All Files(*.*)'...
},'Select data');

[X1,X2,X3,X4,X5] = ...
	textread(filename,'%f%f%f%f%s','delimiter',','); 	%Get data

X = [X1 X2 X3 X4];	%目前貌似还没什么用
[m,n] = size(X);

DataLabel = zeros(m,1);
for i = 1:size(X5,1)
	if (strcmp(X5(i),'Iris-setosa'))
		DataLabel(i) = 1 ;
	elseif(strcmp(X5(i),'Iris-versicolor'))
		DataLabel(i) = 2 ;
	else
		DataLabel(i) = 3 ;
	end
end


%二维结果
[MyCenter1,ClusterLabel12] = FindCluster([X1 X2],3,DataLabel);
[MyCenter2,ClusterLabel13] = FindCluster([X1 X3],3,DataLabel);
[MyCenter3,ClusterLabel14] = FindCluster([X1 X4],3,DataLabel);
[MyCenter4,ClusterLabel23] = FindCluster([X2 X3],3,DataLabel);
[MyCenter5,ClusterLabel24] = FindCluster([X2 X4],3,DataLabel);
[MyCenter6,ClusterLabel34] = FindCluster([X3 X4],3,DataLabel);

hold on;
figure(1),title('2D');
subplot(231)
MyPlot2(X1,X2,DataLabel,MyCenter1),xlabel('X1'),ylabel('X2')
subplot(232)
MyPlot2(X1,X3,DataLabel,MyCenter2),xlabel('X1'),ylabel('X3')
subplot(233)
MyPlot2(X1,X4,DataLabel,MyCenter3),xlabel('X1'),ylabel('X4')
subplot(234)
MyPlot2(X2,X3,DataLabel,MyCenter4),xlabel('X2'),ylabel('X3')
subplot(235)
MyPlot2(X2,X4,DataLabel,MyCenter5),xlabel('X2'),ylabel('X4')
subplot(236)
MyPlot2(X3,X4,DataLabel,MyCenter6),xlabel('X3'),ylabel('X4')


%三维结果
% [MyCenter7,ClusterLabel123]  = FindCluster([X1,X2,X3],3,DataLabel);
% [MyCenter8,ClusterLabel124]  = FindCluster([X1,X2,X4],3,DataLabel);
% [MyCenter9,ClusterLabel134]  = FindCluster([X1,X3,X4],3,DataLabel);
% [MyCenter10,ClusterLabel234] = FindCluster([X2,X3,X4],3,DataLabel);

% hold on;
% figure(2),title('3D');
% subplot(221)
% MyPlot3(X1,X2,X3,DataLabel,MyCenter7)
% xlabel('X1'),ylabel('X2'),zlabel('X3');
% subplot(222)
% MyPlot3(X1,X2,X4,DataLabel,MyCenter8)
% xlabel('X1'),ylabel('X2'),zlabel('X4');
% subplot(223)
% MyPlot3(X1,X3,X4,DataLabel,MyCenter9)
% xlabel('X1'),ylabel('X3'),zlabel('X4');
% subplot(224)
% MyPlot3(X2,X3,X4,DataLabel,MyCenter10)
% xlabel('X2'),ylabel('X3'),zlabel('X4');



%聚类精度评价
%二维结果
ClusterLabel = [ClusterLabel12,ClusterLabel13,ClusterLabel14...
				ClusterLabel23,ClusterLabel24,ClusterLabel34]
ClusterAccuracy = Accuracy(DataLabel,ClusterLabel)  % 输出每种分类的精度

%三维结果
%略
这里要说明的一点是，我的函数的运行效率比较低，建议运行时只单独进行二维的分类或三维的分类。。运算量比较多，在二维空间上大概需要5秒钟才能出结果。。。


FindCluster.m

%函数功能 ： 输入数据集、聚类中心个数与样本标签
%	 得到聚类中心与聚类样本标签

function [ClusterCenter,ClusterLabel] = FindCluster(MyData,ClusterCounts,DataLabel)
[m,n] = size(MyData);

ClusterLabel = zeros(m,1);		%用于存储聚类标签

% MyLabel = unique(DataLabel,'rows');
% for i = 1:size(MyLabel,2);	
% 	LabelIndex(1,i) = i;		%为数据标签创建索引
% end

%已知数据集的每个样本的中心
OriginCenter = zeros(ClusterCounts,n);
for j = 1:ClusterCounts
	DataCounts = 0;
	for i = 1:m
		%按照数据标签，计算样本中心
		if DataLabel(i) == j
			OriginCenter(j,:) = OriginCenter(j,:) + MyData(i,:);
			DataCounts = DataCounts + 1;
		end
	end
	OriginCenter(j,:) = OriginCenter(j,:) ./ DataCounts;
end
%按照第一列对样本中心排序
SortCenter1 = sortrows(OriginCenter,1);

FalseTimes = 0;
CalcuateTimes = 0;
%此循环用于纠正分类错误的情况
while (CalcuateTimes < 15)
	ClusterCenter = zeros(ClusterCounts,n);		%初始化聚类中心
	for i = 1:ClusterCounts
		ClusterCenter(i,:) = MyData( randi(m,1),:);	%随机选取一个点作为中心
	end

	%此循环用于寻找聚类中心
	%目前还未解决该循环陷入死循环的问题，所以设置一个参数来终止循环
	kk = 0;
	while (kk < 15)
		Distance   = zeros(1,ClusterCounts);	%存储单个样本到每个聚类中心的距离
		DataCounts = zeros(1,ClusterCounts);	%记录每个聚类的样本数目
		NewCenter  = zeros(ClusterCounts,n);
		for i = 1:m
			for j = 1:ClusterCounts
				Distance(j) = norm(MyData(i,:) - ClusterCenter(j,:));
			end
			%index返回最小距离的索引，亦即聚类中心的标号
			[~,index] = min(Distance);
			ClusterLabel(i) = index;
		end 

		k = 0;
		for j = 1:ClusterCounts
			for i = 1:m
				%按照标记，对样本进行分类,并计算聚类中心
				if ClusterLabel(i) == j
					NewCenter(j,:) = NewCenter(j,:) + MyData(i,:);
					DataCounts(j)  = DataCounts(j) + 1;
				end
			end
			NewCenter(j,:) = NewCenter(j,:) ./ DataCounts(j);
			%若新的聚类中心与上一个聚类中心相同，则认为算法收敛
			if norm(NewCenter(j,:) - ClusterCenter(j,:)) < 0.1
				k = k + 1;
			end
		end

		ClusterCenter = NewCenter;
		%判断是否全部收敛
		if k == ClusterCounts
			break;
		end
		kk = kk + 1 ;
	end

	%再次判断每个聚类是否分类正确，若分类错误则进行惩罚
	t = ClusterCounts;
	SortCenter2 = sortrows(ClusterCenter,1);
	for i = 1:ClusterCounts
		if norm(SortCenter1(i,:) - SortCenter2(i,:)) > 0.5
			t = t - 1;
			FalseTimes = FalseTimes + 1;
			break;
		end
	end

	if t == ClusterCounts
		break;
	end
	CalcuateTimes = CalcuateTimes + 1;
end

% FalseTimes
% CalcuateTimes
% t

% kk
% DataCounts
% OriginCenter
% NewCenter

%理论上每个聚类的标签应是123排列的，但实际上，由于每个聚类中心都是随机选取的，
%实际分类的顺序可能是213，132等，所以需要对分类标签进行纠正，这对之后的精度评
%价带来了方便。
%对分类标签进行纠正：
%算法原理：从第一个已知的样本中心开始，寻找离其最近的聚类中心，然后将归类于该
%		   聚类中心的样本的聚类标签更换为i
for i = 1:ClusterCounts
	for j = 1:ClusterCounts
		if norm(OriginCenter(i,:) - ClusterCenter(j,:)) < 0.6
			for p = 1:m
				if ClusterLabel(p) == j
					ClusterLabel(p) = 2 * ClusterCounts + i;
				end
			end
			%break;
		end
	end
end
ClusterLabel = ClusterLabel - 2 * ClusterCounts;
%Temp = [MyData(:,:),ClusterLabel]

至此已经完成了聚类中心的计算。该方法基本解决了因随机选取初始中心而导致最后聚类中心明显错误的情况，但缺点在于循环太多，时间复杂度O（？），而且目前还未解决偶尔会进入死循环的问题。若判明进入了死循环，还是重新运行下程序吧 = = 。。。


画图的函数（- -）：

MyPlot2.m

%函数功能 : 输入样本，样本标签及求出的聚类中心，显示二维图像，实现数据可视化

function MyPlot2(X1,X2,DataLabel,ClusterCenter)
[m,~] = size(X1);

hold on;
for i = 1:m
	if(DataLabel(i) == 1)
		plot(X1(i),X2(i),'*r')
	elseif(DataLabel(i) == 2)
		plot(X1(i),X2(i),'*g')
	else
		plot(X1(i),X2(i),'*b')
	end
end

% xlabel(who('X1'));
% ylabel(who('X2'));
%PS ： 我想在坐标轴中根据输入的形参的矩阵的名字，转换成字符串，来动态输出
%	   坐标轴名称，不知道该怎么做？用上面注释的语句不行。。
[n,~] = size(ClusterCenter);

for i = 1:n
	plot(ClusterCenter(i,1),ClusterCenter(i,2),'ok')
end

grid on;
MyPlot3.m 

function MyPlot3(X1,X2,X3,DataLabels,ClusterCenter)
[m,~] = size(X1);

hold on;
for i = 1:m
	if(DataLabels(i) == 1)
		plot3(X1(i),X2(i),X3(i),'.r')
	elseif(DataLabels(i) == 2)
		plot3(X1(i),X2(i),X3(i),'.g')
	else
		plot3(X1(i),X2(i),X3(i),'.b')
	end
end

% xlabel('X1');
% ylabel('X2');
% zlabel('X3');

[n,~] = size(ClusterCenter);

for i = 1:n
	plot3(ClusterCenter(i,1),ClusterCenter(i,2),ClusterCenter(i,3),'ok')
end

view([1 1 1]);
grid on;
Accuracy.m 聚类结果精度评价 

%函数功能：根据聚类结果进行精度评价
%精度评价，返回每一种分类的精度值(正确率)
function ClusterAccuracy = Accuracy(DataLable,CLusterLabel)
[m,n] = size(CLusterLabel);
ClusterAccuracy = zeros(1,n);

%理论上的聚类标签应为 1,2,3，但实际上可能变成了 213 , 132等，导致计算失误
%因此需要对分类标签进行纠正，而这一步骤已经在FindCluster函数中完成了
for i = 1:n
	%原理：假设某样本在已知数据集中属于第一类，而其聚类后的也同样被分到了第一类，那么它们的标签
	%都是1，这样相减后结果就为0，表明已经分类正确，否则不为0，分类错误
	Temp(:,i) = DataLable - CLusterLabel(:,i);
end

for j = 1:n
	for i = 1:m
		if Temp(i,j) == 0
			ClusterAccuracy(1,j) = ClusterAccuracy(1,j) + 1;
		end
	end
end

ClusterAccuracy = ClusterAccuracy ./ m;


运行结果展示（二维结果）：
我的K均值算法的matlab实现
我的K均值算法的matlab实现
可以看到，在二维的情况下，比较X 3: 花瓣长度(厘米）和X 4：花瓣宽度(厘米）的精度更高（94.67%），也就是说只比较两种特征时，比较 花瓣长度和花瓣宽度，区分三种花的效果更好，分类结果更可靠。


三维分类结果：

我的K均值算法的matlab实现我的K均值算法的matlab实现

可以看到，在三维的情况下，比较X2：萼片宽度(厘米），X3: 花瓣长度(厘米）和X4：花瓣宽度(厘米）的精度更高（93.33%），也就是说比较三种特征时，取X2，X3，X4，区分三种花的效果更好，分类结果更可靠。



最后：本程序目前还存在诸多不足，比如时间复杂度高，效率较低；目前对Matlab语言还不是很熟，写的程序也比较C-Style，各位看官若是有改进的建议，欢迎留言，或者直接联系我（联系方式附在文最末），大家一起探讨:D（手动比心）……



附iris.data数据集：

5.1,3.5,1.4,0.2,Iris-setosa
4.9,3.0,1.4,0.2,Iris-setosa
4.7,3.2,1.3,0.2,Iris-setosa
4.6,3.1,1.5,0.2,Iris-setosa
5.0,3.6,1.4,0.2,Iris-setosa
5.4,3.9,1.7,0.4,Iris-setosa
4.6,3.4,1.4,0.3,Iris-setosa
5.0,3.4,1.5,0.2,Iris-setosa
4.4,2.9,1.4,0.2,Iris-setosa
4.9,3.1,1.5,0.1,Iris-setosa
5.4,3.7,1.5,0.2,Iris-setosa
4.8,3.4,1.6,0.2,Iris-setosa
4.8,3.0,1.4,0.1,Iris-setosa
4.3,3.0,1.1,0.1,Iris-setosa
5.8,4.0,1.2,0.2,Iris-setosa
5.7,4.4,1.5,0.4,Iris-setosa
5.4,3.9,1.3,0.4,Iris-setosa
5.1,3.5,1.4,0.3,Iris-setosa
5.7,3.8,1.7,0.3,Iris-setosa
5.1,3.8,1.5,0.3,Iris-setosa
5.4,3.4,1.7,0.2,Iris-setosa
5.1,3.7,1.5,0.4,Iris-setosa
4.6,3.6,1.0,0.2,Iris-setosa
5.1,3.3,1.7,0.5,Iris-setosa
4.8,3.4,1.9,0.2,Iris-setosa
5.0,3.0,1.6,0.2,Iris-setosa
5.0,3.4,1.6,0.4,Iris-setosa
5.2,3.5,1.5,0.2,Iris-setosa
5.2,3.4,1.4,0.2,Iris-setosa
4.7,3.2,1.6,0.2,Iris-setosa
4.8,3.1,1.6,0.2,Iris-setosa
5.4,3.4,1.5,0.4,Iris-setosa
5.2,4.1,1.5,0.1,Iris-setosa
5.5,4.2,1.4,0.2,Iris-setosa
4.9,3.1,1.5,0.1,Iris-setosa
5.0,3.2,1.2,0.2,Iris-setosa
5.5,3.5,1.3,0.2,Iris-setosa
4.9,3.1,1.5,0.1,Iris-setosa
4.4,3.0,1.3,0.2,Iris-setosa
5.1,3.4,1.5,0.2,Iris-setosa
5.0,3.5,1.3,0.3,Iris-setosa
4.5,2.3,1.3,0.3,Iris-setosa
4.4,3.2,1.3,0.2,Iris-setosa
5.0,3.5,1.6,0.6,Iris-setosa
5.1,3.8,1.9,0.4,Iris-setosa
4.8,3.0,1.4,0.3,Iris-setosa
5.1,3.8,1.6,0.2,Iris-setosa
4.6,3.2,1.4,0.2,Iris-setosa
5.3,3.7,1.5,0.2,Iris-setosa
5.0,3.3,1.4,0.2,Iris-setosa
7.0,3.2,4.7,1.4,Iris-versicolor
6.4,3.2,4.5,1.5,Iris-versicolor
6.9,3.1,4.9,1.5,Iris-versicolor
5.5,2.3,4.0,1.3,Iris-versicolor
6.5,2.8,4.6,1.5,Iris-versicolor
5.7,2.8,4.5,1.3,Iris-versicolor
6.3,3.3,4.7,1.6,Iris-versicolor
4.9,2.4,3.3,1.0,Iris-versicolor
6.6,2.9,4.6,1.3,Iris-versicolor
5.2,2.7,3.9,1.4,Iris-versicolor
5.0,2.0,3.5,1.0,Iris-versicolor
5.9,3.0,4.2,1.5,Iris-versicolor
6.0,2.2,4.0,1.0,Iris-versicolor
6.1,2.9,4.7,1.4,Iris-versicolor
5.6,2.9,3.6,1.3,Iris-versicolor
6.7,3.1,4.4,1.4,Iris-versicolor
5.6,3.0,4.5,1.5,Iris-versicolor
5.8,2.7,4.1,1.0,Iris-versicolor
6.2,2.2,4.5,1.5,Iris-versicolor
5.6,2.5,3.9,1.1,Iris-versicolor
5.9,3.2,4.8,1.8,Iris-versicolor
6.1,2.8,4.0,1.3,Iris-versicolor
6.3,2.5,4.9,1.5,Iris-versicolor
6.1,2.8,4.7,1.2,Iris-versicolor
6.4,2.9,4.3,1.3,Iris-versicolor
6.6,3.0,4.4,1.4,Iris-versicolor
6.8,2.8,4.8,1.4,Iris-versicolor
6.7,3.0,5.0,1.7,Iris-versicolor
6.0,2.9,4.5,1.5,Iris-versicolor
5.7,2.6,3.5,1.0,Iris-versicolor
5.5,2.4,3.8,1.1,Iris-versicolor
5.5,2.4,3.7,1.0,Iris-versicolor
5.8,2.7,3.9,1.2,Iris-versicolor
6.0,2.7,5.1,1.6,Iris-versicolor
5.4,3.0,4.5,1.5,Iris-versicolor
6.0,3.4,4.5,1.6,Iris-versicolor
6.7,3.1,4.7,1.5,Iris-versicolor
6.3,2.3,4.4,1.3,Iris-versicolor
5.6,3.0,4.1,1.3,Iris-versicolor
5.5,2.5,4.0,1.3,Iris-versicolor
5.5,2.6,4.4,1.2,Iris-versicolor
6.1,3.0,4.6,1.4,Iris-versicolor
5.8,2.6,4.0,1.2,Iris-versicolor
5.0,2.3,3.3,1.0,Iris-versicolor
5.6,2.7,4.2,1.3,Iris-versicolor
5.7,3.0,4.2,1.2,Iris-versicolor
5.7,2.9,4.2,1.3,Iris-versicolor
6.2,2.9,4.3,1.3,Iris-versicolor
5.1,2.5,3.0,1.1,Iris-versicolor
5.7,2.8,4.1,1.3,Iris-versicolor
6.3,3.3,6.0,2.5,Iris-virginica
5.8,2.7,5.1,1.9,Iris-virginica
7.1,3.0,5.9,2.1,Iris-virginica
6.3,2.9,5.6,1.8,Iris-virginica
6.5,3.0,5.8,2.2,Iris-virginica
7.6,3.0,6.6,2.1,Iris-virginica
4.9,2.5,4.5,1.7,Iris-virginica
7.3,2.9,6.3,1.8,Iris-virginica
6.7,2.5,5.8,1.8,Iris-virginica
7.2,3.6,6.1,2.5,Iris-virginica
6.5,3.2,5.1,2.0,Iris-virginica
6.4,2.7,5.3,1.9,Iris-virginica
6.8,3.0,5.5,2.1,Iris-virginica
5.7,2.5,5.0,2.0,Iris-virginica
5.8,2.8,5.1,2.4,Iris-virginica
6.4,3.2,5.3,2.3,Iris-virginica
6.5,3.0,5.5,1.8,Iris-virginica
7.7,3.8,6.7,2.2,Iris-virginica
7.7,2.6,6.9,2.3,Iris-virginica
6.0,2.2,5.0,1.5,Iris-virginica
6.9,3.2,5.7,2.3,Iris-virginica
5.6,2.8,4.9,2.0,Iris-virginica
7.7,2.8,6.7,2.0,Iris-virginica
6.3,2.7,4.9,1.8,Iris-virginica
6.7,3.3,5.7,2.1,Iris-virginica
7.2,3.2,6.0,1.8,Iris-virginica
6.2,2.8,4.8,1.8,Iris-virginica
6.1,3.0,4.9,1.8,Iris-virginica
6.4,2.8,5.6,2.1,Iris-virginica
7.2,3.0,5.8,1.6,Iris-virginica
7.4,2.8,6.1,1.9,Iris-virginica
7.9,3.8,6.4,2.0,Iris-virginica
6.4,2.8,5.6,2.2,Iris-virginica
6.3,2.8,5.1,1.5,Iris-virginica
6.1,2.6,5.6,1.4,Iris-virginica
7.7,3.0,6.1,2.3,Iris-virginica
6.3,3.4,5.6,2.4,Iris-virginica
6.4,3.1,5.5,1.8,Iris-virginica
6.0,3.0,4.8,1.8,Iris-virginica
6.9,3.1,5.4,2.1,Iris-virginica
6.7,3.1,5.6,2.4,Iris-virginica
6.9,3.1,5.1,2.3,Iris-virginica
5.8,2.7,5.1,1.9,Iris-virginica
6.8,3.2,5.9,2.3,Iris-virginica
6.7,3.3,5.7,2.5,Iris-virginica
6.7,3.0,5.2,2.3,Iris-virginica
6.3,2.5,5.0,1.9,Iris-virginica
6.5,3.0,5.2,2.0,Iris-virginica
6.2,3.4,5.4,2.3,Iris-virginica
5.9,3.0,5.1,1.8,Iris-virginica


博主联系方式 ： QQ：1765928683

邮箱：kelrays@qq.com


2018 / 01 / 27   补充：

我想到了一些提高运行效率的方法，那就是尽量减少for循环，因为matlab对for循环的处理效率是很低的！

可拱参考的优化思路：

1） 使用parfor；

2）使用find；

3）向量化，即使用向量操作代替循环;

4）bsxfun（），sum（）， ‘./’ '.*'

5）待补充


2018 / 05 / 03   更新：

程序之所以跑得慢，是因为画图函数用了for循环，把for循环去掉直接用plot，会快很多








标签：

相关文章