一、K-Means聚类及其优化算法

1.1 传统K-Means聚类算法

1.1.1 原理

• 算法流程：
  输入：样本集合$D=x_1,x_2,...,x_n$，聚类簇数$k$， 最大迭代次数$N$；
  输出：簇划分$C=C_1,C_2,...,C_k$；
  1）样本集中随机选择$k$个样本作为初始的聚类中心点：${\mu }_1,{\mu }_2,...,{\mu }_k$
  2）$fori=1toN$：
    a）将簇$C_t$初始化为空，$t=1,2,...,k$；
    b）对于样本$x_i,i=1,2,...,n$，计算和聚类中心${\mu }_j,j=1,2,..,k$之间的距离：$d_{ij}=\Vert x_i-{\mu }_j{\Vert }_2^2$，将$x_i$标记为最小的$d_{ij}$所对应的类别${\lambda }_j,j=1,2,...,k$，此时更新 C λ j = C λ j ∪ { x i } C_{\lambda_j}=C_{\lambda_j}\cup \{x_i\} Cλj​​=Cλj​​∪{xi​}；
    c）对于$j=1,2,...,k$，计算$C_j$中聚类中心点：${\mu }_j=\frac{1}{|C_j|}{\sum }_{x_i\in C_j}{x}_i$，其中$i=1,2,...,n$；
    e）直到达到最大迭代次数或者聚类中心点不再发生变化，转到步骤3）
  3）输出簇划分$C=C_1,C_2,...,C_k$

1.1.2 K_Means算法优缺点

1. 优点：
   1）原理简单，易容易；
   2）可解释性较强；
   3）时间复杂度是是$O(n\times k\times t)$，近似于线性，而且适合挖掘大规模数据集。其中n代表数据集中样本的数量，N代表着算法迭代的次数，k代表着簇的数目。
2. 缺点：
   1）易受初始聚类中心的选择和离群点影响，结果不稳定；
   2）无法很好地解决簇分布差别较大的情况（如一类是另一类数量的100倍）；
   3）聚类簇数$k$需事先给定，改值与真实分布未必吻合；
   4）受噪声影响较大。

1.1.3 K_Means算法调优及k值选取

1. 数据预处理：归一化，统一单位；离群点处理。
2. $k$值选择：
   1）手肘法：拐点就是$k$的最佳值，下图中x轴为聚类数$k$的值，y轴为误差平方和$SSE$，$SSE={\sum }_{i=1}^n{\sum }_{x_i\in C_j}||x_i-{\mu }_j|{|}^2,j=1,2,...,k$，拐点为4，所以4为聚类簇数。
   
   实现：

#导入库
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans   

#导入数据
df=pd.read_csv(r'E:\data analysis\test\cluster.txt',header=None,sep='\s+')
print(df.head())
x=df.iloc[:,:-1]  
y=df.iloc[:,-1]

#肘方法看k值
d=[]
for i in range(1,11):    #k取值1~11，做kmeans聚类，看不同k值对应的簇内误差平方和
    km=KMeans(n_clusters=i,init='k-means++',random_state=0)
    km.fit(x)
    d.append(km.inertia_)  #inertia簇内误差平方和

plt.plot(range(1,11),d,marker='o')
plt.xlabel('number of clusters')
plt.ylabel('distortions')
plt.show()

2）Gap Statistic方法：找到最大的Gap Statistic值对应的$k$值
$$Gap(k)=E(\log D_k)-\log D_k$$
其中$E(\log D_k)$为随机样本损失的期望，是$\log (D_k)$的期望，$\log D_k$为实际样本的损失，$D_k$为损失函数，是随机样本的损失与实际样本的损失之差。
$E(\log D_k)$的计算方法如下：在样本所在的区域内按均匀分布随机地产生和原始数据一样多的随机样本，并对于随机样本做K-Means聚类，得到一个$D_k$，重复多次计算得出$E(\log D_k)$的近似值。
如果实际样本对应的最佳簇数为$k$， 那么$\log D_k$会很小，即实际样本损失会很小，从而随机样本损失与实际样本损失之差相应达到最大，即$Gap(k)$取得最大值时对应的$k$值就是最佳的簇数。
3）采用核函数
核k均值算法：是核聚类方法的一种。通过非线性映射，将输入空间中的数据点映射到高维特征空间中，并在新的特征空间中聚类。
非线性映射增加了数据点线性可分的概率，从而经典的聚类算法失效时，引入核函数可以达到更准确的聚类结果。

1.2 K-Means++算法

1.2.1 原理

