聚类方法

聚类是针对给定的样本，依据它们特征的相似度或距离，将其归并到若干个“类”或“簇”的数据分析问题。一个类是给定样本集合的一个子集。直观上，相似的样本聚集在相同的类，不相似的样本分散在不同的类。这里，样本之间的相似度或距离起着重要作用。

聚类的目的是通过得到的类或簇来发现数据的特点或对数据进行处理，在数据挖掘、模式识别等领域有着广泛的应用。聚类属于无监督学习，因为只是根据样本的相似度或距离将其进行归类，而类或簇事先并不知道。

聚类算法有很多，在这里我们介绍最常用的两种聚类算法：层次聚类（hierarchical clustering）和$k$均值聚类（$k$-means clustering）。层次聚类又有聚合（自下而上）和分裂（自上而下）两种方法。聚合法开始将每个样本各自分到一个类，之后将相距最近的两类合并，建立一个新的类，重复此操作直到满足停止条件，得到层次化的类别；分裂法开始将所有样本分到一个类，之后将已有类中相距最远的样本分到两个新的类，重复此操作直到满足停止条件，得到层次化的类别。$k$均值聚类是基于中心的聚类方法，通过迭代，将样本分到$k$个类中，使得每个样本与其所属的中心或均值最近；得到$k$个“平坦的”、“非层次化的”类别，构成对空间的划分。

一、聚类的基本概念

1. 相似度或距离

聚类的对象是观测数据，或样本集合。假设有$n$个样本，每个样本由$m$个属性的特征向量组成。样本集合可以用矩阵$X$表示

$$X={\left[\begin{array}{c}{x}_{ij}\end{array}\right]}_{m\times n}=\left[\begin{array}{cccc}{x}_{11}& x_{12}& \cdots & x_{1n}\\ x_{21}& x_{22}& \cdots & x_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ x_{m1}& x_{m2}& \cdots & x_{mn}\end{array}\right]\text{\hspace{1em}(1.1)}$$

矩阵的第$j$列表示第$j$个样本，$j=1,2,\cdots ,n$；第$i$行表示第$i$个属性，$i=1,2,\cdots ,m$；矩阵元素$x_{ij}$表示第$j$个样本的第$i$个属性值，$i=1,2,\cdots ,m$；$j=1,2,\cdots ,n$。

聚类的核心概念是相似度（similarity）或距离（distance），有多种相似度或距离的定义。因为相似度直接影响聚类的结果，所以其选择是聚类的根本问题。具体哪种相似度更合适取决于应用问题的特性。

（1）闵可夫斯基距离（Minkowski distance）

在聚类中，可以将样本集合看作是向量空间中点的集合，以该空间的距离表示样本之间的相似度。常用的距离有闵可夫斯基距离，特别是欧氏距离。林可夫斯基距离越大相似度越小，距离越小相似度越大。

定义1.1    给定样本集合$X$，$X$是$m$维实数向量空间${\mathbf{R}}^m$中点的集合，其中$x_i,x_j\in X$，$x_i={\left(\begin{array}{c}{x}_{1i},x_{2i},\cdots ,x_{mi}\end{array}\right)}^T$，$x_j={\left(\begin{array}{c}{x}_{1j},x_{2j},\cdots ,x_{mj}\end{array}\right)}^T$，样本$x_i$与样本$x_j$的闵可夫斯基距离（Minkowski distance）定义为

$$d_{ij}={\left(\sum _{k=1}^m|x_{ki}-x_{kj}{|}^p\right)}^{\frac{1}{p}}\text{\hspace{1em}(1.2)}$$

这里$p⩾1$。当$p=2$时称为欧氏距离（Euclidean distance），即

$$d_{ij}={\left(\sum _{k=1}^m|x_{ki}-x_{kj}{|}^2\right)}^{\frac{1}{2}}\text{\hspace{1em}(1.3)}$$

当$p=1$时称为曼哈顿距离（Manhattan distance），即

$$d_{ij}=\sum _{k=1}^m|x_{ki}-x_{kj}|\text{\hspace{1em}(1.4)}$$

