上一篇文章，我们讲的期望最大化(EM)算法是一种非常强大的算法，应用于数据科学的许多场景。k-means是该算法非常简单且易于理解的一个应用。

k-means聚类

k均值聚类算法(k-means)将样本集合划分为k个子集，也就是将n个样本划分到k个类别中，每个样本到类别的中心距离最近。

EM角度的理解

如果从EM角度来理解的话，k均值聚类算法的隐变量是聚类的中心，模型的参数是每个数据点属于哪个分类。k-means算法描述如下：

（1）选择初始聚类中心

（2）E步：将样本分配到距离最近的中心形成簇

（3）M步：将簇的中心设置为簇内各个样本的均值

（4）重复（2）（3）直到收敛

E步就是估计隐含类别y的期望值，M步调整其他参数使得在给定类别y的情况下，极大似然估计P(x,y)能够达到极大值。然后在其他参数确定的情况下，重新估计y，周而复始，直至收敛。

算法流程

k-means可以采用欧式距离的平方作为样本之间的距离度量：
$d(x_i, x_j) = \sum^m_{k=1}(x_{ki}-x_{kj})^2=||x_i-x_j||^2_2$d(xi​,xj​)=k=1∑m​(xki​−xkj​)2=∣∣xi​−xj​∣∣22​
然后定义样本与所属类的中心之间的距离的总和为损失函数：
$W(C) = \sum^k_{l=1}\sum_{C(i)=1}||x_i-\bar x_l||^2_2$W(C)=l=1∑k​C(i)=1∑​∣∣xi​−xˉl​∣∣22​
$\bar x_l=(\bar x_{1l},\bar x_{2l},...,\bar x_{ml})^T$xˉl​=(xˉ1l​,xˉ2l​,...,xˉml​)T是第$l$l个类的均值或者中心，$n_l=\sum_{i=1}^n I(C(i)=l)$nl​=∑i=1n​I(C(i)=l)，$I(C(i)=l)$I(C(i)=l)是指示函数，第i个样本属于类别$l$l时取1，否则取0。W© 衡量了相同类别的样本的相似程度。

k-means算法就是求出最佳的k，使得各个类别内的样本相似度尽可能的高：
$C^* = \arg\min\limits_{C}W(C) =\arg\min\limits_{C}\sum^k_{l=1}\sum_{C(i)=1}||x_i-\bar x_l||^2_2$C∗=argCmin​W(C)=argCmin​l=1∑k​C(i)=1∑​∣∣xi​−xˉl​∣∣22​
虽然上述目标函数的最优化可以达到聚类的目的，但是n个样本划分到k个类别，可能的划分方案的个数是指数级别的，难以求解。

为此，k均值聚类算法采用迭代的方式达到聚类的效果，每次迭代包含两个步骤：

1. 随机生成或者选择k个样本作为聚类中心，将每个样本分配到最近的中心，得到一个聚类结果。
2. 计算每个类内所有样本的均值，作为这个类别新的聚类中心。
3. 重复1，2直到聚类结果不再发生变化。

迭代过程如下面的动图所示：

具体过程如下：

首先，对于给定的k个中心$(m_1, m_2,...,m_k)$(m1​,m2​,...,mk​)，求一个划分C，每个样本划分到一个类中，使得样本和所属类的中心之间的距离总和最小，即下面的式子极小化：
$W(C)=\min\limits_C\sum_{l=1}^k\sum_{C(i)=l}||x_i-m_l||^2_2$W(C)=Cmin​l=1∑k​C(i)=l∑​∣∣xi​−ml​∣∣22​
然后，对给定的划分C或者说在类中心确定的情况下，再求各个类的中心$(m_1, m_2,...,m_k)$(m1​,m2​,...,mk​)，使得各个样本到这个中心的距离之和最小，一般要求的这个中心就是各个类别内样本的均值。这个过程描述如下：
$W(C)=\min\limits_{m_1,...,m_k}\sum_{l=1}^k\sum_{C(i)=l}||x_i-m_l||^2_2$W(C)=m1​,...,mk​min​l=1∑k​C(i)=l∑​∣∣xi​−ml​∣∣22​
对于包含了$n_l$nl​个样本的$C_l$Cl​，更新其均值$m_l$ml​。也就是说，具体求新的中心的过程如下：
$m_l = \frac{1}{n_l}\sum_{C(i)=l}x_i, l=1,2,...k$ml​=nl​1​C(i)=l∑​xi​,l=1,2,...k
重复上面的步骤，直到划分不再改变，就得到了聚类的结果。

特点

k均值聚类算法的总体特点：