• 由于K-Means算法受初始聚类中心的影响较大，所以采用K-Means++算法来进行解决，基本思想是：初始聚类中心点之间的距离尽可能的远。
• 算法流程：
  1）在输入的数据集中随机选择一个样本点作为初始的聚类中心${\mu }_1$；
  2）对于样本中的每一个点$x_i$，计算它与已选择的聚类中心点之间的距离$D(x_i)=\Vert x_i-{\mu }_r{\Vert }_2^2,r=1,2,..,k_{selected}$；
  3）选择一个新的数据点作为新的聚类中心，选择的原则是：$D(x)$较大的点，被选取作为新的聚类中心的概率较大；
  4）重复步骤2）和3）直到选择出$k$个聚类中心点；
  5）利用这$k$个聚类中心点去执行标准的K-Means聚类算法。

1.2.2 代码实现

• KMeans参数：
  n_clusters：整形，聚类数，默认值是8；
  max_iter：整形，缺省值=300，执行一次k-means算法所进行的最大迭代数；
  init：三个可选值：’k-means++’（默认值，一种特殊的方法选定初始质心从而能加速迭代过程的收敛）， ‘random’（随机从训练数据中选取初始质心），或者传递一个ndarray向量（形如 (n_clusters, n_features)并给出初始质心）；
  n_jobs：整形数，计算所用的进程数；
  random_state：整形或 numpy.RandomState 类型，可选
  用于初始化质心的生成器（generator）。如果值为一个整数，则确定一个seed。此参数默认值为numpy的随机数生成器；
• KMeans属性：
  cluster_centers_：向量，[n_clusters, n_features] (聚类中心的坐标)；
  labels_: 每个点的分类标签；
  inertia_：float型，每个点到其簇类中心的距离之和。
• KMeans方法：
  fit(X)：计算k-means聚类。
  fit_predict(X)：　计算聚类中心并给每个样本预测类别。
  fit_transform(X)：计算簇并 transform X to cluster-distance space。
  get_params([deep])：取得估计器的参数。
  predict(X)：给每个样本估计最接近的簇。
  score(X)：计算聚类误差
  **set_params(params)：为这个估计器手动设定参数。
  transform(X[,y])：将X转换为群集距离空间。在新空间中，每个维度都是到集群中心的距离。 请注意，即使X是稀疏的，转换返回的数组通常也是密集的。

from sklearn.cluster import KMenas
import numpy as np
import matplotlib.pyplot as plt
# 要分类的数据点
x = np.array([ [1,2],[1.5,1.8],[5,8],[8,8],[1,0.6],[9,11] ])

# 把上面数据点分为两组（非监督学习）
clf = KMeans(n_clusters=2).fit(x)
centers = clf.cluster_centers_ # 两组数据点的中心点
labels = clf.labels_   # 每个数据点所属分组
print(centers)
print(labels)

for i in range(len(labels)):
    plt.scatter(x[i][0], x[i][1], c=('r' if labels[i] == 0 else 'b')) #画出样本点，参数c是颜色
plt.scatter(centers[:,0],centers[:,1],marker='*', s=100) #画出聚类中心，s是散点图点的大小，默认为20
# 预测
predict = [[2,1], [6,9]]
label = clf.predict(predict)
for i in range(len(label)):
    plt.scatter(predict[i][0], predict[i][1], c=('r' if label[i] == 0 else 'b'), marker='x') #marker：指定散点图点的形状，默认为空心圆
pyplot.show()

*是两组数据的”中心点”；x是预测点分组。

1.3 K-Means距离计算优化elkan K-Means

• 传统的K-Means聚类算法中，每轮迭代时，要计算所有的样本到类中心的距离。elkan K-Means在距离方面进行了改进。
• elkan K-Means是利用两边之和大于等于第三边，两边之差小于第三边的性质来减少距离计算。
• 1）对于一个样本点和两个聚类中心${\mu }_1$和${\mu }_2$，我们预先计算出两个聚类中心之间的距离为$D({\mu }_1,{\mu }_2)$，如果计算发现$2D(x,{\mu }_1)\le D({\mu }_1,{\mu }_2)$，我们即可知道$D(x,{\mu }_1)\le D(x,{\mu }_2)$，因此可以省去计算$x$到${\mu }_2$之间的距离。
• 2）对于一个样本点和两个聚类中心${\mu }_1$和${\mu }_2$，我们可以得到 D ( x , μ 2 ) ≥ m a x { 0 , D ( x , μ 1 ) − D ( μ 1 , μ 2 ) } D(x, \mu_2)\ge max\{0, D(x, \mu_1)-D(\mu_1, \mu_2)\} D(x,μ2​)≥max{0,D(x,μ1​)−D(μ1​,μ2​)}
• 利用上边的两个规律，elkan K-Means比起传统的K-Means迭代速度有很大的提高。但是如果我们的样本的特征是稀疏的，有缺失值的话，这个方法就不使用了，此时某些距离无法计算，则不能使用该算法。