当$p=\infty$时称为切比雪夫距离（Chebyshev distance），取各个坐标数值差的绝对值的最大值，即

$$d_{ij}=\underset{k}{\max }|x_{ki}-x_{kj}|\text{\hspace{1em}(1.5)}$$

（2）马哈拉诺比斯距离（Mahalanobis distance）

马哈拉诺比斯距离（Mahalanobis distance），简称马氏距离，也是另一种常用的相似度，考虑各个分量（特征）之间的相关性并与各个分量的尺度无关。马哈拉诺比斯距离越大相似度越小，距离越小相似度越大。

定义1.2    给定一个样本集合$X$，$X={\left[\begin{array}{c}{x}_{ij}\end{array}\right]}_{m\times n}$，其协方差矩阵记作$S$。样本$x_i$与样本$x_j$之间的马哈拉诺比斯距离$d_{ij}$定义为

$$d_{ij}={\left[{\left(x_i-x_j\right)}^T{S}^{-1}\left(x_i-x_j\right)\right]}^{\frac{1}{2}}\text{\hspace{1em}(1.6)}$$

其中

$$x_i={\left(\begin{array}{c}{x}_{1i},x_{2i},\cdots ,x_{mi}\end{array}\right)}^T,x_j={\left(\begin{array}{c}{x}_{1j},x_{2j},\cdots ,x_{mj}\end{array}\right)}^T\text{\hspace{1em}(1.7)}$$

当$S$为单位矩阵时，即样本数据的各个分量互相独立且各个分量的方差为1时，由式$(1.6)$知马氏距离就是欧氏距离，所以马氏距离是欧氏距离的推广。

（3）相关系数（correlation coefficient）

样本之间的相似度也可以用相关系数（correlation coefficient）来表示。相关系数的绝对值越接近于1，表示样本越相似；越接近于0，表示样本越不相似。

定义1.3    样本$x_i$与样本$x_j$之间的相关系数定义为

$$r_{ij}=\frac{\sum _{k=1}^m\left(x_{ki}-{\bar{x}}_i\right)\left(x_{kj}-{\bar{x}}_j\right)}{{\left[\sum _{k=1}^m{\left(x_{ki}-{\bar{x}}_i\right)}^2\sum _{k=1}^m{\left(x_{kj}-{\bar{x}}_j\right)}^2\right]}^{\frac{1}{2}}}\text{\hspace{1em}(1.8)}$$

其中

$${\bar{x}}_i=\frac{1}{n}\sum _{k=1}^n{x}_{ki},{\bar{x}}_j=\frac{1}{n}\sum _{k=1}^n{x}_{kj}$$

（4）夹角余弦（cosine）

样本之间的相似度也可以用夹角余弦（cosine）来表示。夹角余弦越接近于1，表示样本越相似；越接近于0，表示样本越不相似。

定义1.4    样本$x_i$与样本$x_j$之间的夹角余弦定义为

$$s_{ij}=\frac{\sum _{k=1}^m{x}_{ki}{x}_{kj}}{{\left[\sum _{k=1}^m{x}_{ki}^2\sum _{k=1}^m{x}_{kj}^2\right]}^{\frac{1}{2}}}\text{\hspace{1em}(1.9)}$$

由上述定义看出，用距离度量相似度时，距离越小样本越相似；用相关系数时，相关系数越大样本越相似。注意不同相似度度量得到的结果并不一定一致。例如下图：

图1.1  距离与相关系数的关系

从图上可以看出，如果从距离的角度看，$A$和$B$比$A$和$C$更相似；但从相关系数的角度看，$A$和$C$比$A$和$B$更相似。所以，进行聚类时，选择适合的距离或相似度非常重要。

2. 类或簇

通过聚类得到的类或簇，本质是样本的子集。如果一个聚类方法假定一个样本只能属于一个类，或类的交集为空集，那么该方法称为硬聚类（hard clustering）方法。否则，如果一个样本可以属于多个类，或类的交集不为空集，那么该方法称为软聚类（soft clustering）方法。在这里我们只考虑硬聚类方法。

