西瓜书习题试答-第09章-聚类

试答系列：“西瓜书”-周志华《机器学习》习题试答
系列目录
[第01章：绪论]
[第02章：模型评估与选择]
[第03章：线性模型]
[第04章：决策树]
[第05章：神经网络]
[第06章：支持向量机]
第07章：贝叶斯分类器
第08章：集成学习
第09章：聚类
第10章：降维与度量学习
第11章：特征选择与稀疏学习
第12章：计算学习理论(暂缺)
第13章：半监督学习
第14章：概率图模型
(后续章节更新中...)

---

目录

• 9.1 试证明：p≥1时，闵可夫斯基距离满足距离度量的四条基本性质；0≤p<1时，闵可夫斯基距离不满足直递性，但满足非负性、同一性、对称性；p趋向无穷大时，闵可夫斯基距离等于对应分量的最大绝对距离，即
• 9.2 同一样本空间中的集合X与Z之间的距离可通过“豪斯多夫距离”(Hausdorff distance)计算：
• 9.3 试析k均值算法能否找到最小化式(9.24)的最优解。
• 9.4 试编程实现k均值算法，设置三组不同的k值、三组不同初始中心点，在西瓜数据集4.0上进行实验比较，并讨论什么样的初始中心有利于取得好结果。
• 9.5 基于DBSCAN的概念定义，若\(x\)为核心对象，由\(x\)密度可达的所有样本构成的集合为\(X\)。试证明：\(X\)满足连接性(3.39)与最大性(9.40)。
• 9.6 试析AGNES算法使用最小距离和最大距离的区别。
• 9.7 聚类结果中若每个簇都有一个凸包(包含簇样本的凸多面体)，且这些凸包不相交，则称为凸聚类。试析本章介绍的哪些聚类算法只能产生凸聚类，哪些能产生非凸聚类。
• 9.8 试设计一个聚类性能度量指标，并与9.2节中的指标比较。（暂缺）
• 9.9 试设计一个能用于混合属性的非度量距离。（暂缺）
• 9.10 试设计一个能自动确定聚类数的改进k均值算法，编程实现并在西瓜数据集4.0上运行。
• 附：编程代码
  • 习题9.4 (Python)
  • 习题9.10 (Python)

---

9.1 试证明：p≥1时，闵可夫斯基距离满足距离度量的四条基本性质；0≤p<1时，闵可夫斯基距离不满足直递性，但满足非负性、同一性、对称性；p趋向无穷大时，闵可夫斯基距离等于对应分量的最大绝对距离，即

\[\lim_{p\rightarrow \infty}(\sum_{u=1}^n|x_{iu}-x_{ju}|^p)^{\frac{1}{p}}=\max_u |x_{iu}-x_{ju}| \]

证明：

• 非负性、同一性、对称性
  对于所有的p≥0,由于距离公式中的绝对值形式，决定了最终闵可夫斯基距离的\(dist_p(x_i,x_j)\)的非负性、同一性、对称性：\(|x_{iu}-x_{ju}|≥0\)，当且仅当\(x_{iu}=x_{ju}\)时等于0，而且\(|x_{iu}-x_{ju}|=|x_{ju}-x_{iu}|\)。

