摘要：本文详细介绍Kmeans聚类算法的原理和程序实现。首先介绍利用该算法的原理及理解，详细介绍基于MATLAB设计一个自定义的Kmeans函数过程，然后利用该函数对UCI的数据集进行聚类以测试聚类结果。后续章节将介绍的主要部分有：

Kmeans算法的原理及理解
编程实现
聚类结果评价
类簇中心点的选取

点击下载：本文Kmeans算法M函数及测试完整文件
 

1. 前言

作为无监督聚类算法中的代表——K均值聚类（Kmeans）算法，该算法的主要作用是将相似的样本自动归到一个类别中。所谓的监督算法，就是输入样本没有对应的输出或标签。聚类（clustering）试图将数据集中的样本划分为若干个通常是不相交的子集，每个子集称为一个“簇（cluster）”，聚类既能作为一个单独过程，用于找寻数据内在的分布结构，也可作为分类等其他学习任务的前去过程。——《Machine Learning》

聚类算法也许是机器学习中“新算法”出现最多、最快的领域，一个重要的原因是聚类不存在客观标准，给定数据集总能从某个角度找到以往算法未覆盖的某种标准从而设计出新算法。Kmeans算法十分简单易懂而且非常有效，但是合理的确定K值和K个初始类簇中心点对于聚类效果的好坏有很大的影响。众多的论文基于此都提出了各自行之有效的解决方案，新的改进算法仍然不断被提出，此类文章大家可以在Web Of Science中搜索。
尽管Kmeans算法在MATLAB、Python等语言的工具箱函数中都有自带的函数可供调用，但作为机器学习的研究者新来说要设计出新的算法，有时就得“定制”自己的Kmeans函数了。自己动手编写无疑也更加能理解算法的具体过程，接下来就让我们进入正题吧

2. Kmeans算法的原理与理解
Kmeans算法是最常用的聚类算法，主要思想是:在给定K值和K个初始类簇中心点的情况下，把每个点(亦即数据记录)分到离其最近的类簇中心点所代表的类簇中，所有点分配完毕之后，根据一个类簇内的所有点重新计算该类簇的中心点(取平均值)，然后再迭代的进行分配点和更新类簇中心点的步骤，直至类簇中心点的变化很小，或者达到指定的迭代次数。
2.1 基本原理
假定给定数据样本X，包含了n个对象$X=\left \{ X_{1} ,X_{2},X_{3},...,X_{n}\right \}$,其中每个对象都具有m个维度的属性。Kmeans算法的目标是将n个对象依据对象间的相似性聚集到指定的k个类簇中，每个对象属于且仅属于一个其到类簇中心距离最小的类簇中。对于Kmeans，首先需要初始化k个聚类中心$\left \{ C_{1} ,C_{2},C_{3},...,C_{k}\right \},1<k\leq n$,然后通过计算每一个对象到每一个聚类中心的欧式距离，如下式所示$dis(X_{i},C_{j})=\sqrt{\sum_{t=1}^{m}(X_{it}-C_{jt})^{2}}$

上式中，$X_{i}$表示第i个对象$1\leq i\leq n$,$C_{j}$表示第j个聚类中心的$1\leq j\leq k$，$X_{it}$表示第i个对象的第t个属性，$1\leq t\leq m$，$C_{jt}$表示第j个聚类中心的第t个属性。
依次比较每一个对象到每一个聚类中心的距离，将对象分配到距离最近的聚类中心的类簇中，得到k个类簇$\left \{ S_{1},S_{2},S_{3},...,S_{k} \right \}$


Kmeans算法用中心定义了类簇的原型，类簇中心就是类簇内所有对象在各个维度的均值，其计算公式如下$C_{t}=\frac{\sum_{X_{i}\in S_{l}}X_{i}}{\left | S_{l} \right |}$

式中，$C_{l}$表示第l个聚类的中心，$1\leq l\leq k$，$\left | S_{l} \right |$表示第l个类簇中对象的个数，$X_{i}$表示第l个类簇中第i个对象，$1\leq i\leq\left | S_{l} \right |$。
2.2 算法流程

输入：样本集$D=\left \{ x_{1},x_{2},x_{3},...,x_{m} \right \}$；聚类簇数k.

过程：

1：从D中随机选择k个样本作为初始均值向量$\left \{ \mu _{1},\mu _{2},\mu _{3},...,\mu _{k} \right \}$

2：repeat

3：      令$C_{i}=\varnothing (1\leqslant i\leqslant k)$

4：     for j=1,2,…,m  do

5：           计算样本$x_{j}$与各均值向量$\mu_{i}(1\leqslant i\leqslant k)$的距离：$d_{ji}=\left \| x_{j}-\mu_{i} \right \|_{2}$;

6：            根据距离最近的均值向量确定$x_{j}$的簇标记：$\lambda _{j}=arg min_{i\in \left \{ 1,2,3,...,k \right \}}d_{ji}$；

7：            将样本$x_{j}$划入相应的簇：$C_{\lambda_{j}}=C_{\lambda_{j}}\cup \left \{ x_{j}  \right \}；$

8：        end for

 9：        for i=1,2,…,k do

