西瓜书习题试答-第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.1 试证明：p≥1时，闵可夫斯基距离满足距离度量的四条基本性质；0≤p<1时，闵可夫斯基距离不满足直递性，但满足非负性、同一性、对称性；p趋向无穷大时，闵可夫斯基距离等于对应分量的最大绝对距离，即
\]
证明：
-
非负性、同一性、对称性
对于所有的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)计算：
\]
其中
\]
试证明：豪斯多夫距离满足距离度量的四条基本性质。
证明：首先，我们来直观的理解一下豪斯多夫距离。
下面表示两个点的欧式距离时，为了简便，右下角的角标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\)，那么有：
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 聚类结果中若每个簇都有一个凸包(包含簇样本的凸多面体)，且这些凸包不相交，则称为凸聚类。试析本章介绍的哪些聚类算法只能产生凸聚类，哪些能产生非凸聚类。
答：
- 原型聚类(k-means，LVQ，高斯混合)所得结果均是凸聚类。比如，在k-means的聚类分界面线段必然为某两个聚类中心的中垂线，考虑如下图所示情况：
聚类中心1的聚边界为非凸包，右上角存在凹口。设其中一个分界线段为中心点1和中心点2的中垂线，那么上图中蓝色区域距离中心点2更近，不可能划分到聚类1。因此，k-means只能产生凸聚类。 - 基于密度聚类的DBSCAN能够产生非凸聚类，因为该算法基于样本的分布密度对密度相连的样本进行聚类，最终的结果通过聚类簇内的样本来表征。如果聚类簇内的样本分布呈非凸型分布，那么聚类结果也将是非凸包。比如，教材图9.10中的情形。
- 属于层次聚类的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()