• 直递性
  基于闵可夫斯基（Minkowski）不等式，可以证得p≥1时满足直递性，0≤p≤1时不满足。

  定理：闵可夫斯基（Minkowski）不等式
  对任意p≥1以及\(x,y\in R^n\),有

  \[|x+y|_p≤|x|_p+|y|_p \]

  其中,\(|x|_p=(\sum|x_i|_p)^{1/p}=distm_k(x,0)\)。
  如果，\(x_1,…,x_n,y_1,…,y_n>0，p<1\)，则"≤"可以变为"≥"。

  令\(x=x_i-x_k，y=x_k-x_j\)，即得：p≥1时的直递性满足，p<1时，直递性不满足。
  p≥1时的闵可夫斯基（Minkowski）不等式的证明过程如下：

  首先需要证明一个引理和Holder不等式，然后再证明闵可夫斯基不等式。

  引理 \(a^λb^{1-λ}≤λa+(1-λ)b\)，其中a,b≥0, 0≤λ≤1。
  证明思路：两边取对数，由于对数ln函数是凸函数，利用Jensen不等式可证得。

  定理(Holder不等式)
  对任意的\(1\leq p,q\leq \infty,\frac{1}{p}+\frac{1}{q}=1\)，以及\(x,y\in R^n\)，有

  \[\sum_{i=1}^n|x_i y_i|\leq|x|_p |x|_q \]

  当p=q=2时，Holder不等式退化为柯西-施瓦茨不等式(文字可表述为：两个向量的内积不大于两个模长之积)。
  证明：由上面的引理有

  \[\frac{|x_i||y_i|}{|x|_p|y|_p}=\left(\frac{|x_i|^p}{|x|_p^p}\right)^{1/p}\cdot\left(\frac{|y_i|^q}{|y|_q^q}\right)^{1/q}\leq\frac{1}{p}\frac{|x_i|^p}{|x|_p^p}+\frac{1}{q}\frac{|y_i|^q}{|y|_q^q} \]

  不等式两边对i求和便有：

  \[\frac{1}{|x|_p|y|_p}\sum_{i=1}^n |x_i y_i|\leq \frac{1}{p}+\frac{1}{q}=1 \]

  定理得证。

  定理(Minkowski不等式)
  内容见前文
  证明：只需考虑1<p<∞的情况，p=1或者∞的情形易证。当1<p<∞时有

  \[\begin{aligned} |x+y|_p^p&=\sum_i|x_i+y_i|^p\\ &=\sum_i |x_i+y_i|^{p-1}|x_i+y_i|\\ &\leq\sum_i |x_i+y_i|^{p-1}(|x_i|+|y_i|) \end{aligned}\]

  由Holder不等式

  \[\begin{aligned} \sum_i |x_i+y_i|^{p-1}|x_i|&\leq|(x+y)^{p-1}|_q|x|_p\\ &=\left(\sum_i |x_i+y_i|^{(p-1)q}\right)^{1/q}|x|_p\\ &=\left(\sum_i |x_i+y_i|^p\right)^{1/q}|x|_p\\ &=|x+y|_p^{p/q}|x|_p \end{aligned}\]

  第三行的利用了关系(p-1)q=p，同理有：

  \[\sum_i |x_i+y_i|^{p-1}|y_i|\leq|x+y|_p^{p/q}|y|_p \]

  结合以上三个不等式有：

  \[|x+y|_p^p\leq|x+y|_p^{p/q}(|x|_p+|y|_p) \]

  不等式两边同除\(|x+y|_p^{p/q}\)即得闵可夫斯基不等式。

  以上证明过程参考了这篇博文，对于p<1时的情况，暂未找到证明方法。

• \(p\rightarrow\infty\)时的闵可夫斯基距离
  设第u个特征差(坐标差)：\(|x_{iu}-x_{ju}|=Δu\)，最大特征差为\(\max_uΔ_u=ΔM\)，则有

  \[\begin{aligned} \lim_{p\rightarrow\infty}dist_{mk}(x_i,x_j)&=\lim_{p\rightarrow\infty}[\sum_u (\Delta_u)^p]^{1/p}\\ &=\Delta M\cdot \lim_{p\rightarrow\infty}[\sum_u (\frac{\Delta_u}{\Delta M})^p]^{1/p}\\ &=\Delta M\cdot [\sum_u Ⅱ(\Delta_u = \Delta M)]^0\\ &=\Delta M\\ &=\max_u |x_{iu}-x_{ju}| \end{aligned}\]

  其中\(∑_uⅡ(\Delta_u=\Delta_M)\)代表特征差等于最大特征差的个数，其结果至少有1个。

9.2 同一样本空间中的集合X与Z之间的距离可通过“豪斯多夫距离”(Hausdorff distance)计算：