用$G$表示类或簇（cluster），用$x_i$，$x_j$表示类中的样本，用$n_G$表示$G$中样本的个数，用$d_{ij}$表示样本$x_i$与样本$x_j$之间的距离。类或簇有多种定义，下面给出几个常见的定义。

定义1.5    设$T$为给定的正数，若集合$G$中任意两个样本$x_i$，$x_j$，有

$$d_{ij}⩽T$$

则称$G$为一个类或簇。

定义1.6    设$T$为给定的正数，若对集合$G$的任意样本$x_i$，一定存在$G$中的另一个样本$x_j$，使得

$$d_{ij}⩽T$$

则称$G$为一个类或簇。

定义1.7    设$T$为给定的正数，若对集合$G$中任意一个样本$x_i$，$G$中的另一个样本$x_j$满足

$$\frac{1}{n_G-1}\sum _{x_j\in G}{d}_{ij}⩽T$$

其中$n_G$为$G$中样本的个数，则称$G$为一个类或簇。

定义1.8    设$V$和$T$为给定的两个正数，如果集合$G$中任意两个样本$x_i$，$x_j$的距离$d_{ij}$满足

$$\frac{1}{n_G(n_G-1)}\sum _{x_i\in G}\sum _{x_j\in G}{d}_{ij}⩽T$$

$$d_{ij}⩽V$$

则称$G$为一个类或簇。

以上四个定义，第一个定义最常用，并且由它可推出其他三个定义。

类的特征可以通过不同角度来刻画，常用的特征有下面三种：

（1）类的均值${\bar{x}}_G$，又称为类的中心

$${\bar{x}}_G=\frac{1}{n_G}\sum _{i=1}^{n_G}{x}_i⩽T\text{\hspace{1em}(1.10)}$$

式中$n_G$是类$G$的样本个数。

（2）类的直径（diameter）$D_G$

类的直径$D_G$是类中任意两个样本之间的最大距离，即

$$D_G=\underset{x_i,x_j\in G}{\max }{d}_{ij}\text{\hspace{1em}(1.11)}$$

（3）类的样本散步矩阵（scatter matrix）$A_G$与样本协方差矩阵（covariance matrix）$S_G$

类的样本散布矩阵$A_G$为

$$A_G=\sum _{i=1}^{n_G}(x_i-{\bar{x}}_G)(x_i-{\bar{x}}_G{)}^T\text{\hspace{1em}(1.12)}$$

样本协方差矩阵$S_G$为

$$\begin{array}{rl}{S}_G& =\frac{1}{m-1}{A}_G\\ & =\frac{1}{m-1}\sum _{i=1}^{n_G}(x_i-{\bar{x}}_G)(x_i-{\bar{x}}_G{)}^T\end{array}\text{\hspace{1em}(1.13)}$$

其中$m$为样本的维数（样本属性的个数）。

3. 类与类之间的距离

下面考虑类$G_p$与类$G_q$之间的距离$D(p,q)$，也称为连接（linkage）。类与类之间的距离也有多种定义。

设类$G_p$包含$n_p$个样本，$G_q$包含$n_q$个样本，分别用${\bar{x}}_p$和${\bar{x}}_q$表示$G_p$和$G_q$的均值，即类的中心。

（1）最短距离或単连接（single linkage）

定义类$G_p$的样本与$G_q$的样本之间的最短距离为两类之间的距离

D p q = min ⁡ { d i j ∣ x i ∈ G p , x j ∈ G q } (1.14) D_{pq} = \min \left\{ d_{ij} \mid x_i \in G_p, x_j \in G_q \right\} \tag{1.14} Dpq​=min{dij​∣xi​∈Gp​,xj​∈Gq​}(1.14)

（2）最长距离或完全连接（complete linkage）

定义类$G_p$的样本与$G_q$的样本之间的最长距离为两类之间的距离