1. 是基于划分的聚类方法
2. 类别k事先指定
3. 以欧式距离平方作为距离的度量方式
4. 以中心或者样本的均值表示类别
5. 以样本和所属的类的中心的距离之和作为目标函数
6. 算法是迭代算法，不能保证得到全局最优。

k均值聚类算法属于启发式方法，不能保证达到全局最优，受到初始中心的很大影响。

关于初始中心的选择，可以先用层次聚类对样本进行聚类，得到k个类别后，再从k个类别中选取一个与中心距离最近的点。

k值选择

关于类别数k的选择，在实际应用中需要尝试不同的k值，检验聚类结果的质量。聚类结果可以用类别的平均直径来衡量。

如上图所示，一般地，类别数量变大时，类别平均直径会变小，当类别数量超过某个值后，平均直径不再改变，这个值就是最优的k值。

局限性

不过，k-means聚类也有比较明显的局限性。当各个样本的分布接近圆形(球形)的时候，效果还不错：

但是，当各个样本的分布不那么“圆”的时候，聚类效果就打折扣了：

因此，我们需要一种不同的方法来为数据点分配聚类。因此，我们将不再使用基于距离的模型，而是使用基于分布的模型。效果如下：

下面要介绍的高斯混合模型就是基于分布的模型。

高斯混合模型

GMM的问题描述

高斯混合模型(GMM)假设存在一定数量的高斯分布，每个分布代表一个簇。因此，高斯混合模型倾向于将属于单一分布的数据点聚在一起。

高斯混合模型的数学公式:
$P(y|\theta)=\sum\limits^{K}_{k=1}\alpha_k\phi(y|\theta_k)$P(y∣θ)=k=1∑K​αk​ϕ(y∣θk​)
其中，$\alpha_k$αk​是系数，$\alpha_k\ge0$αk​≥0，$\sum\limits^{K}_{k=1}\alpha_k=1$k=1∑K​αk​=1, $\phi(y|\theta_k)$ϕ(y∣θk​) 是高斯分布密度，$\theta_k=(\mu_k,\sigma_k^2)$θk​=(μk​,σk2​)
$\phi(y|\theta_k)=\frac{1}{\sqrt{2\pi}\sigma_k}\exp\left(-\frac{(y-\mu_k)^2}{2\sigma_k^2}\right)$ϕ(y∣θk​)=2π​σk​1​exp(−2σk2​(y−μk​)2​)
上式表示第k个分模型。

GMM的数学公式在形式上可以看成是高斯分布概率密度函数(PDF)的加权求和：
$\sum(权重\times 分布密度)=分布$∑(权重×分布密度)=分布

书中描述的是一维的高斯混合模型，d维的形式如下[^2]，被称作多元正态分布，也叫多元高斯分布
$\phi(y|\theta_k)=\frac{1}{\sqrt{(2\pi)^d|\Sigma|}}\exp\left(-\frac{(y-\mu_k)^T\Sigma^{-1}(y-\mu_k)}{2}\right)$ϕ(y∣θk​)=(2π)d∣Σ∣​1​exp(−2(y−μk​)TΣ−1(y−μk​)​)
其中，$\Sigma\in \R^{n\times n}$Σ∈Rn×n是协方差矩阵。

已知观测数据$y_1, y_2, \dots , y_N$y1​,y2​,…,yN​，由高斯混合模型生成：
$P(y|\theta)=\sum_{k=1}^K\alpha_k\phi(y|\theta_k) =\prod_{j=1}^N\sum_{k=1}^K\alpha_k\phi(y_j|\theta_k)$P(y∣θ)=k=1∑K​αk​ϕ(y∣θk​)=j=1∏N​k=1∑K​αk​ϕ(yj​∣θk​)
其中， $\theta=(\alpha_1,\alpha_2,\dots,\alpha_K;\theta_1,\theta_2,\dots,\theta_K)$θ=(α1​,α2​,…,αK​;θ1​,θ2​,…,θK​)，我们要做的是使用EM算法估计GMM的参数$\theta$θ。

1，明确隐变量

观测数据$y_j, j=1,2,\dots,N$yj​,j=1,2,…,N的产生方式：

1. 依概率$\alpha_k$αk​选择第$k$k个高斯分布分模型$\phi(y|\theta_k)$ϕ(y∣θk​);
2. 依第$k$k个分模型的概率分布$\phi(y|\theta_k)$ϕ(y∣θk​)生成观测数据$y_j$yj​