\[dist_H (X,Z)=max(dist_h (X,Z),dist_h(Z,X))\tag{9.44} \]

其中

\[dist_h(X,Z)=\max_{x\in X}\min_{z\in Z}||x-z||_2\tag{9.45} \]

试证明：豪斯多夫距离满足距离度量的四条基本性质。
证明：首先，我们来直观的理解一下豪斯多夫距离。

下面表示两个点的欧式距离时，为了简便，右下角的角标2略去不写，比如将\(|x-y|_2\)表示为\(|x-y|\)。
如上图所示，关于\(dist_h(X,Z)\)，我们把X想象成一个发光体(就像太阳)，把Z想象成一个受光体(就像地球)。太阳上的某个发光点x发射的光线只要到达地球上任意一个点，我们即认为完成了光线传播，在Z上最先接收到光线的点\(z^*\)称之为受光点。于是，\(z^*=argmin_z|x-z|\)，将距离\(|x-z^*|\)称作为x点到Z的传播距离，那么\(dist_h(X,Z)=\max_x\min_z |x-z|\)表示X到Z的最远光线传播距离。
而\(dist_h(Z,X)\)考察的是Z到X的光线传播距离，此时，Z是发光体，X是受光体。\(dist_h(X,Z)\)和\(dist_h(Z,X)\)未必相等，两者中较大者即为豪斯多夫距离\(dist_H(X,Z)\)，代表着X和Z集合之间的最远光线传播距离。

豪斯多夫距离的非负性和对称性显而易见，下面只证明同一性和直递性。
同一性：如果两个集合的豪斯多夫距离为零，则两个集合相等。
\(dist_H(X,Z)=0\)→ \(dist_h(X,Z)=0\) →任意\(x\in X，min_z|x-z|=0\) →任意\(x\in X\)，都存在\(z=x\) →\(X\subseteq Z\)，同理有\(Z\subseteq X\)，因此X=Z。
直递性：(证明过程参考了这篇博文，在理解原文基础上进行了修改)

观察上图，将X，Y，Z故意画成不同的形状，为了方便绘图和讨论，这里假设集合内元素连续分布，其结果同样适用于离散元素分布的情况。\(dist_h(X,Z)\)对应的距离为\(|x_1-z_1|\)（\(x_1\)和\(z_1\)分别为发光点和受光点，后面类似），\(dist_h(X,Y)\)对应的距离为\(|x_2-y_1|\)，\(dist_h(Y,Z)\)对应的距离为\(|y_2-z_2|\)，可以证明\(dist_h(X,Z)≤dist_h(X,Y)+ dist_h(Y,Z)\)：
设\(x_1\)点到Y受光点为\(y_3\)，亦即\(x_1\)到Y的最短距离为\(|x_1-y_3|\)；类似的，\(y_3\)到Z的受光点为\(z_3\)，那么有：

\[\begin{aligned} dist_h(X,Z)&=|x_1-z_1|\\ &\leq|x_1-z_3|\\ &\leq|x_1-y_3|+|y_3-z_3|\\ &\leq|x_2-y_1|+|y_2-z_2|\\ &=dist_h(X,Y)+dist_h(Y,Z) \end{aligned}\]

上式中，第2行的不等式是由于\(|x_1-z_1|\)是\(x_1\)到Z的最近距离，第3行是应用了距离的三角不等式性质，第4行是根据定义，\(dist_h(A,B)\)是A到B的最远光线传播距离。
同理，可以证明 \(dist_h(Z,X)≤dist_h(Z,Y)+ dist_h(Y,X)\)。
不失一般性，假设\(dist_h(Z,X)≥dist_h(X,Z)\)，则有：
\(dist_H(Z,X)=dist_h(Z,X)≤dist_h(Z,Y)+ dist_h(Y,X)≤dist_H(Z,Y)+ dist_H(Y,X)\)
豪斯多夫距离的直递性得证。

9.3 试析k均值算法能否找到最小化式(9.24)的最优解。