1.4 大样本优化Mini Batch K-Means

• 在统的K-Means算法中，要计算所有的样本点到所有的质心的距离。如果样本量非常大，比如达到10万以上，特征有100以上，此时用传统的K-Means算法非常的耗时，就算加上elkan K-Means优化也依旧。在大数据时代，这样的场景越来越多。此时Mini Batch K-Means应运而生。
• Mini Batch，也就是用样本集中的一部分的样本来做传统的K-Means，这样可以避免样本量太大时的计算难题，算法收敛速度大大加快。当然此时的代价就是我们的聚类的精确度也会有一些降低。一般来说这个降低的幅度在可以接受的范围之内。
• 在Mini Batch K-Means中，我们会选择一个合适的批样本大小batch size，我们仅仅用batch size个样本来做K-Means聚类。那么这batch size个样本怎么来的？一般是通过无放回的随机采样得到的。
• 为 了增加算法的准确性，我们一般会多跑几次Mini Batch K-Means算法，用得到不同的随机采样集来得到聚类簇，选择其中最优的聚类簇。

1.5 K-Means与KNN区别与联系

• 区别：
  1）K-Means是无监督学习的聚类算法，没有样本输出；而KNN是有监督学习的分类算法，有对应的类别输出。
  2）KNN基本不需要训练，对测试集里面的点，只需要找到在训练集中最近的k个点，用这最近的k个点的类别来决定测试点的类别。而K-Means则有明显的训练过程，找到k个类别的最佳质心，从而决定样本的簇类别。
• 联系：
  1）两个算法都要找出和某一个点最近的点，即两者都利用了最近邻(nearest neighbors)的思想。

二、ISODATA聚类算法（迭代自组织数据分析法）：k值不确定时，可以使用此法

• 在K-Means中，聚类数k值需要人为指定，并在算法过程中无法更改，当遇到高维。海量的数据集时，一般无法准确的估计出k值。
• ISODATA法针对此问题进行 了改进：当属于某个类别的样本数过少时，把该类别去除；当属于某个类别的样本数过多、分散程度较大时，把该类别分为2个类别。
• ISODATA算法在Ｋ-Means算法基础上增加两个操作：
  １）分裂操作：增加聚类中心数
  ２）合并操作：减少聚类中心数
• ISODATA算法过程：
  １）预期的聚类中心数$k_0$，在ISODATA算法中，聚类数可以改变，$k_0$是用户指定的参考值，聚类中心数目的变动范围也由其决定，变动范围为$[\frac{1}{2}{k}_0,2k_0]$；
  2）每个类所要求的最小样本数目为$N_{min}$：分裂后样本数小于$N_{min}$时，就不对该类别进行分类操作；
  3）最大方差$\sigma$：用于控制类中样本分散的程度。当样本分散程度大于$\sigma$时，且分裂后满足条件2），进行分裂操作；
  4）两个聚类中心之间允许最小距离为$D_{min}$：若两个聚类中心之间的距离小于$D_{min}$时，则对两个类进行合并。
• 若希望样本不划分到单一类中，可使用模糊聚类或者高斯混合聚类。

三、高斯混合聚类（GMM）

3.1 原理

• 假设：每个簇的数据分布都符合高斯分布，当前数据呈现的分布就是各个簇高斯分布叠加在一起的结果。

• 对于n维样本空间的随机变量x，若服从高斯分布，则分布函数为：
  $$P(x|\mu ,\Sigma )=\frac{1}{(2\pi {)}^{\frac{n}{2}}|\Sigma {|}^{\frac{1}{2}}}{e}^{-\frac{1}{2}(x-\mu {)}^T{\Sigma }^{-1}(x-\mu )}$$
  定义高斯混合分布为：
  $$P_M(x)=\sum _{i=1}^K{\alpha }_{i}P(x|{\mu }_i,{\Sigma }_i)\text{\hspace{1em}(1)}$$
  每个单独的分模型都是标准的高斯模型，其均值和方差为${\mu }_i,{\Sigma }_i$都是待估计参数。
  ${\alpha }_i$是“混合系数”，满足${\sum }_{i=1}^K{\alpha }_i=1$