D p q = max ⁡ { d i j ∣ x i ∈ G p , x j ∈ G q } (1.15) D_{pq} = \max \left\{ d_{ij} \mid x_i \in G_p, x_j \in G_q \right\} \tag{1.15} Dpq​=max{dij​∣xi​∈Gp​,xj​∈Gq​}(1.15)

（3）中心距离

定义类$G_p$与类$G_q$的中心${\bar{x}}_p$与${\bar{x}}_q$之间的距离为两类之间的距离

$$D_{pq}=d_{{\bar{x}}_p{\bar{x}}_q}\text{\hspace{1em}(1.16)}$$

（4）平均距离

定义类$G_p$与类$G_q$任意两个样本之间距离的平均值为两类之间的距离

$$D_{pq}=\frac{1}{n_p{n}_q}\sum _{x_i\in G}\sum _{x_j\in G}{d}_{ij}\text{\hspace{1em}(1.17)}$$

三、$k$均值聚类

$k$均值聚类是基于样本集合划分的聚类算法。$k$均值聚类将样本集合划分为$k$个子集，构成$k$个类，将$n$个样本分到$k$个类中，每个样本到其所属类的中心的距离最小。每个样本只能属于一个类，所以$k$均值聚类是硬聚类。下面分别介绍$k$均值聚类的模型、策略、算法，讨论算法的特性及相关问题。

1. 模型

给定$n$个样本的集合 X = { x 1 , x 2 , ⋯   , x n } X = \left\{ x_1, x_2, \cdots, x_n \right\} X={x1​,x2​,⋯,xn​}，每个样本由一个特征向量表示，特征向量的维数是$m$。$k$均值聚类的目标是将$n$个样本分到$k$个不同的类或簇中，这里假设$k<n$。$k$个类$G_1,G_2,\cdots ,G_k$形成对样本集合$X$的划分，其中$G_i\cap G_j=\varnothing$，$\bigcup _{i=1}^k{G}_i=X$。用$C$表示划分，一个划分对应着一个聚类结果。

划分$C$是一个多对一的函数。事实上，如果把每个样本用一个整数 i ∈ { 1 , 2 , ⋯   , n } i \in \{ 1, 2, \cdots, n \} i∈{1,2,⋯,n}表示，每个类也用一个整数 l ∈ { 1 , 2 , ⋯   , k } l \in \{ 1, 2, \cdots, k \} l∈{1,2,⋯,k}表示，那么划分或聚类可以用函数$l=C(i)$表示，其中 i ∈ { 1 , 2 , ⋯   , n } i \in \{ 1, 2, \cdots, n \} i∈{1,2,⋯,n}， l ∈ { 1 , 2 , ⋯   , k } l \in \{ 1, 2, \cdots, k \} l∈{1,2,⋯,k}。所以$k$均值聚类的模型是一个从样本到类的函数。

2. 策略

$k$均值聚类归结为样本集合$X$的划分，或者从样本到类的函数的选择问题。$k$均值聚类的策略是通过损失函数的最小化选取最优的划分或函数$C^{\ast }$。

首先，采用欧氏距离平方（squared Euclidean distance）作为样本之间的距离$d(x_i,x_j)$

$$\begin{array}{rl}d(x_i,x_j)& =\sum _{k=1}^m{\left(x_{ki}-x_{kj}\right)}^2\\ & =\Vert x_i-x_j{\Vert }^2\end{array}\text{\hspace{1em}(3.1)}$$

然后，定义样本与其所属类的中心之间的距离的总和为损失函数，即

$$W(C)=\sum _{l=1}^k\sum _{C(i)=l}\Vert x_i-{\bar{x}}_l{\Vert }^2\text{\hspace{1em}(3.2)}$$

式中${\bar{x}}_l={\left({\bar{x}}_{1l},{\bar{x}}_{2l},\cdots ,{\bar{x}}_{ml}\right)}^T$是第$l$个类的均值或中心，$n_l=\sum _{i=1}^{n}I(C(i)=l)$，$I(C(i)=l)$是指数函数，取值为1或0。函数$W(C)$也称为能量，表示相同类中的样本相似的程度。

$k$均值聚类就是求解最优化问题：