10:               计算新均值向量：$\mu_{i}^{'}=\frac{1}{\left | C_{i}  \right |}\sum _{x\in C_{i}}x$;

11：             if  $\mu_{i}^{'}\neq \mu_{i}$  then

12：                 将当前均值向量 $\mu_{i}$更新为$\mu_{i}^{'}$

13：             else

14：                 保持当前均值不变

15：                 end if

16：                 end for

17：until  当前均值向量均未更新

输出：簇划分$C=\left \{ C_{1} ,C_{2},...,C_{k} \right \}$



以上算法流程引自周志华《机器学习》，从流程来看K-means算法计算步骤基本上可以概括为两个部分：（1）计算每一个对象到类簇中心的距离；（2）根据类簇内的对象计算新的簇类中心。

3. 编程实现
为了方便应用我们将其编写为一个M函数KMeans()，首先需要确定函数的输入输出。这里输入参数为：data，K，iniCentriods，iterations（其中data为输入的不带标号的数据集数据矩阵，大小为numOfDatanumOfAttributes，K为数据分的类簇数目，iniCentriods为自行指定的初始聚类中心矩阵，大小为KnumOfAttributes，iterations为算法迭代次数。）

输出参数为：Idx，centroids，Distance（Idx为返回的分类标号， centroids为每一类的中心，Distance为类内总距离）
根据前面2.2节中的算法流程编写Kmeans算法的MATLAB程序如下
%% Kmeans算法
% 输入：
% data 输入的不带分类标号的数据
% K 数据一共分多少类
% iniCentriods 自行指定初始聚类中心
% iterations 迭代次数

% 输出：
% Idx 返回的分类标号
% centroids 每一类的中心
% Distance 类内总距离

 
function [Idx,centroids,Distance]=KMeans(data,K,iniCentriods,iterations)
[numOfData,numOfAttr]=size(data); % numOfData是数据个数，numOfAttr是数据维数
centroids=iniCentriods;
%% 迭代
for iter=1:iterations
    pre_centroids=centroids;% 上一次求得的中心位置
    
    tags=zeros(numOfData,K);
    %% 寻找最近中心，更新中心
    for i=1:numOfData
        D=zeros(1,K);% 每个数据点与每个聚类中心的标准差
        Dist=D;
        
        % 计算每个点到每个中心点的标准差
        for j=1:K
            Dist(j)=norm(data(i,:)-centroids(j,:),2);
        end
        
        [minDistance,index]=min(Dist);% 寻找距离最小的类别索引
        tags(i,index)=1;% 标记最小距离所处的位置（类别）
    end
    
    
    %% 取均值更新聚类中心点
    for i=1:K
        if sum(tags(:,i))~=0
            % 未出现空类，计算均值作为下一聚类中心
            for j=1:numOfAttr
                centroids(i,j)=sum(tags(:,i).*data(:,j))/sum(tags(:,i));
            end
        else % 如果出现空类，从数据集中随机选中一个点作为中心
            randidx = randperm(size(data, 1));
            centroids(i,:) = data(randidx(1),:);
            tags(randidx,:)=0;
            tags(randidx,i)=1;
        end
    end
    
   
    if sum(norm(pre_centroids-centroids,2))<0.001  % 不断迭代直到位置不再变化
        break;
    end
    
    
end

%% 计算输出结果
Distance=zeros(numOfData,1);
Idx=zeros(numOfData,1);
for i=1:numOfData
    D=zeros(1,K);% 每个数据点与每个聚类中心的标准差
    Dist=D;
    % 计算每个点到每个中心点的标准差
    for j=1:K
        Dist(j)=norm(data(i,:)-centroids(j,:),2);
    end
    
    [distance,idx]=min(Dist);% 寻找距离最小的类别索引
    distance=Dist(idx);
    
    Distance(i)=distance;
    Idx(i)=idx;
end
Distance=sum(Distance,1);% 计算类内总距离
end

在以上代码中其最主要部分在于第18至58行，进行寻找最近中心和求取均值更新聚类中心点。值得注意的是，在聚类过程中可能会出现空类即代码第44行那样，为保证算法的继续运行，从数据集中随机选取一个点作为中心。
4. 聚类结果评价
为了验证编写的Kmeans函数的性能，这里对想用的UCI数据集Iris数据集进行聚类并计算聚类的准确率，Iris数据集可以在http://archive.ics.uci.edu/ml/index.php上下载得到。首先读取Iris数据集，自行指定初始聚类中心调用前面编写的KMeans函数进行聚类，然后计算聚类的准确率，其代码如下
clear 
data=load('Iris.txt');
data=data(:,2:end);

matrix=[5.9016,2.7484,4.3935,1.4339;6.8500,3.0737,5.7421,2.0711;5.0060,3.4280,1.4620,0.2460];
[Idx,C,distance]=KMeans(data,3,matrix,500);
Distance=sum(distance)

c1=Idx(1:50,1);c2=Idx(51:100,1);c3=Idx(101:150,1);
accuracy=(sum(c1==mode(Idx(1:50,1)))+sum(c2==mode(Idx(51:100,1)))+sum(c3==mode(Idx(101:150,1))))/150

为方便使用Iris数据集经过了一些整理，这里将最后一列的带字符串的标签Iris-setosa，Iris-versicolor，Iris-virginica分别用数字1,2,3代替并移到了第一列，所以第三行选取的是从第二列至最后一列的数据。第5行中的matrix是查阅论文得到的一个初始聚类中心，正好用来比对聚类结果。第6行则调用KMeans()函数进行聚类，得到聚类标号和类内距离。对每类的类内距离求和即得到总的距离Distance，如第7行。准确率的计算有点麻烦，因为不能直接用KMeans计算后得到的标号跟原数据集中的标号对比计算准确率，KMeans只需要也只能将那些“相似”的数据点聚集到一类中，而给这一类数据的标号却是可能跟原数据集不同的。
这里采用一个简单的方法，从原数据集的标签可以看出第1-50个数据点为一类（Iris-setosa），第51-100为一类（Iris-versicolor），第101-150为一类（Iris-virginica），因此只需确定每50个数据点中的聚类标号是不是一致。取它们之中数目最多的标号作为正确的个数，最终比上数据集的总数即为准确率。以上代码运行结果如下所示。
 
5. 类簇中心点的选取
KMeans算法本身思想比较简单，但是合理的确定K值和K个初始类簇中心点对于聚类效果的好坏有很大的影响。最简单的确定初始类簇中心点的方法是随机产生数据大小范围内的K个点作为初始的簇类中心点。随机产生初始点并进行测试的程序代码如下
clear
data=load('Iris.txt');
data=data(:,2:end);
K=3;


%% 产生随机初始点
[numOfData,numOfAttr]=size(data);   % numOfData是数据个数，numOfAttr是数据维数

centroids=zeros(K,numOfAttr);       % 随机初始化，最终迭代到每一类的中心位置
maxAttr=zeros(numOfAttr);        % 每一维最大的数
minAttr=zeros(numOfAttr);        % 每一维最小的数
for i=1:numOfAttr
    maxAttr(i)=max(data(:,i));    % 每一维最大的数
    minAttr(i)=min(data(:,i));    % 每一维最小的数
    for j=1:K
        centroids(j,i)=maxAttr(i)+(minAttr(i)-maxAttr(i))*rand();  % 随机初始化，选取每一维[min max]中初始化
    end
end

[Idx,C,distance]=KMeans(data,K,centroids,500);% 调用KMeans
Distance=sum(distance)% 计算类内距离之和

%% 计算准确率
c1=Idx(1:50,1);c2=Idx(51:100,1);c3=Idx(101:150,1);
Accuracy=(sum(c1==mode(Idx(1:50,1)))+sum(c2==mode(Idx(51:100,1)))+sum(c3==mode(Idx(101:150,1))))/numOfData

可以多运行几次以上代码，可以看出由于初始点事随机选取的每次运行得到的结果有所差异。这也是基本Kmeans算法的一个缺点，随着众多改进算法的提出Kmeans算法的这一问题也得到改善，深入了解的朋友可以查阅相关论文。
6. 结束语
本博文的完整MATLAB程序文件与整理好的数据集文件已经上传，下载即可运行。下载地址如下
点击下载：本文Kmeans算法M函数及测试完整文件
【公众号获取】

    本人微信公众号已创建，扫描以下二维码并关注公众号“AI技术研究与分享”，后台回复“KM20180516”即可获取全部资源文件信息。

由于编者能力有限，代码即使经过了多次校对，也难免会有疏漏之处。希望您能热心指出其中的错误，以便下次修改时能以一个更完美更严谨的样子，呈现在大家面前。同时如果有更好的实现方法也请您不吝赐教。

在学习聚类算法得时候并没有涉及到评估指标，主要原因是聚类算法属于非监督学习，并不像分类算法那样可以使用训练集或测试集中得数据计算准确率、召回率等。那么如何评估聚类算法得好坏呢？好的聚类算法,一般要求类簇具有：
高的类内 (intra-cluster) 相似度
低的类间 (inter-cluster) 相似度
对于聚类算法大致可分为 2类度量标准：
内部评估的方法：通过一个单一的量化得分来评估算法好坏；该类型的方法
外部评估的方法：通过将聚类结果与已经有“ground truth”分类进行对比。要么通过人类进行手动评估，要么通过一些指标在特定的应用场景中进行聚类用法的评估。不过该方法是有问题的，如果真的有了label，那么还需要聚类干嘛，而且实际应用中，往往都没label；另一方面，这些label只反映了数据集的一个可能的划分方法，它并不能告诉你存在一个不同的更好的聚类算法。
内部评价指标
当一个聚类结果是基于数据聚类自身进行评估的，这一类叫做内部评估方法。如果某个聚类算法聚类的结果是类间相似性低，类内相似性高，那么内部评估方法会给予较高的分数评价。不过内部评价方法的缺点是：
那些高分的算法不一定可以适用于高效的信息检索应用场景；
这些评估方法对某些算法有倾向性，如k-means聚类都是基于点之间的距离进行优化的，而那些基于距离的内部评估方法就会过度的赞誉这些生成的聚类结果。
这些内部评估方法可以基于特定场景判定一个算法要优于另一个，不过这并不表示前一个算法得到的结果比后一个结果更有意义。这里的意义是假设这种结构事实上存在于数据集中的，如果一个数据集包含了完全不同的数据结构，或者采用的评价方法完全和算法不搭，比如k-means只能用于凸集数据集上，许多评估指标也是预先假设凸集数据集。在一个非凸数据集上不论是使用k-means还是使用假设凸集的评价方法，都是徒劳的。
SSE(和方差)
该统计参数计算的是拟合数据和原始数据对应点的误差的平方和，计算公式如下：
SSE越接近于0，说明模型选择和拟合更好，数据预测也越成功。
#断崖碎石图选取最优K值import pandas as pd  from sklearn.cluster import KMeans  import matplotlib.pyplot as plt  '利用SSE选择k'  SSE = []  # 存放每次结果的误差平方和  for k in range(1,9):      estimator = KMeans(n_clusters=k)  # 构造聚类器      estimator.fit(df[['calories','sodium','alcohol','cost']])      SSE.append(estimator.inertia_)  N = range(1,9)  plt.xlabel('k')  plt.ylabel('SSE')  plt.plot(N,SSE,'o-')  plt.show()
轮廓系数 Silhouette Coefficient
轮廓系数适用于实际类别信息未知的情况。对于单个样本，设a是与它同类别中其他样本的平均距离，b是与它距离最近不同类别中样本的平均距离，其轮廓系数为：
对于一个样本集合，它的轮廓系数是所有样本轮廓系数的平均值。轮廓系数的取值范围是[-1,1]，同类别样本距离越相近不同类别样本距离越远，分数越高。缺点：不适合基高密度的聚类算法DBSCAN。
from sklearn import metricsfrom sklearn.metrics import pairwise_distancesfrom sklearn import datasetsdataset = datasets.load_iris()X = dataset.datay = dataset.target import numpy as npfrom sklearn.cluster import KMeanskmeans_model = KMeans(n_clusters=3, random_state=1).fit(X)labels = kmeans_model.labels_metrics.silhouette_score(X, labels, metric='euclidean')
Calinski-Harabaz Index
在真实的分群label不知道的情况下，Calinski-Harabasz可以作为评估模型的一个指标。Calinski-Harabasz指标通过计算类中各点与类中心的距离平方和来度量类内的紧密度，通过计算各类中心点与数据集中心点距离平方和来度量数据集的分离度，CH指标由分离度与紧密度的比值得到。从而，CH越大代表着类自身越紧密，类与类之间越分散，即更优的聚类结果。
其中m为训练样本数，k是类别个数，是类别之间协方差矩阵，是类别内部数据协方差矩阵，为矩阵的迹。也就是说，类别内部数据的协方差越小越好，类别之间的协方差越大越好，这样的Calinski-Harabasz分数会高。同时，数值越小可以理解为：组间协方差很小，组与组之间界限不明显。
优点
当 cluster (簇)密集且分离较好时，分数更高，这与一个标准的 cluster(簇)有关。
得分计算很快与轮廓系数的对比，最大的优势：快！相差几百倍！毫秒级。
缺点
凸的簇的 Calinski-Harabaz index(Calinski-Harabaz 指数)通常高于其他类型的 cluster(簇)，例如通过 DBSCAN 获得的基于密度的 cluster(簇)。所以不适合基于密度的聚类算法，DBSCAN。
import numpy as npfrom sklearn.cluster import KMeanskmeans_model = KMeans(n_clusters=3, random_state=1).fit(X)labels = kmeans_model.labels_print(metrics.calinski_harabaz_score(X, labels))
Compactness(紧密性)(CP)
CP计算每一个类各点到聚类中心的平均距离CP越低意味着类内聚类距离越近。著名的 K-Means 聚类算法就是基于此思想提出的。缺点：没有考虑类间效果。
Separation(间隔性)(SP)
SP计算各聚类中心两两之间平均距离，SP越高意味类间聚类距离越远。缺点：没有考虑类内效果。
Davies-Bouldin Index(戴维森堡丁指数)(分类适确性指标)(DB)(DBI)
DB计算任意两类别的类内距离平均距离(CP)之和除以两聚类中心距离求最大值。DB越小意味着类内距离越小同时类间距离越大。该指标的计算公式：
其中n是类别个数，是第i个类别的中心，是类别i中所有的点到中心的平均距离； 
 中心点和之间的距离。算法生成的聚类结果越是朝着类内距离最小(类内相似性最大)和类间距离最大(类间相似性最小)变化，那么Davies-Bouldin指数就会越小。缺点：因使用欧式距离所以对于环状分布聚类评测很差。
from sklearn import datasets from sklearn.cluster import KMeans from sklearn.metrics import davies_bouldin_score from sklearn.datasets.samples_generator import make_blobs   # loading the dataset X, y_true = make_blobs(n_samples=300, centers=4,                         cluster_std=0.50, random_state=0)   # K-Means kmeans = KMeans(n_clusters=4, random_state=1).fit(X)   # we store the cluster labels labels = kmeans.labels_   print(davies_bouldin_score(X, labels))
Dunn Validity Index (邓恩指数)(DVI)
DVI计算任意两个簇元素的最短距离(类间)除以任意簇中的最大距离(类内)。DVI越大意味着类间距离越大同时类内距离越小。
其中 
 表示类别,之间的距离； 
 表示类别内部的类内距离：
类间距离 可以是任意的距离测度，例如两个类别的中心点的距离；
类内距离 可以以不同的方法去测量，例如类别kk中任意两点之间距离的最大值。
因为内部评估方法是搜寻类内相似最大，类间相似最小，所以算法生成的聚类结果的Dunn指数越高，那么该算法就越好。缺点：对离散点的聚类测评很高、对环状分布测评效果差。
import pandas as pd from sklearn import datasets from jqmcvi import base  # loading the dataset X = datasets.load_iris() df = pd.DataFrame(X.data)  # K-Means from sklearn import cluster k_means = cluster.KMeans(n_clusters=3) k_means.fit(df) #K-means training y_pred = k_means.predict(df)  # We store the K-means results in a dataframe pred = pd.DataFrame(y_pred) pred.columns = ['Type']  # we merge this dataframe with df prediction = pd.concat([df, pred], axis = 1)  # We store the clusters clus0 = prediction.loc[prediction.Species == 0] clus1 = prediction.loc[prediction.Species == 1] clus2 = prediction.loc[prediction.Species == 2] cluster_list = [clus0.values, clus1.values, clus2.values]  print(base.dunn(cluster_list))
外部评价指标
在外部评估方法中，聚类结果是通过使用没被用来做训练集的数据进行评估。例如已知样本点的类别信息和一些外部的基准。这些基准包含了一些预先分类好的数据，比如由人基于某些场景先生成一些带label的数据，因此这些基准可以看成是金标准。这些评估方法是为了测量聚类结果与提供的基准数据之间的相似性。然而这种方法也被质疑不适用真实数据。
纯度(Purity)
纯度(Purity)是一种简单而透明的评估手段，为了计算纯度(Purity)，我们把每个簇中最多的类作为这个簇所代表的类，然后计算正确分配的类的数量，然后除以N。形式化表达如下：
其中：
是聚类的集合，表示第k个聚类的集合。
是文档集合，表示第J个文档。
表示文档总数。
上述过程即给每个聚类簇分配一个类别,且这个类别的样本在该簇中出现的次数最多，然后计算所有 K 个聚类簇的这个次数之和再归一化即为最终值。Purity值在0～1之间 ,越接近1表示聚类结果越好。
如图认为x代表一类文档，o代表一类文档，方框代表一类文档。如上图的purity = ( 3+ 4 + 5) / 17 = 0.71，其中第一类正确的有5个，第二个4个，第三个3个，总文档数17。
当簇的数量很多的时候，容易达到较高的纯度——特别是，如果每个文档都被分到独立的一个簇中，那么计算得到的纯度就会是1。因此，不能简单用纯度来衡量聚类质量与聚类数量之间的关系。另外Purity无法用于权衡聚类质量与簇个数之间的关系。
def purity(result, label):    # 计算纯度     total_num = len(label)    cluster_counter = collections.Counter(result)    original_counter = collections.Counter(label)     t = []    for k in cluster_counter:        p_k = []        for j in original_counter:            count = 0            for i in range(len(result)):                if result[i] == k and label[i] == j: # 求交集                    count += 1            p_k.append(count)        temp_t = max(p_k)        t.append(temp_t)        return sum(t)/total_num
标准化互信息(NMI)
互信息(Normalized Mutual Information)是用来衡量两个数据分布的吻合程度。也是一有用的信息度量，它是指两个事件集合之间的相关性。互信息越大，词条和类别的相关程度也越大。NMI (Normalized Mutual Information) 即归一化互信息：
其中,表示互信息(Mutual Information),为熵，当 log 取 2 为底时，单位为 bit，取 e 为底时单位为 nat。
其中,
可以分别看作样本 (document) 属于聚类簇, 属于类别, 同时属于的概率。第二个等价式子则是由概率的极大似然估计推导而来。
互信息 
 表示给定类簇信息的前提条件下,类别信息的增加量，或者说其不确定度的减少量。直观地，互信息还可以写出如下形式：
互信息的最小值为 0, 当类簇相对于类别只是随机的, 也就是说两者独立的情况下,对于未带来任何有用的信息.如果得到的与关系越密切, 那么 
 值越大。如果完整重现了, 此时互信息最大：
当K=N时,即类簇数和样本个数相等，MI 也能达到最大值。所以 MI 也存在和纯度类似的问题,即它并不对簇数目较大的聚类结果进行惩罚,因此也不能在其他条件一样的情况下,对簇数目越小越好的这种期望进行形式化。NMI 则可以解决上述问题,因为熵会随着簇的数目的增长而增大。当K=N时, 
 会达到其最大值 
 , 此时就能保证 NMI 的值较低。之所以采用
作为分母是因为它是 
 的紧上界, 因此可以保证
。
示例：
gnd 是 ground truth 的意思，grps 表示聚类后的 groups. 问题：计算序列 gnd 和 grps 的 NMI.
先计算联合概率分布
计算边际分布
计算熵和互信息
计算 NMI
代码实现：
def NMI(result, label):    # 标准化互信息     total_num = len(label)    cluster_counter = collections.Counter(result)    original_counter = collections.Counter(label)        # 计算互信息量    MI = 0    eps = 1.4e-45 # 取一个很小的值来避免log 0        for k in cluster_counter:        for j in original_counter:            count = 0            for i in range(len(result)):                if result[i] == k and label[i] == j:                    count += 1            p_k = 1.0*cluster_counter[k] / total_num            p_j = 1.0*original_counter[j] / total_num            p_kj = 1.0*count / total_num            MI += p_kj * math.log(p_kj /(p_k * p_j) + eps, 2)        # 标准化互信息量    H_k = 0    for k in cluster_counter:        H_k -= (1.0*cluster_counter[k] / total_num) * math.log(1.0*cluster_counter[k] / total_num+eps, 2)    H_j = 0    for j in original_counter:        H_j -= (1.0*original_counter[j] / total_num) * math.log(1.0*original_counter[j] / total_num+eps, 2)            return 2.0 * MI / (H_k + H_j)
sklearn中自带的方法：
from sklearn.metrics.cluster import normalized_mutual_info_scoreprint(normalized_mutual_info_score([0, 0, 1, 1], [0, 0, 1, 1]))
调整互信息AMI( Adjusted mutual information)
已知聚类标签与真实标签，互信息(mutual information)能够测度两种标签排列之间的相关性，同时忽略标签中的排列。有两种不同版本的互信息以供选择，一种是Normalized Mutual Information(NMI),一种是Adjusted Mutual Information(AMI)。
假设U与V是对N个样本标签的分配情况，则两种分布的熵(熵表示的是不确定程度)分别为：
其中：
是从U中随机选取的对象到类的概率
从V中随机选取的对象到类的概率
U与V之间的互信息(MI)定义为：
其中
是随机选择的对象落入两个类的概率和。
调整互信息(Adjusted mutual information)定义为：
MI的期望可以用以下公式来计算。在这个方程式中， 
 为元素的数量， 
 为元素的数量：
利用基于互信息的方法来衡量聚类效果需要实际类别信息，MI与NMI取值范围为[0,1]，AMI取值范围为[-1,1]，它们都是值越大意味着聚类结果与真实情况越吻合。
优点
随机(统一)标签分配的AMI评分接近0)
有界范围 [0, 1]: 接近 0 的值表示两个主要独立的标签分配，而接近 1 的值表示重要的一致性。此外，正好 0 的值表示 purely(纯粹) 独立标签分配，正好为 1 的 AMI 表示两个标签分配相等(有或者没有 permutation)。
对簇的结构没有作出任何假设: 可以用于比较聚类算法
缺点：
与 inertia 相反， MI-based measures 需要了解 ground truth classes，而在实践中几乎不可用，或者需要人工标注或手动分配(如在监督学习环境中)。然而，基于 MI-based measures (基于 MI 的测量方式)也可用于纯无人监控的设置，作为可用于聚类模型选择的 Consensus Index (共识索引)的构建块。
NMI 和 MI 没有调整机会。
from sklearn import metricslabels_true = [0, 0, 0, 1, 1, 1]labels_pred = [0, 0, 1, 1, 2, 2] print(metrics.adjusted_mutual_info_score(labels_true, labels_pred))
Rand index兰德指数
兰德指数 (Rand index, RI), 将聚类看成是一系列的决策过程,即对文档集上所有N(N-1)/2个文档 (documents) 对进行决策。当且仅当两篇文档相似时,我们将它们归入同一簇中。
Positive:
TP 将两篇相似文档归入一个簇 (同 – 同)
TN 将两篇不相似的文档归入不同的簇 (不同 – 不同)
Negative:
FP 将两篇不相似的文档归入同一簇 (不同 – 同)
FN 将两篇相似的文档归入不同簇 (同- 不同) (worse)
RI 则是计算「正确决策」的比率(精确率, accuracy)：
RI取值范围为[0,1]，值越大意味着聚类结果与真实情况越吻合。
def contingency_table(result, label):        total_num = len(label)        TP = TN = FP = FN = 0    for i in range(total_num):        for j in range(i + 1, total_num):            if label[i] == label[j] and result[i] == result[j]:                TP += 1            elif label[i] != label[j] and result[i] != result[j]:                TN += 1            elif label[i] != label[j] and result[i] == result[j]:                FP += 1            elif label[i] == label[j] and result[i] != result[j]:                FN += 1    return (TP, TN, FP, FN) def rand_index(result, label):    TP, TN, FP, FN = contingency_table(result, label)    return 1.0*(TP + TN)/(TP + FP + FN + TN)
调整兰德系数 (Adjusted Rand index)
对于随机结果，RI并不能保证分数接近零。为了实现“在聚类结果随机产生的情况下，指标应该接近零”，调整兰德系数(Adjusted rand index)被提出，它具有更高的区分度：
ARI取值范围为[-1,1]，值越大意味着聚类结果与真实情况越吻合。从广义的角度来讲，ARI衡量的是两个数据分布的吻合程度。
优点：
对任意数量的聚类中心和样本数，随机聚类的ARI都非常接近于0
取值在［－1，1］之间，负数代表结果不好，越接近于1越好
对簇的结构不需作出任何假设：可以用于比较聚类算法。
缺点：
与 inertia 相反，ARI 需要 ground truth classes 的相关知识，ARI需要真实标签，而在实践中几乎不可用，或者需要人工标注者手动分配(如在监督学习环境中)。然而，ARI 还可以在 purely unsupervised setting (纯粹无监督的设置中)作为可用于 聚类模型选择(TODO)的共识索引的构建块。
from sklearn import metricslabels_true = [0, 0, 0, 1, 1, 1]labels_pred = [0, 0, 1, 1, 2, 2] print(metrics.adjusted_rand_score(labels_true, labels_pred)) 
F值方法
这是基于上述RI方法衍生出的一个方法，我们可以 FN 罚更多,通过取 
 中的大于 1, 此时实际上也相当于赋予召回率更大的权重：
RI方法有个特点就是把准确率和召回率看得同等重要，事实上有时候我们可能需要某一特性更多一点，这时候就适合F值方法。
def precision(result, label):    TP, TN, FP, FN = contingency_table(result, label)    return 1.0*TP/(TP + FP) def recall(result, label):    TP, TN, FP, FN = contingency_table(result, label)    return 1.0*TP/(TP + FN) def F_measure(result, label, beta=1):    prec = precision(result, label)    r = recall(result, label)    return (beta*beta + 1) * prec * r/(beta*beta * prec + r)
Fowlkes-Mallows scores
Fowlkes-Mallows Scores(FMI) FMI是成对的precision(精度)和recall(召回)的几何平均数。取值范围为 [0,1]，越接近1越好。定义为：
代码实现：
from sklearn import metricslabels_true = [0, 0, 0, 1, 1, 1]labels_pred = [0, 0, 1, 1, 2, 2] print(metrics.fowlkes_mallows_score(labels_true, labels_pred))
调和平均V-measure
说V-measure之前要先介绍两个指标：
同质性(homogeneity)：每个群集只包含单个类的成员。
完整性(completeness)：给定类的所有成员都分配给同一个群集。
同质性和完整性分数基于以下公式得出：
其中 
 是给定给定簇赋值的类的条件熵，由以下公式求得：
 是类熵，公式为：
其中，n是样本总数，和分别属于类c和类k的样本数，而是从类c划分到类k的样本数量。条件熵H(K|C)和类熵H(K)，根据以上公式对称求得。
V-measure是同质性homogeneity和完整性completeness的调和平均数，V-measure取值范围为 [0,1]，越大越好，但当样本量较小或聚类数据较多的情况，推荐使用AMI和ARI。公式：
代码实现：
from sklearn import metricslabels_true = [0, 0, 0, 1, 1, 1]labels_pred = [0, 0, 1, 1, 2, 2] print(metrics.homogeneity_score(labels_true, labels_pred))print(metrics.completeness_score(labels_true, labels_pred))print(metrics.v_measure_score(labels_true, labels_pred))
优点：
分数明确：从0到1反应出最差到最优的表现；
解释直观：差的调和平均数可以在同质性和完整性方面做定性的分析；
对簇结构不作假设：可以比较两种聚类算法如k均值算法和谱聚类算法的结果。
缺点：
以前引入的度量在随机标记方面没有规范化，这意味着，根据样本数，集群和先验知识，完全随机标签并不总是产生相同的完整性和均匀性的值，所得调和平均值V-measure也不相同。特别是，随机标记不会产生零分，特别是当簇的数量很大时。
当样本数大于一千，聚类数小于10时，可以安全地忽略该问题。对于较小的样本量或更大数量的集群，使用经过调整的指数(如调整兰德指数)更为安全。
这些指标要求的先验知识，在实践中几乎不可用或需要手动分配的人作注解者(如在监督学习环境中)。
Jaccard 指数
该指数用于量化两个数据集之间的相似性，该值得范围为0-1.其中越大表明两个数据集越相似：
该指数和近年来的IOU计算方法一致
Dice 指数
该指数是基于jaccard指数上将TP的权重置为2倍。

上一篇博客中简单介绍了K均值聚类算法，在本篇博客中介绍一下关于谱聚类算法，简单谈一谈自己的心得。简单介绍一下谱聚类算法


谱聚类算法建立在谱图理论基础上，与传统的聚类算法相比，它具有能在任意形状的样本空间上聚类且收敛于全局最优解的优点。该算法首先根据给定的样本数据集定义一个描述成对数据点相似度的亲合矩阵,并且计算矩阵的特征值和特征向量 ， 然后选择合适 的特征向量聚类不同的数据点。谱聚类算法最初用于计算机视觉 、VLS I 设计等领域， 最近才开始用于机器学习中，并迅速成为国际上机器学习领域的研究热点。谱聚类算法建立在图论中的谱图理论基础上，其本质是将聚类问题转化为图的最优划分问题，是一种点对聚类算法，对数据聚类具有很好的应用前景。


详细了解谱聚类算法的原理需要一定的图论学和矩阵论的知识，详细的原理及介绍在这篇博客中有介绍。

核心思想

将所有数据看成一个个空间中的点，每个点之间有边连接，每条边有不同的权重（距离远的点之间的边权重较小，距离近的点之间的边权重较低），通过选择合适的分割，使得每个子图内的边权重之和尽可能大，且子图与子图间的边权重之和尽可能小，基于此标准将数据进行分类

具体介绍

举例而言，张一，张二，张三，王一，王二，王三，李一，李二，李三是同一学校的学生，张一、二、三相互认识，王一、二、三相互认知，李一、二、三相互认识，张一、王一、李一相互认识，其余人互不认识，理论上，通过关系网，他们之间都有直接或间接的关系，如下图



切断王一、张一、李一的联系，就将其分为了三类，谱聚类的原理大致如此，只是没有那么简单而已

图的定义

一般来说，我们用点与边来定义一个图，记点的集合为V，边的集合为E，则我们称G（V,E）是一个图，V就是我们数据中的一行行数据（没有标签等），即



共有n条数据，每条数据有m个属性。紧接着我们引入权重与度的定义，记是点到点边上的权重，这里我们是无向图，所以有



紧接着我们引入度的概念，点的度是所有与之相连的边的权重之和，即



度矩阵与邻接矩阵

由上面的定义我们可以得到一个n*n的度矩阵D和邻接矩阵W





相似矩阵

以上我们只给出了权重的定义，但是如何定量的描述权重，是谱聚类算法中的一个核心，定量刻画权重的方法有很多，比如ϵ-邻近法，K邻近法和全连接法，最常用的还是全连接法，全连接法有可以细分成很多种，这里只阐述最常用的一种——高斯核函数RBF，所以对权重有如下定量描述



其中与是数据点，是需要给定的参数（不同的值为对最后的聚类效果有影响）

用此种方法得到的相似矩阵与邻接矩阵一样。

值得注意的是：引入相似矩阵（样本点距离度量的相似矩阵 ）的目的是为了得到邻接矩阵，不同的方法得到的邻接矩阵是不一样的，所以一般是需要给定相似矩阵的

拉普拉斯矩阵

基于上述的度矩阵和邻接矩阵我们可以得到一个新的矩阵——拉普拉斯矩阵（目的是将数据进行降维处理）



由于D与W均是对称矩阵，所以拉普拉斯矩阵有一些很好的性质

L是一个对称矩阵
	
L的所有特征值均是实数
	
L是一个半正定矩阵，所以特征值均大于0，且最小特征值为
具体的证明不在赘述，有兴趣的可以看看

切图聚类

切图是谱聚类算法的核心，切图的方法有很多，但最常用的还是Ncut切图，因为原理比较复杂，这里就不在说明了，感兴趣的可以自行查找，切图的目的是为了分类，所以需要进行以下操作。

1、对拉普拉斯矩阵进行标准化



值得注意的是：是指对角线上的元素进行运算（因为其余均为0元素）

2、计算标准化后的拉普拉斯矩阵的特征值,对数据进行降维处理

选取最小的k个特征值对于的特征向量，将其放入同一个矩阵中，组成一个n*k维的矩阵F(因为每一个特征向量都是一个n维的列向量)

3、利用K均值聚类对矩阵F进行聚类处理

矩阵F的每一行对应一条数据，每条数据有k个属性（从n个属性缩减到k个属性，这就是一种降维处理），然后可以利用K均值聚类算法进行聚类处理。（这里降维后的维度可以不一定是k，也可以是k1，即可以任意指定一个维度，并不需要与聚类后的维度一致，但是为了方便起见，我们一般都让其保持一致）

到这里基本上谱聚类算法的主要内容已经阐述完毕，接下来总结一下算法思路

算法流程

输入相似矩阵S,高斯核函数的参数，以及所分类数k
	
计算度矩阵D和邻接矩阵W
	
计算拉普拉斯矩阵L
	
标准化矩阵L，求最小k个特征值对应的特征向量组成的矩阵
	
基于K均值聚类进行聚类
代码实现

"""
谱聚类算法
"""

# 读取数据
from numpy.linalg import *
import numpy as np
import pandas as pd
from sklearn.cluster import SpectralClustering, KMeans
import math
data = pd.read_excel(
    "D:\\wangyang\\聚类算法\\算法程序\\谱聚类算法\\data.xlsx", names=range(7), header=None)  # 添加列索引
n, m = data.shape  # 获取数据的行数与列数
sigma = 0.5
k = 3
# 计算邻接矩阵


def Adjacency_matrix(data):
    """
    计算邻接矩阵
    data：数据集，n*m的矩阵
    不足：两个嵌套循环导致的计算量还是很大的，所以还需要优化
    """
    W = np.zeros((n, n))  # 初始化邻接矩阵
    for i in range(n):
        x_i = np.array(data.loc[i])
        for j in range(n):
            x_j = np.array(data.loc[j])
            w_ij = np.linalg.norm(x_i - x_j, ord=2)  # 计算向量的2-范数
            W[i][j] = math.exp(-w_ij**2 / (2 * sigma**2))
    return W


def Degree_matrix(W):
    """
    计算度矩阵
    输入：W是邻接矩阵，n*n
    """
    D = np.zeros((n, n))
    for i in range(n):
        D[i][i] = sum(W[i])
    return D


def diagpower_matrix(X, k):
    """
    对矩阵对角线上的元素进行运算
    输入：X是一个方阵，类型是dataframe
    k是幂次方，一个实数
    """
    nn, mm = X.shape
    for i in range(nn):
        X[i][i] = X[i][i]**k
    return X


W = Adjacency_matrix(data)
D = Degree_matrix(W)
# 计算拉普拉斯矩阵
L = D - W
L = diagpower_matrix(D, -0.5) * L * diagpower_matrix(D, -0.5)  # 矩阵标准化
eigenvalues, Feature_vector = np.linalg.eig(L)
# k个最小特征值对应的特征向量组成的n*k的矩阵
data2 = Feature_vector[np.argsort(eigenvalues)[:k]].T
result = KMeans(n_clusters=k, max_iter=100)  # 构造聚类器
result.fit(data2)  # 聚类
label = result.labels_  # 获取聚类标签
centroids = result.cluster_centers_  # 获取聚类中心
# 写入文件
label = pd.DataFrame(label)
label.to_excel("D:\\wangyang\\聚类算法\\算法程序\\谱聚类算法\\label.xlsx",
               index=False, header=False)
centroids = pd.DataFrame(centroids)
centroids.to_excel("D:\\wangyang\\聚类算法\\算法程序\\谱聚类算法\\centroids.xlsx",
                   index=False, header=False)


相关数据集请点击下载，直接下载运行即可，由于只是为了实现上述功能以及作者水平有限，代码的运行效率不是很高，还有一些参数选取的不合适，所以最后的分类效果不是特别理想。不过本博客的初衷只是为了介绍谱聚类算法的原理并复现，对精度要求不是很在意，其实如果需要python中自带的sklearn库中就已经自带了谱聚类算法，感兴趣的可以点击看看

层次聚类（Hierarchical Clustering算法）

层次聚类算法又称为树聚类算法，它根据数据之间的距离，透过一种层次架构方式，反复将数据进行聚合，创建一个层次以分解给定的数据集。 

常用于一维数据的自动分组

层次聚类方法 

hclust(dist)

dist 样本的距离矩阵
距离矩阵的计算方式 

dist(data)

data 样本数据
层次聚类的代码实现：

pColumns <- c('Sepal.Length', 'Sepal.Width', 'Petal.Length', 'Petal.Width');

data <- iris[ , pColumns]

data.e <- dist(data)

data.m <- as.matrix(data.e)

model <- hclust(data.e)

plot(model)

result <- cutree(model, k=3)

table(iris[, 5], result)
 result
              1  2  3
  setosa     50  0  0
  versicolor  0 23 27
  virginica   0 49  1

我们可以看到，层次聚类对这份数据的聚类得到的结果并不是太好！