• 对于训练集$D=x_1.x_2,...,x_m$，令随机变量$Z_j\in 1,2,...,K$，表示生成样本$x_j$的高斯混合成分，且先验概率为$P(Z_j=i)={\alpha }_i,i=1,2,,,.K$。
  $Z_j$的后验概率为：
  $$\begin{array}{rl}{P}_M(Z_j=i|x_j)& =\frac{P_M(x_j|Z_j=i)P(Z_j=i)}{P_M(x_j)}\\ & =\frac{{\alpha }_{i}P(x_j|{\mu }_i,{\Sigma }_i)}{{\sum }_{l=1}^K{\alpha }_{l}P(x_j|{\mu }_l,{\Sigma }_l)}\\ & ={\gamma }_{ji}\end{array}\text{\hspace{1em}(2)}$$
  $P_M(Z_j=i|x_j)$为给定样本$x_j$的条件下，此样本由第$i$个高斯混合成分生成的后验概率。
  样本集$D$有$K$个簇， C = { C 1 , C 2 , . . . , C K } C=\{C_1, C_2, ..., C_K\} C={C1​,C2​,...,CK​}，每个样本$x_j$的簇划分为：
  $$\begin{array}{r}{\lambda }_j=\underset{i=1,2,...,K}{\mathrm{argmax}}{\gamma }_{ji}将样本x_j划分到后验概率最大的那个类中\end{array}\text{\hspace{1em}(3)}$$
  对参数${\alpha }_i,{\mu }_i,{\Sigma }_i$的求解，采用极大似然法。
  $$\begin{array}{rl}LL(D)& =\sum _{j=1}^m\log P_M(x_j)\\ & =\sum _{j=1}^m\log (\sum _{i=1}^K{\alpha }_{i}P(x_j|{\mu }_i,{\Sigma }_i))\end{array}\text{\hspace{1em}(4)}$$
  EM算法：
  式（4）对${\mu }_i$求偏导数，有：
  $$\frac{\partial LL(D)}{\partial {\mu }_i}=0⟹\sum _{j=1}^m\frac{{\alpha }_{i}P(x_j|{\mu }_i,{\Sigma }_i)}{{\sum }_{l=1}^K{\alpha }_{l}P(x_j|{\mu }_l,{\Sigma }_l)}(x_j-{\mu }_i)=0⟹\sum _{j=1}^m{\gamma }_{ji}(x_j-{\mu }_i)=0\text{\hspace{1em}(5)}$$
  所以可得：
  $${\mu }_i=\frac{{\sum }_{j=1}^m{\gamma }_{ji}{x}_j}{{\sum }_{j=1}^m{\gamma }_{ji}}\text{\hspace{1em}(6)}$$
  即各混合成分的均值，通过样本的加权平均来估计，样本权重是每个样本属于该成分的后验概率。
  同理：
  $$\frac{\partial LL(D)}{\partial {\Sigma }_i}=0⟹{\Sigma }_i=\frac{{\sum }_{j=1}^m{\gamma }_{ji}(x_j-{\mu }_i)(x_j-{\mu }_i{)}^T}{{\sum }_{j=1}^m{\gamma }_{ji}}\text{\hspace{1em}(7)}$$
  对于混合系数${\alpha }_i$，为了使$LL(D)$最大，需满足${\alpha }_i\ge 0,{\sum }_{i=1}^K=1$，采用拉格朗日乘数法，即
  $$LL(D)+\lambda (\sum _{i=1}^K-1),{\lambda }_i为拉格朗日乘子\text{\hspace{1em}(8)}$$
  对${\alpha }_i$求偏导等于0，可得：
  $$\sum _{j=1}^m\frac{P(x_j|{\mu }_i,{\Sigma }_i)}{{\sum }_{l=1}^K{\alpha }_{l}P(x_j|{\mu }_l,{\Sigma }_l)}+\lambda =0\text{\hspace{1em}(9)}$$
  令:
  $$\varphi =\frac{P(x_j|{\mu }_i,{\Sigma }_i)}{{\sum }_{l=1}^K{\alpha }_{l}P(x_j|{\mu }_l,{\Sigma }_l)}\text{\hspace{1em}(10)}$$
  由式(9)可得：
  $$\sum _{j=1}^m\sum _{i=1}^K{\alpha }_i\varphi +\sum _{i=1}^K{\alpha }_i\lambda =0\text{\hspace{1em}(11)}$$
  其中：
  $$\sum _{i=1}^K{\alpha }_i\varphi =1和\sum _{i=1}^K{\alpha }_i=1\text{\hspace{1em}(12)}$$
  所以式(11)可以化解为：
  $$m+\lambda =0\text{\hspace{1em}(13)}$$
  所以对式(9)同时乘以${\alpha }_i$可得：
  $$\sum _{i=1}^m{\gamma }_{ji}+\lambda {\alpha }_i=0⟹\sum _{i=1}^m{\gamma }_{ji}-m{\alpha }_i=0\text{\hspace{1em}(14)}$$
  所以可得：
  $${\alpha }_i=\frac{1}{m}\sum _{j=1}^m{\gamma }_{ji}\text{\hspace{1em}(15)}$$
  ${\alpha }_i$由样本属于该成分的平均后验概率确定。