但是，反映观测数据$y_j$yj​来自第$k$k个分模型的数据是未知的, $k=1,2,\dots,K$k=1,2,…,K 这里用**隐变量$\gamma_{jk}$γjk​**表示：
$\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}$γjk​的维度为$(j\times k)$(j×k)，每一行都是一种one-hot的形式，表示对于观测数据$j$j，在$\gamma_{jk}$γjk​的第j行的第k列为1，该行其余的列都为0。

完全数据为$(y_j,\gamma_{j1},\gamma_{j2},\dots,\gamma_{jK},k=1,2,\dots,N)$(yj​,γj1​,γj2​,…,γjK​,k=1,2,…,N)，完全数据的似然函数为：
$\begin{aligned} P(y,\gamma|\theta)=&\prod_{j=1}^NP(y_j,\gamma_{j1},\gamma_{j2},\dots,\gamma_{jK}|\theta)\\ =&\prod_{k=1}^K\prod_{j=1}^N\left[\alpha_k\phi(y_j|\theta_k)\right]^{\gamma_{jk}}\\ =&\prod_{k=1}^K\alpha_k^{n_k}\prod_{j=1}^N\left[\phi(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^2}\right)\right]^{\gamma_{jk}} \qquad (1)\\ \end{aligned}$P(y,γ∣θ)====​j=1∏N​P(yj​,γj1​,γj2​,…,γjK​∣θ)k=1∏K​j=1∏N​[αk​ϕ(yj​∣θk​)]γjk​k=1∏K​αknk​​j=1∏N​[ϕ(yj​∣θk​)]γjk​k=1∏K​αknk​​j=1∏N​[2π​σk​1​exp(−2σ2(yj​−μk​)2​)]γjk​(1)​
其中$n_k=\sum_{j=1}^N\gamma_{jk}, \sum_{k=1}^Kn_k=N$nk​=∑j=1N​γjk​,∑k=1K​nk​=N。$n_k$nk​是所有行第k列的和，可以理解为属于第k个高斯分布的样本的个数。

下面求完全数据对数似然函数：

我们假设上面似然函数的式子(1)中，$\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$∏j=1N​[2π​σk​1​exp(−2σk2​(yj​−μk​)2​)]γjk​=ωk​

则:
$\begin{aligned}\log P(y,\gamma|\theta) &= log\prod_{k=1}^K\alpha_k^{n_k}\omega_k\\&= \sum_{k=1}^K(n_klog\alpha_k + log\omega_k) \\\end{aligned}$logP(y,γ∣θ)​=logk=1∏K​αknk​​ωk​=k=1∑K​(nk​logαk​+logωk​)​
其中$log\omega_k$logωk​展开为：
$\begin{aligned}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{aligned}$logωk​​=logj=1∏N​[2π​σk​1​exp(−2σk2​(yj​−μk​)2​)]γjk​=j=1∑N​γjk​(log2π​1​−logσk​−2σk2​(yj​−μk​)2​)​
所以完全数据的对数似然函数为：
$\begin{aligned}\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{aligned}$logP(y,γ∣θ)=k=1∑K​{nk​logαk​+j=1∑N​γjk​[log(2π​1​)−logσk​−2σ21​(yj​−μk​)2]}​

2、EM算法的E步：确定Q函数