$$\begin{array}{rl}{C}^{\ast }& =\arg \underset{C}{\min }W(C)\\ & =\arg \underset{C}{\min }\sum _{l=1}^k\sum _{C(i)=l}\Vert x_i-{\bar{x}}_l{\Vert }^2\end{array}\text{\hspace{1em}(3.3)}$$

相似的样本被聚到同类时，损失函数值最小，这个目标函数的最优化能达到聚类的目的。但是，这是一个组合优化问题，$n$个样本分到$k$类，所有可能分法的数目是第二类斯特林数：

$$S(n,k)=\frac{1}{k!}\sum _{l=1}^k(-1{)}^{k-l}(\genfrac{}{}{0ex}{}{k}{l})(k-l{)}^n\text{\hspace{1em}(3.4)}$$

第二类斯特林（Stirling）数    将$n$个不同的球放入$m$个无差别的盒子中，要求盒子非空，可以得到其方案数公式为

$$S(n,m)=\frac{1}{m!}\sum _{i=0}^m(-1{)}^{m-i}(\genfrac{}{}{0ex}{}{m}{i})(m-i{)}^n$$

证明    采用容斥原理的方法来证明。
先把盒子区分开，即认为给盒子加上标号顺序，最后答案再除以$m!$即可。
既然盒子有区别，那么先不考虑空盒，对于每个元素有$m$种选择，方法数是：$m^n$。但是这里显然包含了空盒的情况，即空盒数量是大于等于0的。那么不妨设在这$m^n$种方法中空盒数正好等于$h$的方法是$T(h)$，于是显然有

$$m!\cdot S(n,m)=T(0)\text{\hspace{1em}(3.5)}$$

$$T(m)=0\text{\hspace{1em}(3.6)}$$

以及

$$\sum _{h=0}^{m}T(h)=m^n\text{\hspace{1em}(3.7)}$$

那么现在我们考虑空盒子数量大于等于$k$的种类个数，我们先将这$k$个空盒取出来，有$C_m^k$种取法，那么这$n$个球每个球可以随意选择剩下的$(m-k)$个盒子中的一个，所以一共有$C_m^k(m-k{)}^n$种情况。
但是这里面有重复的情况，比如说空盒子数量正好等于$h$时$(h⩾k)$，不妨设这$h$个空盒子的编号为$1,2,\cdots ,h$，$C_m^k(m-k{)}^n$里面可以是先选出了编号为$1,2,\cdots ,k$这些盒子作为空盒子，然后在$(m-k{)}^n$里面选择整好将球没有放在编号为$k+1,k+2,\cdots ,h$的盒子中；也可以是先选出了编号为$1,2,\cdots ,k-1,k+1$这些盒子作为空盒子，然后在$(m-k{)}^n$里面选择整好将球没有放在编号为$k,k+2,k+3,\cdots ,h$的盒子中，这样便产生了完全相同的两种放法，而且都被计算在了$C_m^k(m-k{)}^n$里面，这表明，$C_m^k(m-k{)}^n$里面有重复的方案。
那么我们考虑在$C_m^k(m-k{)}^n$里面空盒子数量正好等于$h$时$(h⩾k)$重复的次数，我们发现，这里产生重复的原因就在于先选出来的$k$个空盒子不同，但是最后总共$h$个空盒子的编号是相同的。由于先选哪些都是等价的，所以重复的次数就是这$h$个空盒子中可以先选出哪$k$个空盒子的方案数，即为$C_h^k$，也就是说本来应该只算一次的空盒子数量正好等于$h$的这$T(h)$种方案在$C_m^k(m-k{)}^n$里面被算了$C_h^k$次。
所以有

$$C_m^k(m-k{)}^n=\sum _{h=k}^m{C}_h^{k}T(h)\text{\hspace{1em}(3.8)}$$

注意到对任意$k\in {\mathbb{N}}^{+}$，有

$$0=(1-1{)}^k=C_k^0-C_k^1+C_k^2+\cdots +(-1{)}^k{C}_k^k$$

所以有