总结：
1）E-step：根据当前参数，计算每个样本属于高斯成分的后验概率${\gamma }_{ji}$
2）M-step：根据${\gamma }_{ji}$，更新参数${\mu }_i,{\Sigma }_i,{\alpha }_i$
3）停止条件：达到最大迭代次数或者似然函数$LL(D)$增长很小或不增长。

3.2 代码实现

import numpy as np
from sklearn.mixture import GaussianMixture
X = np.array([[1, 2], [1, 4], [1, 0], [10, 2], [10, 4], [10, 0]])
gm = GaussianMixture(n_components=2, random_state=0).fit(X)
gm.means_ #输出每个高斯混合模型的中心点
array([[10.,  2.],
       [ 1.,  2.]])
gm.predict([[0, 0], [12, 3]]) #对样本进行预测
array([1, 0])

参考网址：https://scikit-learn.org/stable/modules/generated/sklearn.mixture.GaussianMixture.html

四、DBSCAN聚类

4.1 原理

4.1.1 DBSCAN密度定义

• DBSCAN（Density-Based Spatial Clustering of Applications with Noise，具有噪声的基于密度的聚类方法）：是基于一组领域来描述样本集的紧密程度的，参数($ϵ,MinPts)$来猫叔领域的样本分布紧密程度，其中$ϵ$描述了某一样本的领域距离的阈值，$MinPts$描述了某一样本的距离为$ϵ$的领域中样本个数的阈值。

• 假设样本集是$D=(x_1,x_2,...,x_m)$，则DBSCAN具体的密度描述定义如下：
  1）$ϵ$-邻域：对于$x_j\in D$，其$ϵ$-邻域包含样本集$D$中与$x_j$的距离不大于$ϵ$的子样本集，即 N ϵ ( x j ) = { x i ∈ D ∣ d i s t a n c e ( x i , x j ) ≤ ϵ } N_{\epsilon}(x_j) = \{x_i \in D | distance(x_i,x_j) \leq \epsilon\} Nϵ​(xj​)={xi​∈D∣distance(xi​,xj​)≤ϵ}，这个子样本集的个数记为$|N_{ϵ}(x_j)|$；
  2）核心对象：对任一样本$x_j\in D$，如果其$ϵ$-邻域对应的$|N_{ϵ}(x_j)|$至少包含$MinPts$个样本，即如果$|N_{ϵ}(x_j)|\ge MinPts$，则$x_j$是核心对象；
  3）密度直达：如果$x_i$位于$x_j$的$ϵ$-邻域内，且$x_j$是核心对象，则称$x_i$由$x_j$密度直达。反之不一定成立，即此时不能说$x_j$由$x_i$密度直达，除非$x_i$也是核心对象，且$x_j$位于$x_i$的$ϵ$-邻域内；
  4）密度可达：对于$x_i$和$x_j$，如果存在样本序列$p_1,p_2,...,p_T$，满足$p_1=x_i,p_T=x_j$，且$p_{t+1}$由$p_t$密度直达，则称$x_j$由$x_i$密度可达，即密度可达有传递性，此时序列中传递样本$p_1,p_2,...,p_{T-1}$均为核心对象，只有核心对象才能使其他样本密度直达。注意密度可达也不满足对称性；
  5）密度相连：对于$x_i$和$x_j$，如果存在核心对象样本$x_k$，使得$x_i$和$x_j$均由$x_k$密度可达，则称$x_i$和$x_j$密度相连，注意密度相连满足对称性；

• 下图可以很容易看出理解上述定义，图中MinPts=5，红色的点都是核心对象，因为其$ϵ$-邻域至少有5个样本。黑色的样本是非核心对象。所有核心对象密度直达的样本在以红色核心对象为中心的超球体内，如果不在超球体内，则不能密度直达。图中用绿色箭头连起来的核心对象组成了密度可达的样本序列。在这些密度可达的样本序列的$ϵ$-邻域内所有的样本相互都是密度相连的。
  

4.1.2 DBSCAN密度聚类思想