答：k均值算法的运行结果依赖于初始选择的聚类中心，找到的结果是局部最优解，未必是全局最优解。

9.4 试编程实现k均值算法，设置三组不同的k值、三组不同初始中心点，在西瓜数据集4.0上进行实验比较，并讨论什么样的初始中心有利于取得好结果。

答：编程代码附后。
对于k值，考虑三种情况：2，3，4。
对于初始中心点，从直觉上判断，初始中心点应该尽量分散一些较好。结合西瓜数据集4.0的具体数据分布，通过手动选取方法来确定三组不同的初始点：
相互靠拢的几个点：5,7,17,18
相互分散，靠近数据区域外周的几个点：29,15,10,25
相互分散，位于数据区域中部的几个点：27,17,12,21。

运行结果如下(小黑点是数据点，蓝色圆圈是初始中心，红色十字点是最终得到的聚类中心，红色虚线是聚类边界)：

上图中第一列为初始聚类中心比较靠近，第二列的初始中心比较分散而靠近外周，第三列的初始中心分散而居中。
讨论：在k=2和3时，不同的聚类中心所得结果相差不明显，所得结果相近；k=4时，初始中心选择分散时，最终的平方误差(按9.24式计算)更小一些。因此，就当前所尝试的有限情况而言，貌似可以得到结论：当k较小时，不同的初始中心选取差别不大，当k较大时，选取，选择较分散的初始中心更有利于取得较好结果.

9.5 基于DBSCAN的概念定义，若\(x\)为核心对象，由\(x\)密度可达的所有样本构成的集合为\(X\)。试证明：\(X\)满足连接性(3.39)与最大性(9.40)。

证明：这个结论显然的，设\(x_i, x_j \in X\)，由于\(X\)是由所有\(x\)密度可达的样本构成，因此\(x_i, x_j\)都可由\(x\)密度可达，于是\(x_i, x_j\)密度可连，连接性成立。
若\(x_j\)由\(x_i\)密度可达，而\(x_j\)又由\(x\)密度可达，于是\(x_j\)由\(x\)密度可达，因此，\(x_j\in X\)，最大性成立。

9.6 试析AGNES算法使用最小距离和最大距离的区别。

答：
参考1 (from 百度文库)：

• 最小距离又叫单链接，容易造成一种叫做Chaining的效果，两个cluster明明从“大局”上离得比较远，但是由于其中个别的点距离比较近就被合并了，并且这样合并之后的- Chaining效应会进一步扩大，最后会得到比较松散的cluster。
• 最大距离又叫全链接，则是单链接的反面极端，其效果刚好相反，限制非常大，两个cluster即使已经很接近了，但是只要有不配合的点存在，就顽固到底，老死不相合并，也不是太好的办法。
• 这两种方法的共同问题就是只考虑了某个有特点的数据，而没有考虑类内数据的整体特点。
• 平均距离或者平均链接，相对能得到合适一点的结果。平均链接的一个变种就是取两两距离的中值(中位数)，与取均值相比更加能够解除个别偏离样本对结果的干扰。

参考2 (from 百度文库)：

• 单链技术可以处理非椭圆形状的簇，但对噪声和离群点很敏感。
• 全链接对噪声和离群点不敏感，但是可能使大的簇破裂，偏好球型簇。
• 平均距离则是一种单链与全链的折中算法。优势是对噪声和极端值影响小，局限是偏好球型簇。

个人理解：
前面两处参考中主要考虑了不同距离规则下受噪声数据的影响大小。下面不考虑噪声的影响，考虑不同距离规则下发生合并的趋势。
考虑下图所示的情况，这里有三个簇，有一个较大的簇和两个较小的簇，假如把它们看成是我国战国时期的秦、韩、魏三个诸侯国。

按最小距离算法，首先合并秦韩；按最大距离算法，首先合并韩魏。
因此，在最小距离算法下，趋向于“远交近攻”的蚕食原则，相互靠近的诸侯国首先进行合并；在最大距离算法下，大国反应较为迟钝，小国家之间优先进行“合纵抗衡”。