$$\begin{array}{rl}& \sum _{i=0}^m(-1{)}^{m-i}(\genfrac{}{}{0ex}{}{m}{i})(m-i{)}^n\\ =& \sum _{i=0}^m\left[(-1{)}^{m-i}\sum _{h=i}^m{C}_h^{i}T(h)\right]\\ =& \sum _{i=0}^m\sum _{h=i}^m(-1{)}^{m-i}{C}_h^{i}T(h)\\ =& \sum _{h=0}^m\sum _{i=0}^h(-1{)}^i{C}_h^{i}T(h)\\ =& \sum _{h=0}^{m}T(h)\sum _{i=0}^h(-1{)}^i{C}_h^i\\ =& T(0)+\sum _{h=1}^{m}T(h)\sum _{i=0}^h(-1{)}^i{C}_h^i\\ =& T(0)+0\\ =& T(0)\end{array}$$

那么由公式$(3.5)$可得

$$S(n,m)=\frac{1}{m!}\sum _{i=0}^m(-1{)}^{m-i}(\genfrac{}{}{0ex}{}{m}{i})(m-i{)}^n$$

那么我们注意到，这个数是指数级的。事实上，$k$均值聚类的最优解求解问题是NP困难问题。所以现实中采用迭代的方法求解。

3. 算法

$k$均值聚类的算法是一个迭代的过程，每次迭代包括两个步骤。首先选择$k$个类的中心，将样本逐个指派到与其最近的中心的类中，得到一个聚类结果；然后更新每个类的样本均值中心，作为类的新的中心；重复以上步骤，直到收敛为止。具体过程如下。

首先，对于给定的中心值$(m_1,m_2,\cdots ,m_k)$，求一个划分$C$，使得目标函数极小化：

$$\underset{C}{\min }\sum _{l=1}^k\sum _{C(i)=l}\Vert x_i-m_l{\Vert }^2\text{\hspace{1em}(3.9)}$$

就是说在类中心确定的情况下，将每个样本分到 一个类中，使样本和其所属类的中心之间的距离总和最小。求解结果，将每个样本指派到与其最近的中心$m_l$的类$G_l$中。

然后，对给定的划分$C$，再求各个类的中心$(m_1,m_2,\cdots ,m_k)$，使得目标函数极小化：

$$\underset{m_1,\cdots ,m_k}{\min }\sum _{l=1}^k\sum _{C(i)=l}\Vert x_i-m_l{\Vert }^2$$

就是说在划分确定的情况下，使样本和其所属类的中心之间的距离总和最小。求解结果，对于每个包含$n_l$个样本的类$G_l$，更新其均值$m_l$：

$$m_l=\frac{1}{n_l}\sum _{C(i)=l}{x}_i,l=1,\cdots ,k$$

重复以上两个步骤，直到划分不再改变，得到聚类结果。如此，我们便可以得到$k$均值聚类算法。

算法2（$k$均值聚类算法）

输入：$n$个样本的集合$X$；

输出：样本集合的聚类$C^{\ast }$。

（1）初始化。令$t=0$，随机选择$k$个样本点作为初始聚类中心$m^{(0)}=(m_1^{(0)},m_2^{(0)},\cdots ,m_k^{(0)})$。

（2）对样本进行聚类。对固定的类中心$m^{(t)}=(m_1^{(t)},m_2^{(t)},\cdots ,m_k^{(t)})$，将每个样本指派到与其最近的中心的类中，构成聚类结果$C^{(t)}$。

（3）计算新的类中心。对聚类结果$C^{(t)}$，计算当前各个类中的样本的均值，作为新的类中心$m^{(t+1)}=(m_1^{(t+1)},m_2^{(t+1)},\cdots ,m_k^{(t+1)})$。

（4）如果迭代收敛或符合停止条件，输出$C^{\ast }=C^{(t)}$。
否则，令$t=t+1$，返回步骤（2）。

那么从上面我们可以看出$k$均值聚类算法的复杂度是$O(mnk)$，其中$m$是样本维数，$n$是样本个数，$k$是类别个数。

四、代码

