EM应用—K均值聚类与混合高斯分布

期望最大化（expectation-maximization，E-M）是一种非常强大的算法，应用于数据科学的很多场景中。k-means 是EM算法的一个非常简单并且易于理解的应用，本文先从K均值聚类讲起，然后引出K均值的缺陷，提出了混合高斯模型对K均值聚类进行优化。混合高斯模型也使用了EM算法。

关于EM算法可以参考文章：EM算法详解

K均值聚类

K均值聚类将样本集合划分到$k$ 个子集，构成$k$ 个类，将$n$ 个样本分到$k$ 个类中，使得每个样本到其所属类的中心距离最小。每个样本只能属于一个类。

模型

我们先来约定以下变量的含义：

(1) X = { x ⃗ 1 , x ⃗ 2 , … x ⃗ n } X = \{\vec x_1,\vec x_2,\ldots\vec x_n\} X={x 1​,x 2​,…x n​} 表示$n$个样本的集合，每个样本${\vec{x}}_n$ 是一个维度为m的特征向量。$k$ 均值聚类的目标是将$n$ 个样本分到$k$ 个不同的类(簇)中，这里假设$k<n$。

(2) $G_1,G_2\dots G_k$ 表示$X$ 的$k$个类。满足：两个类$G_i,G_j$ 之间没有交集，所有的样本都在这$k$ 类里。

聚类模型从函数的角度来说，是指存在一个一对一的函数$C$，使得：
$$l=C({\vec{x}}_i)$$
其中： l ∈ { 1 , 2 … k } l \in \{1,2\ldots k\} l∈{1,2…k} 表示每一个类。上式表示样本${\vec{x}}_i$ 属于第$l$个类别。

我们接下来来说说如何确定这个$C$ 来获得最小的损失函数，从而实现最优划分。

策略

（1）损失函数

确定一个模型之前我们先定义模型的损失函数。

首先，采用欧氏距离平方作为样本之间的距离 $d({\vec{x}}_i,{\vec{x}}_j)$
$$d({\vec{x}}_i,{\vec{x}}_j)=||{\vec{x}}_i-{\vec{x}}_j|{|}^2=\sum _{k=1}^m(x_{ki}-x_{kj}{)}^2$$
其次，定义样本与其所属类的中心之间的距离的总和为损失函数，即
$$W(C)=\sum _{l=1}^k\sum _{C(i)=l}||{\vec{x}}_i-{\vec{x}}_l|{|}^2$$
其中： x ⃗ l = { x ⃗ 1 l , x ⃗ 2 l , … x ⃗ m l } \vec x_l = \{\vec x_{1l},\vec x_{2l},\ldots\vec x_{ml}\} x l​={x 1l​,x 2l​,…x ml​} 是第$l$ 类的均值中心。$n_l={\sum }_{i=1}^{n}I(C(i)=l)$，$I(C(i)=l)$是指示函数，第i个样本属于类别$l$时取1，否则取0。W© 衡量了相同类别的样本的相似程度。函数$W(C)$ 也称为能量，表示相同类中的样本相似程度。

k 均值聚类问题就是求解最优化问题：
$$C^{\ast }=\arg \underset{C}{\min }W(C)=\arg \underset{C}{\min }\sum _{l=1}^k\sum _{C(i)=l}||{\vec{x}}_i-{\vec{x}}_l|{|}^2$$
相似的样本被聚到同类时，损失函数最小。

算法流程

k 均值聚类是一个迭代的过程，每次迭代包括两个步骤：

（1）选择k个类的中心$(m_1,m_2\dots ,m_k)$，将样本逐个指派到与其最近的中心类中得到一个暂时的聚类结果，这样使得样本和其所属类中心的距离最小，也即最小化损失函数：
$$\underset{C}{\min }\sum _{l=1}^k\sum _{C(i)=l}||{\vec{x}}_i-{\vec{m}}_l|{|}^2$$