9.7 聚类结果中若每个簇都有一个凸包(包含簇样本的凸多面体)，且这些凸包不相交，则称为凸聚类。试析本章介绍的哪些聚类算法只能产生凸聚类，哪些能产生非凸聚类。

答：

1. 原型聚类(k-means，LVQ，高斯混合)所得结果均是凸聚类。比如，在k-means的聚类分界面线段必然为某两个聚类中心的中垂线，考虑如下图所示情况：
   
   聚类中心1的聚边界为非凸包，右上角存在凹口。设其中一个分界线段为中心点1和中心点2的中垂线，那么上图中蓝色区域距离中心点2更近，不可能划分到聚类1。因此，k-means只能产生凸聚类。
2. 基于密度聚类的DBSCAN能够产生非凸聚类，因为该算法基于样本的分布密度对密度相连的样本进行聚类，最终的结果通过聚类簇内的样本来表征。如果聚类簇内的样本分布呈非凸型分布，那么聚类结果也将是非凸包。比如，教材图9.10中的情形。
3. 属于层次聚类的AGNES也能产生非凸聚类，比如如下情形中，采用dmin进行聚类，结果会出现包围形式的非凸聚类结果。而当采用dmax时，趋向于产生凸型聚类结果。
   

9.8 试设计一个聚类性能度量指标，并与9.2节中的指标比较。（暂缺）

答：

9.9 试设计一个能用于混合属性的非度量距离。（暂缺）

答：

9.10 试设计一个能自动确定聚类数的改进k均值算法，编程实现并在西瓜数据集4.0上运行。

答：

• 方案1：最小化平方误差(9.24)式来找到最佳k值。
  然而，该方案不可行，极限思考一下，令\(k=m，u_i=x_i, i=1~\sim m\)，此时的平方误差为0，显然，这样就没有什么意义了。
• 方案2：在(9.24)式的基础上增加一个惩罚项作为最小化的目标函数，比如令\(J=erro+0.1*k=\sum_c\sum_{x\in C} |x-u_c|_2^2+0.1*k\)。
  试验一下，结果如下：
  
  观察可见，平方误差erro随k值增加而递减，以此来确定最佳k值将得到\(k_{best}=m\)的结论，与预期相符。
  而目标函数\(J\)曲线则呈上开口抛物线型，当k=4时，目标函数\(J\)最小。上图为某一次运行结果，多次运行结果较稳定，该方案可行，可以用来确定最佳k值。
• 方案3：以衡量聚类性能的内部指标DBI和DI指数来找到最佳k值，见(9.12)和(9.13)式，其中DBI指数越小越好，DI指数越大越好。在西瓜数据集4.0上试验一下。结果如下：
  
  观察上图，DBI曲线大体上呈抛物线型，k=4时最佳，多次运行较稳定，可用于确定最佳k值。
  DI指数大体上呈下开口抛物线，上图所示最佳k值为6，但是多次运行时，DI曲线表现不稳定，不太适合用于确定最佳k值。
  从DBI和DI的表达式可见，DI的计算式中\(d_{min}(C_i,C_j)\)表示两个簇之间的最短样本距离，受个别噪声样本的影响较大，而DBI计算式中\(avg(C_i)\)和\(avg(C_j)\)表示簇内平均距离，是个统计平均量，受噪声影响相对较小一些。

---

附：编程代码

习题9.4 (Python)

# -*- coding: utf-8 -*-
"""
Created on Mon Apr 13 13:00:25 2020

@author: Administrator
"""

import numpy as np
import matplotlib.pyplot as plt

#========首先编写两个函数==============
# k_means：实现k-means聚类算法
# points：用于绘图的辅助函数，根据样本点的坐标计算外围凸多边形的顶点坐标