本代码实现了高斯混合模型参数估计的聚类算法，选取了三个二维高斯分布模型${\mathbf{\mu }}_1=\left(\begin{array}{c}2\\ 1\end{array}\right)$，${\mathbf{\sigma }}_1^2=\left(\begin{array}{cc}0.5& 0\\ 0& 0.2\end{array}\right)$；${\mathbf{\mu }}_2=\left(\begin{array}{c}5\\ 9\end{array}\right)$，${\mathbf{\sigma }}_2^2=\left(\begin{array}{cc}0.3& 0\\ 0& 0.3\end{array}\right)$；${\mathbf{\mu }}_3=\left(\begin{array}{c}10\\ 6\end{array}\right)$，${\mathbf{\sigma }}_3^2=\left(\begin{array}{cc}0.5& 0\\ 0& 0.3\end{array}\right)$，每个模型都随机生成100个点，然后用聚类算法对这300个点进行划分和聚类。

2. $k$均值聚类

初值设定为

$$C^{(0)}=\left(\begin{array}{cc}3.0& 2.0\\ 6.0& 8.0\\ 8.0& 7.0\end{array}\right)$$

进行迭代，最终得到结果如下：

图4.1  k-means算法迭代结果 图4.1  k-means算法对高斯混合模型的实现

完整代码如下：

import numpy as np
import matplotlib.pyplot as plt


def create_data(mean, cov, size):
    f = np.random.multivariate_normal(mean, cov, size=size)
    X = f.T
    return X


# 寻找点a现在应该在哪个类里面
def find(x, c):
    index = 0
    d1 = np.linalg.norm(x - c[0], ord=2)
    d2 = np.linalg.norm(x - c[1], ord=2)
    d3 = np.linalg.norm(x - c[2], ord=2)
    if d1 > d2:
        index = 1
    d = min(d1, d2)
    if d > d3:
        index = 2
    return index


# k-means聚类算法
def k_means(X, c):
    x = []
    lens = len(X)
    for i in range(lens):
        tmp_x = X[i].T
        num = len(tmp_x)
        for j in range(num):
            x.append(tmp_x[j])
    x = np.array(x)
    nums = len(x)
    tmp_c = np.array([[0.0, 0.0], [0.0, 0.0], [0.0, 0.0]])
    while True:
        c_index = np.array([0, 0, 0])
        for i in range(nums):
            index = find(x[i], c)
            tmp_c[index] += x[i]
            c_index[index] += 1
        for i in range(3):
            tmp_c[i, 0] /= c_index[i]
            tmp_c[i, 1] /= c_index[i]
        if (c == tmp_c).all():
            break
        else:
            c = tmp_c
            tmp_c = np.array([[0.0, 0.0], [0.0, 0.0], [0.0, 0.0]])
    return c


if __name__ == '__main__':
    size = 100
    mean1 = np.array([2, 1])
    cov1 = np.array([[0.5, 0], [0, 0.2]])
    mean2 = np.array([5, 9])
    cov2 = np.array([[0.3, 0], [0, 0.3]])
    mean3 = np.array([10, 6])
    cov3 = np.array([[0.5, 0], [0, 0.3]])

    X1 = create_data(mean1, cov1, size)
    X2 = create_data(mean2, cov2, size)
    X3 = create_data(mean3, cov3, size)
    X = np.array([X1, X2, X3])

    c = np.array([[3.0, 2.0], [6.0, 8.0], [8.0, 7.0]])
    c = k_means(X, c)
    print("C = {}".format(c))

    plt.scatter(X1[0], X1[1], c='r', marker='+')
    plt.scatter(X2[0], X2[1], c='g', marker='+')
    plt.scatter(X3[0], X3[1], c='b', marker='+')

    plt.scatter(c[0, 0], c[0, 1], s=200, c='black', marker='*')
    plt.scatter(c[1, 0], c[1, 1], s= 200, c='black', marker='*')
    plt.scatter(c[2, 0], c[2, 1], s=200, c='black', marker='*')

    plt.xlabel("X1")
    plt.ylabel("X2")
    plt.show()