本实验的目的是学习和掌握k-均值聚类算法。k-均值算法是一种经典的无监督聚类和学习算法，它属于迭代优化算法的范畴。本实验在MATLAB平台上，编程实现了k-均值聚类算法，并使用20组三维数据进行测试，比较分类结果。实验中初始聚类中心由人为设定，以便于实验结果的比较与分析。
一、技术论述
1.无监督学习和聚类
在之前设计分类器的时候，通常需要事先对训练样本集的样本进行标定以确定类别归属。这种利用有标记样本集的方法称为“有监督”或“有教师”方法。这一类方法的使用固然十分广泛，也有着很坚实的理论基础，但在实际运用中这类方法经常会遇到以下瓶颈：
- 收集并标记大量样本集是一件相当费时费力的前期工作；
- 现实中存在很多应用，其分类模式的性质会随着时间发生变化，单单使用已标记样本无法满足这类情况；
- 有人希望能逆向解决问题：先用大量未标记的样本集来自动地训练分类器，再人工地标记数据分组的结果，如数据挖掘的大型应用，因为这些应用往往不知道待处理数据的具体情况。
可以看到，无监督学习方法的提出是十分必要的。事实上，在任何一项探索性的工作中，无监督的方法均向我们揭示了观测数据的一些普遍规律。很多无监督方法都可以以独立于数据的方式工作，为后续步骤提供“灵巧的预处理”和“灵巧的特征提取”等有效的前期处理。在无监督的情况下，聚类算法是模式识别研究中著名的一类技术。
2.分类与聚类的差别
分类(Classification)：对于一个分类器，通常需要你告诉它“这个东西被分为某某类”这样一些例子。通常情况下，一个分类器会从它得到的训练数据中进行学习，从而具备对未知数据进行分类的能力，这种提供训练数据的过程通常叫做有监督学习。
聚类(Clustering)：简单地说就是把相似的东西分到一组。聚类的时候，我们并不关心某一类是什么，这里需要实现的目标只是把相似的东西聚到一起。因此，一个聚类算法通常只需要知道如何计算相似度就可以开始工作了。因此，聚类方法通常并不需要使用训练数据进行学习，因此桔类方法属于无监督学习的范畴。
3.常见的分类与聚类算法
所谓分类，就是根据数据样本的特征或属性，划分到已有的类别中。前面使用到的模式分类方法主要有：贝叶斯分类算法(Bayesian classifier) 、PCA主分量分析法、Fisher线性判别分析法、Parzen窗估计法、k-最近邻法(k-nearest neighbor，kNN)、基于支持向量机(SVM)的分类器、人工神经网络(ANN)和决策树分类法等等。
分类作为一种有监督学习方法，要求必须在分类之前明确知道各个类别的必要特征和信息，并且标记所有训练样本都有一个类别与之相对应。但是很多情况下这些条件往往无法满足，尤其是在处理海量数据的时候，数据预处理的代价非常大。
聚类算法中最经典的当属k均值聚类(K-means clustering)算法。该算法又称为“c均值算法”，因为它的目标就是找到c个均值向量作为聚类中心：μ1,μ2,…,μc，实际上k与c是等价的。
以上是对二维随机样本进行聚类的实例。
4.聚类任务的基本步骤
假设所有模式都用一组特征表示，模式被表示为L维的特征矢量。聚类任务需包含以下步骤：
- 特征选择。选择特征尽可能与感兴趣的任务相关。特征之间的信息冗余度要尽可能小。
- 相似性测度。一个基本的保证是所有选择的特征对相似性测度计算的贡献都是均衡的，没有那一个特征是绝对占优的。
- 聚类准则。聚类准则依赖于专家判断的在数据集合内部隐含类的类型解释。聚类准则可以被表示为代价函数和其它类型的规则。
- 聚类算法。在确定相似性测度和聚类准则后，这一步就是要选择一个具体的算法方案将数据集合分解为类结构。
- 结构的有效性。一旦聚类算法获得了结果，需要采用合适的检验方法检验其正确性。
- 结果的解释。应用领域的专家必须结合其它的试验证据和分析解释聚类结果，以便得到正确的结论。
5.k-均值聚类算法
k-均值聚类算法的目标是找到k个均值向量或“聚类中心”。算法的实现步骤如下所示，其中n表示模式的数量，c表示类别的数量，通常的做法是从样本中随机取出c个作为初始的聚类中心。当然，初始的聚类中心也可以通过人为来确定：
该算法的计算复杂度为：
其中d代表样本的维数，T是聚类的迭代次数，一般来说，迭代次数通常远少于样本的数量。
该算法是一种典型的聚类算法，把它归入迭代优化算法的范畴是因为算法规定的c个均值会不断地移动，使得一个平方误差准则函数最小化。在算法的每一步迭代中，每个样本点均被认为是完全属于某一类别。
二、实验结果讨论
实验所使用的样本：
设计步骤主要包括以下几个部分：
编写程序，实现以上所述的k-均值聚类算法。其中，在算法中样本与聚类中心的距离采用欧氏距离的形式。
类别数目和聚类中心初始值选为以下参数进行实验：
再将类别数目和聚类中心初始值改变为以下参数进行实验。
下图得到两次聚类的结果，可以看到当分类类别为2时，初始聚类中心对分类结果影响不大（至少对于样本少的情况是这样的），最终两种情况都能得到一样的最终聚类中心。
下面测试将样本分为三类的情况。将类别数目和聚类中心初始值选为以下参数进行实验：
再将类别数目和聚类中心初始值改变为以下参数进行实验：
下图得到两次聚类的结果，可以看到当分类类别为3时，分类复杂度增加，随着聚类中心的移动，对于同一组测试样本可能有不同的划分结果。虽然初始聚类中心的作用只是将样本初步地分为几个区域，但事实上不同的初始中心会给分类结果带来巨大的差异。
在程序中，使用欧氏距离作为样本到聚类中心的距离，事实上也可以使用其他多种距离度量进行运算，如街区距离（1范数）、棋盘距离（∞范数）等等。
三、完整代码
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% k-均值聚类函数
% 输入参数：
% w：需要分类的样本
% k：分类数
% m：初始聚类中心
% iteration：迭代次数
% 中间参数：
% class_id：存放各个样本属于一类的下标
% dist：计算样本到聚类中心的距离
% 输出参数：
% class_result：存放样本的分类结果
% class_num：存放被分到各类的样本个数
% center：迭代结束时的聚类中心
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [class_result, class_num, center] = kmeans(w, k, m, iteration)
[n,d] = size(w);
class_result = zeros(1,n);
class_num = zeros(1,k);
time = 1;
% 以下步骤计算每个样本到聚类中心的距离
while time < iteration % 迭代次数限制
for i = 1:n
dist = sqrt(sum((repmat(w(i,:), k, 1) - m).^2, 2)); % 欧氏距离
% dist = sum(abs(repmat(x(i,:), k, 1) - nc), 2); % 街区距离
[y,class_id] = min(dist); % 计算样本对三类中哪一类有最小距离并存放在class_id
class_result(i) = class_id;
end
for i = 1:k
% 找到每一类的所有数据，计算平均值，其值作为新的聚类中心
class_id = find(class_result == i);
m(i, :) = mean(w(class_id, :)); % 更新聚类中心
% 统计每一类样本的个数
class_num(i) = length(class_id);
end
time = time + 1;
end
center = m;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 样本分类结果的绘图函数
% 输入参数：
% w：需要分类的样本
% class：聚类后的样本分类结果
% flag：分类类别数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function draw(w, class, center)
[flag x] = size(center);
[n, d] = size(w);
% figure;
if flag == 2 % 若将样本分成两类
for i=1:n
if class(i) == 1
plot3(w(i,1),w(i,2),w(i,3),"r+"); % 显示第一类
hold on;
grid on;
elseif class(i) == 2
plot3(w(i,1),w(i,2),w(i,3),"go"); %显示第二类
hold on;
grid on;
end
end
for j = 1:flag
if j == 1
plot3(center(j,1),center(j,2),center(j,3),"rd"); % 聚类中心
elseif j == 2
plot3(center(j,1),center(j,2),center(j,3),"gd"); % 聚类中心
end
end
end
if flag == 3 % 若将样本分成三类
% 显示分类结果
for i = 1:n
if class(i) == 1
plot3(w(i,1),w(i,2),w(i,3),"r+"); % 显示第一类
hold on;
grid on;
elseif class(i) == 2
plot3(w(i,1),w(i,2),w(i,3),"go"); % 显示第二类
hold on;
grid on;
elseif class(i) == 3
plot3(w(i,1),w(i,2),w(i,3),"b*"); % 显示第三类
hold on;
grid on;
end
end
for j = 1:flag
if j == 1
plot3(center(j,1),center(j,2),center(j,3),"rd"); % 聚类中心
elseif j == 2
plot3(center(j,1),center(j,2),center(j,3),"gd"); % 聚类中心
elseif j == 3
plot3(center(j,1),center(j,2),center(j,3),"bd"); % 聚类中心
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% k-均值聚类的研究与实现主函数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
close all;
% 测试样本
w = [-7.82 -4.58 -3.97;...
-6.68 3.16 2.71;...
4.36 -2.19 2.09;...
6.72 0.88 2.80;...
-8.64 3.06 3.50;...
-6.87 0.57 -5.45;...
4.47 -2.62 5.76;...
6.73 -2.01 4.18;...
-7.71 2.34 -6.33;...
-6.91 -0.49 -5.68;...
6.18 2.18 5.28;...
6.72 -0.93 -4.04;...
-6.25 -0.26 0.56;...
-6.94 -1.22 1.33;...
8.09 0.20 2.25;...
6.81 0.17 -4.15;...
-5.19 4.24 4.04;...
-6.38 -1.74 1.43;...
4.08 1.30 5.33;...
6.27 0.93 -2.78];
[n, d] = size(w);
% 以下是k均值聚类的参数设定（c = 2时）
k = 2;
m = [1 1 1; -1 1 -1]; % 初始聚类中心
% m = [0 0 0; 1 1 -1]; % 初始聚类中心
iteration = 200; % k均值聚类的迭代次数
% 调用kmeans函数进行聚类
[class, class_num, center] = kmeans(w, k, m, iteration);
% 画出样本分类结果
figure;subplot(121);draw(w, class, center);
title("使用第一种初始聚类中心时，k-均值聚类算法分类结果");
% 显示信息
disp(["属于第一类的样本个数为：",num2str(class_num(1))]);
disp(["属于第二类的样本个数为：",num2str(class_num(2))]);
disp("最终的聚类中心为：");
disp(num2str(center));
% 以下是k均值聚类的参数设定（c = 2时）
k = 2;
% m = [1 1 1; -1 1 -1]; % 初始聚类中心
m = [0 0 0; 1 1 -1]; % 初始聚类中心
iteration = 200; % k均值聚类的迭代次数
% 调用kmeans函数进行聚类
[class, class_num, center] = kmeans(w, k, m, iteration);
% 画出样本分类结果
subplot(122);draw(w, class, center);
title("使用第二种初始聚类中心时，k-均值聚类算法分类结果");
% 显示信息
disp(["属于第一类的样本个数为：",num2str(class_num(1))]);
disp(["属于第二类的样本个数为：",num2str(class_num(2))]);
disp("最终的聚类中心为：");
disp(num2str(center));
% 以下是k均值聚类的参数设定（c = 3时）
k = 3;
m = [0 0 0; 1 1 1; -1 0 2]; % 初始聚类中心
% m = [-0.1 0.0 0.1; 0 -0.1 0.1; -0.1 -0.1 0.1]; % 初始聚类中心
iteration = 200; % k均值聚类的迭代次数
% 调用kmeans函数进行聚类
[class, class_num, center] = kmeans(w, k, m, iteration);
% 画出样本分类结果
figure;subplot(121);draw(w, class, center);
title("使用第一种初始聚类中心时，k-均值聚类算法分类结果");
% 显示信息
disp(["属于第一类的样本个数为：",num2str(class_num(1))]);
disp(["属于第二类的样本个数为：",num2str(class_num(2))]);
disp(["属于第三类的样本个数为：",num2str(class_num(3))]);
disp("最终的聚类中心为：");
disp(num2str(center));
% 以下是k均值聚类的参数设定（c = 3时）
k = 3;
% m = [0 0 0; 1 1 1; -1 0 2]; % 初始聚类中心
m = [-0.1 0.0 0.1; 0 -0.1 0.1; -0.1 -0.1 0.1]; % 初始聚类中心
iteration = 200; % k均值聚类的迭代次数
% 调用kmeans函数进行聚类
[class, class_num, center] = kmeans(w, k, m, iteration);
% 画出样本分类结果
subplot(122);draw(w, class, center);
title("使用第二种初始聚类中心时，k-均值聚类算法分类结果");
% 显示信息
disp(["属于第一类的样本个数为：",num2str(class_num(1))]);
disp(["属于第二类的样本个数为：",num2str(class_num(2))]);
disp(["属于第三类的样本个数为：",num2str(class_num(3))]);
disp("最终的聚类中心为：");
disp(num2str(center));