def k_means(X,k,u0=None,MaxIte=30):
    # k-means聚类算法
    # 输入：
    #    X:样本数据，m×n数组
    #    k:聚类簇数目
    #    u0:初始聚类中心
    #    MaxIte:最大迭代次数
    # 输出:
    #    u:最终的聚类中心
    #    C:各个样本所属簇号
    #    erro:按(9.24)式计算的平方方差结果
    #    step:迭代次数
    
    m,n=X.shape   #样本数和特征数
    if u0==None:  #随机选取k个样本作为初始中心点
        u0=X[np.random.permutation(m)[:k],:] 
    u=u0.copy()
    step=0
    while True:
        step+=1
        u_old=u.copy()            #上一次迭代的中心点
        dist=np.zeros([m,k])      #存储各个样本到中心点的距离
        for kk in range(k):       #计算距离
            dist[:,kk]=np.sqrt(np.sum((X-u[kk])**2,axis=1))
        C=np.argmin(dist,axis=1)  #以距离最小的中心点索引号最为簇类号
        for kk in range(k):       #更新聚类中心
            u[kk]=X[C==kk,:].mean(axis=0)  
        if (u==u_old).all() or step>MaxIte:
            break                 #如果聚类中心无变化，或者超过最大
                                  #迭代次数，则退出迭代
    #=====计算平方误差(9.24)
    erro=0
    for kk in range(k):
         erro+=((X[C==kk]-u[kk])**2).sum()
    return u,C,erro,step

def points(X,zoom=1.2):
    # 为了绘制出教材上那种凸聚类簇效果
    # 本函数用于计算凸多边形的各个顶点坐标
    # 输入:
    #     X:簇类样本点的坐标
    #     zoom:缩放因子(最外围样本点向外扩展系数)
    # 输出:
    #     ps:凸多边形的顶点坐标
    
    X=X[:,0]+X[:,1]*1j  #将坐标复数化
    u=np.mean(X)        #聚类中心
    X=X-u               #原点移至聚类中心
    # 寻找凸多边形的各个顶点坐标
    indexs=[]                         #存储顶点坐标的索引号
    indexs.append(np.argmax(abs(X)))  #首先将距离距离中心最远的点作为起始顶点  
    while True:
        p=X[indexs[-1]]     #当前最新确定的顶点
        X1=1E-5+(X-p)/(-p)  #以p点为原点，并且以u-p为x轴(角度为0)
        #上式中加1E-5的小正数是因为p点自己减去自己的坐标有时候会出现
        #(-0+0j)的情况，angle(-0+0j)=-180°,但是希望结果为0
        indexs.append(np.argmax(np.angle(X1)))  #新找到顶点
        if indexs[-1]==indexs[0]:               #如果这些顶点首尾相连了，则停止
            break
    # 将复数坐标还原成直角坐标
    ps=np.c_[np.real(X)[indexs],np.imag(X)[indexs]]
    ps=ps*zoom+[np.real(u),np.imag(u)]
    return ps

#================================================
#                  主程序
#================================================

#==============西瓜数据集4.0======================
D=np.array(
        [[0.697,0.460],[0.774,0.376],[0.634,0.264],[0.608,0.318],[0.556,0.215],
         [0.403,0.237],[0.481,0.149],[0.437,0.211],[0.666,0.091],[0.243,0.267],
         [0.245,0.057],[0.343,0.099],[0.639,0.161],[0.657,0.198],[0.360,0.370],
         [0.593,0.042],[0.719,0.103],[0.359,0.188],[0.339,0.241],[0.282,0.257],
         [0.748,0.232],[0.714,0.346],[0.483,0.312],[0.478,0.437],[0.525,0.369],
         [0.751,0.489],[0.532,0.472],[0.473,0.376],[0.725,0.445],[0.446,0.459]])
m=D.shape[0]
#=============绘制样本数据点及其编号===============
plt.figure()
plt.scatter(D[:,0],D[:,1],marker='o',s=250,c='',edgecolor='k')
for i in range(m):
    plt.text(D[i,0],D[i,1],str(i),
             horizontalalignment='center',
             verticalalignment='center')
plt.show()