• DBSCAN密度聚类定义：由密度可达关系导出的最大密度相连的样本集合，即为我们最终聚类的一个类别（簇）。
• 这个DBSCAN的簇里面可以有一个或者多个核心对象。如果只有一个核心对象，则簇里其他的非核心对象样本都在这个核心对象的$ϵ$-邻域里；如果有多个核心对象，则簇里的任意一个核心对象的$ϵ$-邻域中一定有一个其他的核心对象，否则这两个核心对象无法密度可达。这些核心对象的$ϵ$-邻域里所有的样本的集合组成的一个DBSCAN聚类簇。
• 如何找到簇样本集合呢？DBSCAN使用的方法很简单，它任意选择一个没有类别的核心对象作为种子，然后找到所有这个核心对象能够密度可达的样本集合，即为一个聚类簇。接着继续选择另一个没有类别的核心对象去寻找密度可达的样本集合，这样就得到另一个聚类簇。一直运行到所有核心对象都有类别为止。
• 还需考虑三个问题：
  1）第一个：一些异常样本点或者说少量游离于簇外的样本点，这些点不在任何一个核心对象在周围，在DBSCAN中，我们一般将这些样本点标记为噪音点。
  2）第二个：即如何计算某样本和核心对象样本的距离。在DBSCAN中，一般采用最近邻思想，采用某一种距离度量来衡量样本距离，比如欧式距离。这和KNN分类算法的最近邻思想完全相同。对应少量的样本，寻找最近邻可以直接去计算所有样本的距离，如果样本量较大，则一般采用KD树或者球树来快速的搜索最近邻。
  3）第三种：某些样本可能到两个核心对象的距离都小于$ϵ$，但是这两个核心对象由于不是密度直达，又不属于同一个聚类簇，那么如果界定这个样本的类别呢？一般来说，此时DBSCAN采用先来后到，先进行聚类的类别簇会标记这个样本为它的类别。也就是说DBSCAN的算法不是完全稳定的算法。

4.1.3 DBSCAN密度聚类算法

• 算法总结：
  输入：样本集$D=(x_1,x_2,...,x_m)$，邻域参数$(ϵ,MinPts)$，样本距离度量方式；
  输出：簇划分$C$；
  1）初始化核心对象集合$\Omega =\varnothing$，初始化聚类簇数$k$，初始化未访问样本集合$\Gamma =D$，簇划分$C=\varnothing$；
  2）对于$j=1,2,...,m$：
  　　　　a）通过距离度量方式，找到样本$x_j$的$ϵ$-邻域样本集$N_{ϵ}(x_j)$；
  　　　　b）如果子样本集样本个数满足$|N_{ϵ}(x_j)|\ge MinPts$，将样本$x_j$加入核心对象样本集合： Ω = Ω ∪ { x j } \Omega=\Omega\cup \{x_j\} Ω=Ω∪{xj​}
  3）如果核心对象集$\Omega =\varnothing$，则算法结束，否则转到步骤4）；
  4）在核心对象集合$\Omega$中，随机选择一个核心对象$o$，初始化当前核心对象 Ω c u r = { o } \Omega_{cur} = \{o\} Ωcur​={o}，初始化类别序号$k=k+1$，初始化当前簇样本集合 C k = { o } C_k=\{o\} Ck​={o}，更新未访问样本集合 Γ = Γ − { o } \Gamma=\Gamma-\{o\} Γ=Γ−{o}；
  5）如果当前簇核心对象${\Omega }_{cur}=\varnothing$，则当前聚类簇$C_k$生成完毕，更新簇划分 C = { C 1 , C 2 , . . . , C k } C=\{C_1,C_2,...,C_k\} C={C1​,C2​,...,Ck​}，更新核心对象集合$\Omega =\Omega -C_k$，转入步骤3）。否则更新核心对象集合$\Omega =\Omega -C_k$；
  6）在当前簇核心对象${\Omega }_{cur}$中取出$o^{{}^{\prime }}$，通过邻域距离阈值$ϵ$找出所有的$ϵ$-邻域子样本集$N_{ϵ}(o^{{}^{\prime }})$，令$\Delta =N_{ϵ}(o^{{}^{\prime }})\cap \Gamma$，更新当前簇样本集合$C_k=C_k\cup \Delta$，更新${\Omega }_{cur}={\Omega }_{cur}\cup (\Delta \cap \Omega )-o^{\prime }$，转入步骤5）；
  输出结果为：簇划分 C = { C 1 , C 2 , . . . , C k } C=\{C_1,C_2,...,C_k\} C={C1​,C2​,...,Ck​}

4.1.4 DBSCAN密度聚类算法优缺点