下面的Q 函数是完全数据的对数似然函数$\log P(y,\gamma|\theta)$logP(y,γ∣θ)关于给定观测数据$y$y的当前参数$\theta^{(i)}$θ(i)下，对未观测数据$\gamma$γ的条件概率分布$P(\gamma|y,\theta^{(i)})$P(γ∣y,θ(i))的期望表达式：
$\begin{aligned}Q(\theta,\theta^{(i)})=&E[\log P(y,\gamma|\theta)|y,\theta^{(i)}]\\=&\color{green}E\color{black}\left\{\sum_{k=1}^K\left\{\color{red}n_k\color{black}\log \alpha_k+\color{blue}\sum_{j=1}^N\gamma _{jk}\color{black}\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\} \qquad (1) \\=&\sum_{k=1}^K\left\{\color{green}E\color{black}\left\{\color{red}n_k\color{black}\log \alpha_k+\color{blue}\sum_{j=1}^N\gamma _{jk}\color{black}\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\} \qquad (2) \\=&\sum_{k=1}^K\left\{\color{green}E\color{black}\left\{\color{red}\sum_{j=1}^N\gamma_{jk}\color{black}\log \alpha_k+\color{blue}\sum_{j=1}^N\gamma _{jk}\color{black}\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\} \qquad (3) \\=&\sum_{k=1}^K\left\{\color{red}\sum_{j=1}^{N}(\color{green}E\color{red}\gamma_{jk})\color{black}\log \alpha_k+\color{blue}\sum_{j=1}^N(\color{green}E\color{blue}\gamma _{jk})\color{black}\left[\log \left(\frac{1}{\sqrt{2\pi}}\right)-\log \sigma _k-\frac{1}{2\sigma^2(y_j-\mu_k)^2}\right]\right\}\qquad (4) \\\end{aligned}$Q(θ,θ(i))=====​E[logP(y,γ∣θ)∣y,θ(i)]E{k=1∑K​{nk​logαk​+j=1∑N​γjk​[log(2π​1​)−logσk​−2σ2(yj​−μk​)21​]}}(1)k=1∑K​{E{nk​logαk​+j=1∑N​γjk​[log(2π​1​)−logσk​−2σ2(yj​−μk​)21​]}}(2)k=1∑K​{E{j=1∑N​γjk​logαk​+j=1∑N​γjk​[log(2π​1​)−logσk​−2σ2(yj​−μk​)21​]}}(3)k=1∑K​{j=1∑N​(Eγjk​)logαk​+j=1∑N​(Eγjk​)[log(2π​1​)−logσk​−2σ2(yj​−μk​)21​]}(4)​
说明：

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

(2)到(3)是因为$n_k=\sum_{j=1}^N\gamma_{jk}$nk​=∑j=1N​γjk​。

(3)到(4)是因为Q函数是完全数据的对数似然函数对隐变量的条件概率分布。

上面的式子需要对隐变量$\gamma_{j k}$γjk​求期望（记为$\hat\gamma_{jk}$γ^​jk​），与隐变量无关的项可以视为常数：
$\begin{aligned}\hat \gamma _{jk}&= \color{purple}E(\gamma_{jk}|y_j,\theta)\\ &=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}^KP(\gamma_{jk}=1,y_j|\theta)} \\ &=\frac{P(y_j|\color{red}\gamma_{jk}=1,\theta\color{black})\color{green}P(\gamma_{jk}=1|\theta)}{\sum_{k=1}^KP(y_j|\gamma_{jk}=1,\theta)P(\gamma_{jk}=1|\theta)}\\ &=\frac{\color{green}\alpha_k\color{black}\phi(y_j|\color{red}\theta_k)}{\sum_{k=1}^K\alpha_k\phi(y_j|\theta_k)} \end{aligned}$γ^​jk​​=E(γjk​∣yj​,θ)=P(γjk​=1∣yj​,θ)=P(yj​,θ)P(γjk​=1,yj​,θ)​=P(yj​∣θ)P(θ)P(γjk​=1,yj​∣θ)P(θ)​=P(yj​∣θ)P(γjk​=1,yj​∣θ)​=∑k=1K​P(γjk​=1,yj​∣θ)P(γjk​=1,yj​∣θ)​=∑k=1K​P(yj​∣γjk​=1,θ)P(γjk​=1∣θ)P(yj​∣γjk​=1,θ)P(γjk​=1∣θ)​=∑k=1K​αk​ϕ(yj​∣θk​)αk​ϕ(yj​∣θk​)​​

3. EM算法的E步

迭代的M步是求$Q(\theta,\theta^{(i)})$Q(θ,θ(i))对$\theta$θ的极大值，得到新一轮迭代的模型参数：
$\theta^{(i+1)}=\arg\max_\theta Q(\theta,\theta^{(i)})$θ(i+1)=argθmax​Q(θ,θ(i))
对$Q(\theta,\theta^i)$Q(θ,θi)求导，可以得到未知量$\mu_k, \hat{\sigma_k}^2$μk​,σk​^​2的偏导，使其偏导等于0，求解得到：
$\hat{\mu_k}=\frac{\sum_{j=1}^{N}\hat{\gamma_{jk}}y_i}{\sum_{j=1}^{N}\hat{\gamma_{jk}}} \\ \hat{\sigma_k}^2=\frac{\sum_{j=1}^{N}\hat{\gamma_{jk}}(y_i-\mu_k)^2}{\sum_{j=1}^{N}\hat{\gamma_{jk}}}$μk​^​=∑j=1N​γjk​^​∑j=1N​γjk​^​yi​​σk​^​2=∑j=1N​γjk​^​∑j=1N​γjk​^​(yi​−μk​)2​
对于$\hat{a_k}$ak​^​而言，它是属于各个分布的样本个数的占比：
$\hat{a_k}= \frac{n_k}{N} = \frac{\sum_{j=1}^{N}\hat{\gamma_{jk}} }{N}$ak​^​=Nnk​​=N∑j=1N​γjk​^​​

4. 停止条件

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