#选取三种不同的初始聚类中心点
centers=np.array([[5,7,17,18],     #相互靠近的点
                  [29,15,10,25],   #分散而外周的点
                  [27,17,12,21]])  #分散而中间的点

#======运行k-means聚类，设置三组不同的k值、三组不同初始中心点=======
for i,k in enumerate([2,3,4]):    #不同的k值
    for j in range(3):            #3次不同初始值
        #=====k-means算法
        u0=D[centers[j][:k],:]
        u,C,erro,step=k_means(D,k,u0)
        #=====画图======
        plt.subplot(3,3,i*3+j+1)
        plt.axis('equal')
        plt.title('k=%d,step=%d,erro=%.2f'%(k,step,erro))
        # 画样本点
        plt.scatter(D[:,0],D[:,1],c='k',s=1)
        # 画聚类中心
        plt.scatter(u0[:,0],u0[:,1],marker='o',c='',s=50,edgecolors='b')
        plt.scatter(u[:,0],u[:,1],marker='+',c='r',s=50)#,c='',s=80,edgecolors='g')
        for kk in range(k):
            plt.annotate('',xy=u[kk],xytext=u0[kk],arrowprops=dict(arrowstyle='->'),ha='center')
        # 画聚类边界
        for kk in range(k):
            ps=points(D[C==kk])
            plt.plot(ps[:,0],ps[:,1],'--r',linewidth=1)
plt.show()

习题9.10 (Python)

# -*- coding: utf-8 -*-
"""
Created on Mon Apr 13 13:00:25 2020

@author: Administrator
"""

import numpy as np
import matplotlib.pyplot as plt

#========首先编写两个函数==============
# k_means：实现k-means聚类算法
# points：用于绘图的辅助函数，根据样本点的坐标计算外围凸多边形的顶点坐标

def k_means(X,k,u0=None,MaxIte=30):
    # k-means聚类算法
    # 输入：
    #    X:样本数据，m×n数组
    #    k:聚类簇数目
    #    u0:初始聚类中心
    #    MaxIte:最大迭代次数
    # 输出:
    #    u:最终的聚类中心
    #    C:各个样本所属簇号
    #    erro:按(9.24)式计算的平方方差结果
    #    step:迭代次数
    
    m,n=X.shape   #样本数和特征数
    if u0==None:  #随机选取k个样本作为初始中心点
        u0=X[np.random.permutation(m)[:k],:] 
    u=u0.copy()
    step=0
    #print('D维度',X.shape,'u=',u)
    while True:
        
        step+=1
        u_old=u.copy()            #上一次迭代的中心点
        dist=np.zeros([m,k])      #存储各个样本到中心点的距离
        for kk in range(k):       #计算距离
            dist[:,kk]=np.sqrt(np.sum((X-u[kk])**2,axis=1))
        C=np.argmin(dist,axis=1)  #以距离最小的中心点索引号最为簇类号
        #print(u,C)
        for kk in range(k):       #更新聚类中心
            u[kk]=X[C==kk,:].mean(axis=0)  
        if (u==u_old).all():      #如果聚类中心无变化，则退出迭代
            break
        if step>MaxIte:           #如果超过最大迭代次数，则退出迭代
            print('超过最大迭代次数')
            break
    #=====计算平方误差(9.24)
    erro=0
    for kk in range(k):
         erro+=((X[C==kk]-u[kk])**2).sum()
    return u,C,erro,step