（2）更新每个类的样本均值，作为新的中心。
$$m_l=\frac{1}{n_l}\sum _{C(i)=l}{\vec{x}}_i,l=1,2\dots k$$
重复以上两个步骤，直到收敛。

算法特性

（1）以欧氏距离的平方表示样本之间的距离。以样本和其所属类别的中心距离的总和为最优化的目标函数。

（2）初始中心的选择直接影响聚类结果

（3）不能保证得到全局最优解

（4）类别数$k$ 事先指定，而实际应用中最优的$k$ 值时不知道的。需要不断的测试。类别数与类平均直径的关系如下：

一般情况下类别数较小时，类平均直径较大，类别数增大超过某个值以后，平均直径趋于不变。这个值对应的类别数就是最优值。

K均值与EM算法

K均值聚类的过程中蕴含着EM算法的思想。

简单来说，期望最大化方法包含以下步骤：

(1)猜测一些簇中心点

(2)重复直至收敛。

  a期望步骤(E-step) :将点分配至离其最近的簇中心点。

  b.最大化步骤(M-step) :将簇中心点设置为所有点坐标的平均值。

期望步骤(E-step) 不断更新每个点是属于哪一个簇的期望值，最大化步骤(M-step)计算关于簇中心点的拟合函数值最大化对应坐标，在本例中，通过简单地求每个簇中所有数据点坐标的平均值得到了簇中心点坐标。

K均值聚类的缺陷

如果你的数据长这样，K均值聚类可以实现较好的效果：

但是如果你的数据长这样：

用K均值聚类的结果可能是这样：

不是我们预期的聚类结果。我们想要的聚类效果应该是下面这样的：

这个例子说明: 依靠基于距离模型的K均值聚类不能解决这个问题，K均值聚类忽略了一个问题：即数据来自不同的分布，所以我们在聚类之前应该先考虑数据的分布，在满足数据分布相同的情况下再进行聚类。也就是使用基于分布的模型—高斯混合模型。

高斯混合模型(GMM)

高斯混合模型(Gaussian Mixture Model)模型中的隐变量是不知道数据来源于哪个模型，我们通过EM算法来解决这个问题。

GMM概率图模型

我们先通过GMM的概率图模型来了解一下符合混合高斯模型的数据是如何生成的。

GMM的概率图模型如下：

由概率图模型可知观测数据的产生过程如下：

（1）依概率${\alpha }_k$选择第$k$ 个高斯分布模型$\varphi (y|{\theta }_k)$

（2）依第$k$个分模型的概率分布$\varphi (x|{\theta }_k)$生成观测数据$x_j$

这个过程会使得数据来自不同的分布。

GMM模型

GMM的数学模型长这样:
$$P(y|\theta )=\sum _{k=1}^K{\alpha }_k\varphi (y|{\theta }_k)$$
其中，${\alpha }_k$是系数，${\alpha }_k\ge 0$，$\sum _{k=1}^K{\alpha }_k=1$, $\varphi (y|{\theta }_k)$ 是高斯分布密度，${\theta }_k=(\mu ,{\sigma }^2)$
$$\varphi (y|{\theta }_k)=\frac{1}{\sqrt{2\pi }{\sigma }_k}\exp \left(-\frac{(y-{\mu }_k{)}^2}{2{\sigma }_k^2}\right)$$
上式表示第k个分模型。

从数学模型我们可以看出：GMM是概率分布的加权求和。

由于数据来源于哪个模型，因此我们无法直接用极大似然估计求得模型得最优参数，接下来我们从EM得角度来求解。

GMM的EM算法