• 优点：
  1）可以对任意形状的稠密数据集进行聚类，相对的，K-Means之类的聚类算法一般只适用于凸数据集；
  2）可以在聚类的同时发现异常点，对数据集中的异常点不敏感；
  3）聚类结果没有偏倚，相对的，K-Means之类的聚类算法初始值对聚类结果有很大影响。
• 缺点：
  1）如果样本集的密度不均匀、聚类间距差相差很大时，聚类质量较差，这时用DBSCAN聚类一般不适合；
  2）如果样本集较大时，聚类收敛时间较长，此时可以对搜索最近邻时建立的KD树或者球树进行规模限制来改进；
  3）调参相对于传统的K-Means之类的聚类算法稍复杂，主要需要对距离阈值ϵ，邻域样本数阈值MinPts联合调参，不同的参数组合对最后的聚类效果有较大影响。

4.2 代码实现

import numpy as np
import pandas as pd
from pandas import DataFrame
from sklearn import datasets
import matplotlib.pyplot as plt

X, y = datasets.make_moons(500, noise=0.1, random_state=1)
df = pd.DataFrame(X, columns=['feature1', 'feature2'])
df.plot.scatter('feature1', 'feature2', title='row data', s=100)
plt.show()

'''
核心参数：
eps: float，ϵ-邻域的距离阈值
min_samples ：int，样本点要成为核心对象所需要的 ϵ-邻域的样本数阈值
'''
from sklearn.cluster import dbscan
core_samples, cluster_ids = dbscan(X, eps=0.2, min_samples=20)
df= pd.DataFrame(np.c_[X, cluster_ids], columns=['feature1', 'feature2', 'cluster_id'])
df['cluster_id'] = df['cluster_id'].astype('int')
df.plot.scatter('feature1', 'feature2', s=100,
                c=list(df['cluster_id']), cmap='rainbow', colorbar=False,
                alpha=0.6, title='DBSCAN cluster result')
plt.show()
print(df.head())

https://www.jianshu.com/p/e594c2ce0ac0

五、聚类算法评估—理论

5.1 估计聚类趋势

• 检测数据分布中是否存在非随机的簇结构，如果数据基本是随机的，那么聚类结构毫无意义；
  1）观察聚类误差是否随聚类类别数量的增加而单调变化，若数据是随机的，那么聚类误差随聚类类别数量增加而变化幅度较不显著，并且找不到一个合适的$K$对应数据的真实簇数。；
  2）霍普金斯统计量（Hopkins Statistic）
  从样本中随机找出$n$个点，记为$P_1,P_2,...,P_n$，对其中的每一个点$P_i$，在样本空间中找到一个离它最近的点并估计他们之间的距离$x_i$，从而得到距离$x_1,x_2,...,x_n$，然后从样本可能取值范围内随机生成$n$个点，记为$q_1,q_2,...,q_n$，对每个随机生成的点，找到一个离它最近的样本点并计算他们之间的距离$y_1,y_2,...,y_n$，则统计量H为：
  $$H=\frac{{\sum }_{i=1}^n{y}_i}{{\sum }_{i=1}^n{x}_i+{\sum }_{i=1}^n{y}_i}\text{\hspace{1em}(16)}$$
  若样本是随机分布的，则${\sum }_{i=1}^n{x}_i$和${\sum }_{i=1}^n{y}_i$接近，H接近于0.5，若聚类趋势明显，则${\sum }_{i=1}^n{y}_i>>{\sum }_{i=1}^n{x}_i$，随机生成的样本点距离应远大于实际样本点距离，H接近于1。

5.2 判定聚类簇数

• 确定聚类趋势后，可通过手肘法和Gap Statistic统计量确定簇数。

5.3 确定聚类质量