def Inner_Index(X,u):
    #计算衡量聚类性能的内部指标DBI和DI指数
    m,n=X.shape
    k=u.shape[0]
    dist=np.zeros([m,k])      #存储各个样本到中心点的距离
    for kk in range(k):       #计算距离
        dist[:,kk]=np.sqrt(np.sum((X-u[kk])**2,axis=1))
    C=np.argmin(dist,axis=1)  #以距离最小的中心点索引号最为簇类号
    #===计算avg和diam
    avg=[]                    #存储簇内平均距离
    diam=[]                   #存储簇内最远距离
    for kk in range(k):
        Xkk=X[C==kk]          #第kk个聚类簇内的样本点
        mk=len(Xkk)
        if mk<=1:
            avg.append(0)
            diam.append(0)
            continue
        dist=[]               #存储簇内两个不同样本之间的距离
        for i in range(mk):
            for j in range(i+1,mk):
                dist.append(np.sqrt(((Xkk[i]-Xkk[j])**2).sum()))
        avg.append(np.mean(dist))
        #print('dist=',dist,'mk=',mk)
        diam.append(np.max(dist))
    #===计算DB指数
    DBIij=np.zeros([k,k])    #存储 [avg(Ci)+avg(Cj)]/dcen(ui,uj)
    for i in range(k):
        for j in range(i+1,k):
            DBIij[i,j]=(avg[i]+avg[j])/np.sqrt(((u[i]-u[j])**2).sum())
            DBIij[j,i]=DBIij[i,j]
    DBI=np.mean(np.max(DBIij,axis=1))
    #===计算Dunn指数
    dmin=[]    #存储Ci和Cj之间的最短距离
    for i in range(k):
        for j in range(i+1,k):
            Xi=X[C==i]
            Xj=X[C==j]
            dmin_ij=[]
            for xi in Xi:
                for xj in Xj:
                    dmin_ij.append(np.sqrt(((xi-xj)**2).sum()))
            dmin.append(min(dmin_ij))
    DI=min(dmin)/max(diam)
    
    return DBI,DI


#================================================
#                  主程序
#================================================

#==============西瓜数据集4.0======================
D=np.array(
        [[0.697,0.460],[0.774,0.376],[0.634,0.264],[0.608,0.318],[0.556,0.215],
         [0.403,0.237],[0.481,0.149],[0.437,0.211],[0.666,0.091],[0.243,0.267],
         [0.245,0.057],[0.343,0.099],[0.639,0.161],[0.657,0.198],[0.360,0.370],
         [0.593,0.042],[0.719,0.103],[0.359,0.188],[0.339,0.241],[0.282,0.257],
         [0.748,0.232],[0.714,0.346],[0.483,0.312],[0.478,0.437],[0.525,0.369],
         [0.751,0.489],[0.532,0.472],[0.473,0.376],[0.725,0.445],[0.446,0.459]])
m=D.shape[0]

#======考察平方误差(9.24式)随k值的变化规律=======
# 结果如预期，平方误差随k值增加而减小
# 如果加上一个惩罚项 0.1*k，出现抛物线
erros=[]
ks=range(2,8+1)         #考察k=2~8
for k in ks:    
    erro=[]
    for i in range(5):  #对于每个k值，初始5次聚类，取其中最佳者
        u,C,erro_i,step=k_means(D,k)
        erro.append(erro_i)
    erros.append(min(erro))
plt.figure()
plt.plot(ks,erros,'o-')
plt.plot(ks,0.1*np.array(ks)+np.array(erros),'o-')
plt.xlabel('k')
plt.legend(['erro','erro+0.1*k'])
plt.show()


#======考察DBI评分(9.12式)随k值的变化规律=======
DBIs=[]
DIs=[]
ks=range(2,8+1)         #考察k=2~8
for k in ks:    
    DBI=[]
    DI=[]
    for i in range(5):  #对于每个k值，初始5次聚类，取其中最佳者
        u,C,erro,step=k_means(D,k)
        DBi,Di=Inner_Index(D,u)
        DBI.append(DBi)
        DI.append(Di)
    DBIs.append(min(DBI))
    DIs.append(max(DI))
# 绘制计算结果
ax1=plt.figure().add_subplot(111)
ax2=ax1.twinx()    #双y坐标

ax1.plot(ks,DBIs,'o-k',label='DBI')
ax2.plot(ks,DIs,'o-r',label='DI')

ax1.set_xlabel('k')
ax1.set_ylabel('DBI')
ax2.set_ylabel('DI')

ax1.legend(loc=(0.5,.9))
ax2.legend(loc=(0.5,.8))

plt.show()