已知观测数据$y_1,y_2,\dots ,y_N$，由高斯混合模型生成
$$P(y|\theta )=\sum _{k=1}^K{\alpha }_k\varphi (y|{\theta }_k)$$
其中， $\theta =({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_K;{\theta }_1,{\theta }_2,\dots ,{\theta }_K)$

1、明确隐变量，写出完全数据的对数似然函数。

观测数据$y_j,j=1,2,\dots ,N$的产生过程在上述概率图模型中已经说过, 是已知的

但反映观测数据$y_j$来自第$k$个分模型的数据是未知的, $k=1,2,\dots ,K$ 以**隐变量${\gamma }_{jk}$**表示
注意这里${\gamma }_{jk}$的维度$(j\times k)$
γ j k = { 1 , 第 j 个 观 测 来 自 第 k 个 分 模 型 0 , 否 则 j = 1 , 2 , … , N ; k = 1 , 2 , … , K ; γ j k ∈ { 0 , 1 } \gamma_{jk}= \begin{cases} 1, &第j个观测来自第k个分模型\\ 0, &否则 \end{cases}\\ j=1,2,\dots,N; k=1,2,\dots,K; \gamma_{jk}\in\{0,1\} γjk​={1,0,​第j个观测来自第k个分模型否则​j=1,2,…,N;k=1,2,…,K;γjk​∈{0,1}
${\gamma }_{jk}$ 是一种one-hot 的形式，举例：假设观测数据$y_j$ 来源于第2个模型，则：
$${\gamma }_{j2}=(0,1,0\dots ,0)$$
即：除了第二个位置外，其他位置都为0。

有了观测数据$y_i$ 和未观测数据${\gamma }_{jk}$ ，完全数据如下：
$$(y_j,{\gamma }_{j1},{\gamma }_{j2},\dots ,{\gamma }_{jK},k=1,2,\dots ,N)$$
所以，完全数据的似然函数：
$$\begin{array}{rl}P(y,\gamma |\theta )=& \prod _{j=1}^{N}P(y_j,{\gamma }_{j1},{\gamma }_{j2},\dots ,{\gamma }_{jK}|\theta )\\ =& \prod _{j=1}^N\prod _{k=1}^K{\left[{\alpha }_k\varphi (y_j|{\theta }_k)\right]}^{{\gamma }_{jk}}\\ =& \prod _{k=1}^K{\alpha }_k^{n_k}\prod _{j=1}^N{\left[\varphi (y_j|{\theta }_k)\right]}^{{\gamma }_{jk}}\\ =& \prod _{k=1}^K{\alpha }_k^{n_k}\prod _{j=1}^N{\left[\frac{1}{\sqrt{2\pi }{\sigma }_k}\exp \left(-\frac{(y_j-{\mu }_k{)}^2}{2{\sigma }_k^2}\right)\right]}^{{\gamma }_{jk}}\end{array}$$
其中$n_k={\sum }_{j=1}^N{\gamma }_{jk},{\sum }_{k=1}^K{n}_k=N$

所以完全数据的对数似然函数如下：
$$\begin{array}{rl}\log P(y,\gamma |\theta )& =\prod _{k=1}^K{\alpha }_k^{n_k}\prod _{j=1}^N{\left[\frac{1}{\sqrt{2\pi }{\sigma }_k}\exp \left(-\frac{(y_j-{\mu }_k{)}^2}{2{\sigma }_k^2}\right)\right]}^{{\gamma }_{jk}}\\ & =\sum _{k=1}^K\left\{n_k\log {\alpha }_k+\sum _{j=1}^N{\gamma }_{jk}\left[\log \left(\frac{1}{\sqrt{2\pi }}\right)-\log {\sigma }_k-\frac{1}{2{\sigma }^2}(y_j-{\mu }_k{)}^2\right]\right\}\end{array}$$

我们假设：${\prod }_{j=1}^N{\left[\frac{1}{\sqrt{2\pi }{\sigma }_k}\exp \left(-\frac{(y_j-{\mu }_k{)}^2}{2{\sigma }_k^2}\right)\right]}^{{\gamma }_{jk}}={\omega }_k$

则:
$$\begin{array}{rl}\log P(y,\gamma |\theta )& =log\prod _{k=1}^K{\alpha }_k^{n_k}{\omega }_k\\ & =\sum _{k=1}^K(n_{k}log{\alpha }_k+log{\omega }_k)\end{array}$$
继续求$log{\omega }_k$
$$\begin{array}{rl}log{\omega }_k& =log\prod _{j=1}^N{\left[\frac{1}{\sqrt{2\pi }{\sigma }_k}\exp \left(-\frac{(y_j-{\mu }_k{)}^2}{2{\sigma }_k^2}\right)\right]}^{{\gamma }_{jk}}\\ & =\sum _{j=1}^N{\gamma }_{jk}(log\frac{1}{\sqrt{2\pi }}-log{\sigma }_k-\frac{(y_j-{\mu }_k{)}^2}{2{\sigma }_k^2})\end{array}$$
所以：
$$\begin{array}{r}\log P(y,\gamma |\theta )=\sum _{k=1}^K\left\{n_k\log {\alpha }_k+\sum _{j=1}^N{\gamma }_{jk}\left[\log \left(\frac{1}{\sqrt{2\pi }}\right)-\log {\sigma }_k-\frac{1}{2{\sigma }^2}(y_j-{\mu }_k{)}^2\right]\right\}\end{array}$$
2、EM算法的E步：确定Q函数

把$Q$ 函数表示成参数形式
$$\begin{array}{rl}Q(\theta ,{\theta }^{(i)})=& E[\log P(y,\gamma |\theta )|y,{\theta }^{(i)}]\\ =& E\left\{\sum _{k=1}^K\left\{n_k\log {\alpha }_k+\sum _{j=1}^N{\gamma }_{jk}\left[\log \left(\frac{1}{\sqrt{2\pi }}\right)-\log {\sigma }_k-\frac{1}{2{\sigma }^2(y_j-{\mu }_k{)}^2}\right]\right\}\right\}(1)\\ =& \sum _{k=1}^K\left\{E\left\{n_k\log {\alpha }_k+\sum _{j=1}^N{\gamma }_{jk}\left[\log \left(\frac{1}{\sqrt{2\pi }}\right)-\log {\sigma }_k-\frac{1}{2{\sigma }^2(y_j-{\mu }_k{)}^2}\right]\right\}\right\}(2)\\ =& \sum _{k=1}^K\left\{E\left\{\sum _{j=1}^N{\gamma }_{jk}\log {\alpha }_k+\sum _{j=1}^N{\gamma }_{jk}\left[\log \left(\frac{1}{\sqrt{2\pi }}\right)-\log {\sigma }_k-\frac{1}{2{\sigma }^2(y_j-{\mu }_k{)}^2}\right]\right\}\right\}(3)\\ =& \sum _{k=1}^K\left\{\sum _{j=1}^N(E{\gamma }_{jk})\log {\alpha }_k+\sum _{j=1}^N(E{\gamma }_{jk})\left[\log \left(\frac{1}{\sqrt{2\pi }}\right)-\log {\sigma }_k-\frac{1}{2{\sigma }^2(y_j-{\mu }_k{)}^2}\right]\right\}(4)\end{array}$$
说明：

(1)到(2)是因为期望的基本性质：和的期望等于期望的和。

(2)到(3)是因为$n_k={\sum }_{j=1}^N{\gamma }_{jk}$

(3)到(4)是因为在Q函数是指完全数据的对数似然函数对隐变量（此处是指数据来自于哪个模型未知）的条件概率分布$P({\gamma }_{jk}|Y,{\theta }^{(i)})$ 求期望，所以与隐变量${\gamma }_{jk}$无关的统统看作常数。

因此，我们只需要计算$E({\gamma }_{jk}|y_j,\theta )$， 记为${\hat{\gamma }}_{jk}$

$$\begin{array}{rl}{\hat{\gamma }}_{jk}=& E({\gamma }_{jk}|y_j,\theta )=1\times P({\gamma }_{jk}=1|y_j,\theta )(1)\\ =& \frac{P({\gamma }_{jk}=1,y_j|\theta )}{{\sum }_{k=1}^{K}P({\gamma }_{jk}=1,y_j|\theta )}(2)\\ =& \frac{P(y_j|{\gamma }_{jk}=1,\theta )P({\gamma }_{jk}=1|\theta )}{{\sum }_{k=1}^{K}P(y_j|{\gamma }_{jk}=1,\theta )P({\gamma }_{jk}=1|\theta )}\\ =& \frac{{\alpha }_k\varphi (y_j|{\theta }_k)}{{\sum }_{k=1}^K{\alpha }_k\varphi (y_j|{\theta }_k)}\end{array}$$
上式(1)到(2) 的实现细节：
$$\begin{array}{rl}P({\gamma }_{jk}=1|y_j,\theta )& =\frac{P({\gamma }_{jk}=1,y_j,\theta )}{P(y_j,\theta )}\\ & =\frac{P({\gamma }_{jk}=1,y_j|\theta )P(\theta )}{P(y_j|\theta )P(\theta )}\\ & =\frac{P({\gamma }_{jk}=1,y_j|\theta )}{P(y_j|\theta )}\\ & =\frac{P({\gamma }_{jk}=1,y_j|\theta )}{{\sum }_{k=1}^{K}P({\gamma }_{jk}=1,y_j|\theta )}\end{array}$$

我们定义${\hat{\gamma }}_{jk}$ 是指分模型 k 对观测数据 $y_j$ 的响应度：即当前模型参数下第$j$ 个观测数据来自第 $k$ 个模型的概率。

又因为：
$$\sum _{j=1}^{N}E{\gamma }_{jk}=\sum _{j=1}^{N}P({\gamma }_{jk}=1|y,\theta )=n_k$$
注：当${\gamma }_{jk}=1$ 时，$P({\gamma }_{jk}=1|y,\theta )=1$

将${\hat{\gamma }}_{jk}$ 和${\sum }_{j=1}^{N}E{\gamma }_{jk}$ 带入Q函数，得：
$$Q(\theta ,{\theta }^{(i)})=\sum _{k=1}^K{n}_k\log {\alpha }_k+\sum _{j=1}^N{\hat{\gamma }}_{jk}\left[\log \left(\frac{1}{\sqrt{2\pi }}\right)-\log {\sigma }_k-\frac{1}{2{\sigma }_k^2}(y_j-{\mu }_k{)}^2\right]$$
3、确定EM算法的M步

求函数$Q(\theta ,{\theta }^{(i)})$对$\theta$的极大值，分别求$\sigma ,\mu ,\alpha$
$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})$$
注：$\arg \max$ 就是求Q的极值对应的参数$\theta$，如说是离散的，遍历所有值，最大查找，如果是连续的，偏导为零求极值。

• $\frac{\partial Q}{\partial {\mu }_k}=0,\frac{\partial Q}{\partial {\sigma }^2}=0$ 得到${\hat{\mu }}_k,{\hat{\sigma }}_k^2$
  $$\begin{gathered} \widehat{{\mu }_k}=\frac{{\sum }_{j=1}^N\widehat{{\gamma }_{jk}}{y}_i}{{\sum }_{j=1}^N\widehat{{\gamma }_{jk}}} \\ \widehat{{\sigma }_k}=\frac{{\sum }_{j=1}^N\widehat{{\gamma }_{jk}}(y_i-{\mu }_k{)}^2}{{\sum }_{j=1}^N\widehat{{\gamma }_{jk}}} \end{gathered}$$

• ${\sum }_{k=1}^K{\alpha }_k=1,\frac{\partial Q}{\partial {\alpha }_k}=0$ 得到${\alpha }_k$
  $$\widehat{a_k}=\frac{{\sum }_{j=1}^N\widehat{{\gamma }_{jk}}}{N}$$

4、停止条件

重复以上计算，直到对数似然函数值不再有明显的变化为止。

这就是利用EM算法实现GMM模型的整个过程。