• 在无监督情况下，通过观察簇的分离情况和簇的紧凑情况来评估聚类的质量。
  1）轮廓系数
  给定一个点$p$，该点轮廓系数为：
  S ( p ) = b ( p ) − a ( p ) max ⁡ { a ( p ) , b ( p ) } (17) S(p)=\frac{b(p)-a(p)}{\max\{a(p), b(p)\}}\tag{17} S(p)=max{a(p),b(p)}b(p)−a(p)​(17)
  $a(p)$为$p$与同一簇中其他点$p^{\prime }$之间的平均距离，$a(p)$反映p所属簇中数据的紧凑程度；
  $b(p)$为$p$与另一簇中的点之间的最小的平均距离（在其他簇中，只计算与点p最近的一簇中的点与该点的平均距离），$b(p)$反映该簇与其他临近簇的分离程度；
  $b(p)$越大，$a(p)$越小，对应的簇聚类质量越好；
  因此将所有点对应的轮廓系数$S(p)$求平均来度量聚类结果的质量。
  2）均方根标准偏差（Root-mean-square standard deviation, RMSSTD）- - -衡量聚类的紧凑程度
  $$RMSSTD=(\frac{{\sum }_i{\sum }_{x\in C_i}||x-c_i|{|}^2}{P{\sum }_i(n_i-1)}{)}^{\frac{1}{2}}\text{\hspace{1em}(18)}$$
  其中$C_i$（大写）为第$i$个簇数，$c_i$（小写）为簇中心，$x\in C_i$为第$i$个簇的样本点，$n_i$为第$i$个簇的样本数量，$P$为样本维数，分母对点的维度P做了惩罚，维度越高，则整体的平方距离值越大。
  ${\sum }_i(n_i-1)=n-k$，即样本总数-聚类簇数，通常$k<<n$，因此${\sum }_i(n_i-1)$接近于常数。
  综上，RMSSTD可看做是经过归一化的标准差。
  3）R方（R-Ssquare）- - -衡量聚类的差异度
  $$RS=\frac{{\sum }_{x\in D}||x-c|{|}^2-{\sum }_i{\sum }_{x\in C_i}||x-c_i|{|}^2}{{\sum }_{x\in D}||x-c|{|}^2}\text{\hspace{1em}(19)}$$
  其中$D$为整个数据集，$c$为整个数据集的中心点，RS代表聚类后与聚类前相比，对应平方误差的改进程度。
  4）改进的Hubert $\Gamma$统计量- - -通过数据对的不一致性来评估聚类差异
  $$\Gamma =\frac{2}{n(n-1)}\sum _{x\in D}\sum _{y\in D}d(x,y)d_{x\in C_i,y\in C_j}(c_i,c_j)\text{\hspace{1em}(20)}$$
  其中$\frac{n(n-1)}{2}$为所有点对的数量，$d(x,y)$为点$x$到点$y$之间的距离，$d_{x\in c_i,y\in c_j}(c_i,c_j)$为$x$所在的类中心$c_i$与$y$所在的类中心$c_j$之间是的距离。
  理想情况下：对$d(x,y)$，若$d(x,y)$越小，对应的$d_{x\in c_i,y\in c_j}(c_i,c_j)$也越小，属于同一类时为0，反之越大。所以$\Gamma$值越大，说明聚类结果与样本的原始距离越吻合，聚类质量越高。

六、聚类算法评估—Python实践

6.1 同质性（Homogeneity）

• 如果所有的聚类都只包含属于单个类的成员的数据点，则聚类结果满足同质性。取值范围[0,1]，值越大意味着聚类结果与真实情况越符合。
• Python实现

from sklearn import metrics
print("Homogeneity值: %0.3f" % metrics.homogeneity_score(labels_true, labels_pred) 
# labels_true是真实的类标签、labels_pred是模型预测的类标签

6.2 完整性（Completeness）

• 如果作为给定类的成员的所有数据点是相同集群的元素，则聚类结果满足完整性。取值范围[0,1]，值越大意味着聚类结果与真实情况越符合。
• Python实现

from sklearn import metrics
print("Completeness值: %0.3f" % metrics.completeness_score(labels_true, labels_pred) 
# labels_true是真实的类标签、labels_pred是模型预测的类标签

6.3 同质性和完整性的调和平均（v_meansure_score）

• $V-measure=2\ast \frac{Homogeneity\ast Completeness}{Homogeneity+Completeness}$，取值范围[0,1]，值越大意味着聚类结果与真实情况越符合。
• Python实现

from sklearn import metrics
print("V-measure值: %0.3f" % metrics.v_measure_score(labels_true, labels_pred) 
# labels_true是真实的类标签、labels_pred是模型预测的类标签

6.4 轮廓系数

• 见式子（17）所示。
• Python实现

from sklearn import metrics
import numpy as np
from sklearn.cluster import KMeans
kmeans_model = KMeans(n_clusters=3, random_state=1).fit(X)
labels = kmeans_model.labels_
print("Silhouette Coefficient值: %0.3f"
      % metrics.silhouette_score(X, labels, metric='euclidean', sample_size=1000))                                         

6.5 兰德指数

• Python实现

from sklearn import metrics
print("Adjusted Rand-Index值: %.3f"
      % metrics.adjusted_rand_score(labels_true, labels_pred)) 
# labels_true是真实的类标签、labels_pred是模型预测的类标签               

6.6 互信息

• Python实现

from sklearn import metrics
print("Adjusted mutual info score值: %.3f"
      % metrics.adjusted_mutual_info_score(labels_true, labels_pred))
# labels_true是真实的类标签、labels_pred是模型预测的类标签