高斯混合模型（Guassian Mixed Model, 简称GMM）是一种常见的聚类方法，与K-means类似，同样使用了EM算法进行迭代计算。
高斯混合模型
高斯混合模型假设每个簇的数据都是符合高斯分布的，当前数据呈现的分布就是各个簇的高斯分布叠加在一起的结果。于是高斯混合模型采用线性组合来对数据分布进行拟合，理论上高斯混合模型可以拟合出任意类型的分布。
具体的，我们想要将数据分为K个簇（cluster），可以假设不同簇中的样本各自服从不同的高斯分布，K个簇就有K个高斯分布，将K个高斯分布加起来就得到了GMM模型。对于每一个高斯模型，其均值\(\mu_k\)和方差\(\sigma_k\)都是待估计的参数，此外，每个模型还有一个参数\(\alpha_k\)，表示该模型的权重，该参数满足\(\sum_{k=1}^K \alpha_k =1\)。 \[ \begin{align*} p(x) &=\sum_{k=1}^K P(k)P(x\mid k)\\ &=\sum_{k=1}^K\alpha_k \mathcal{N}({x}\mid{\mu}_k,{\sigma}_k)\tag{2-1}\\ \end{align*} \] 其中， \[ \begin{align*} \mathcal{N}({\bf x}\mid{\mu}_k,{\sigma}_k) = \frac{1}{\sigma_k\sqrt{2\pi}} \exp\left(-\frac{(x-\mu_k)^2}{2\sigma_k^2}\right) \\ \sum_{k=1}^K \alpha_k=1 \end{align*} \] 若\(\bf x\)为向量，则可以写成： \[ \begin{align*} p({\bf x}) &= \sum_{k=1}^K\alpha_k \mathcal{N}({\bf x}\mid{\bf \mu}_k,{\bf \Sigma}_k)\\ \mathcal{N}({\bf x}\mid{\bf \mu}_k,{\bf \Sigma}_k) &= \frac{1}{\sqrt{(2\pi)^K|\Sigma_k|}} \exp\left(-\frac{1}{2}({\bf x} - {\bf \mu}_k)^T\Sigma_k^{-1}({\bf x} - {\bf \mu_k})\right) \\ \sum_{k=1}^K \alpha_k&=1 \end{align*} \] 为了方便讲解，下面我们采用非向量的方式进行讲解。
高斯混合模型是一个生成式的模型，数据可以看成是从多个高斯分布生成出来的。具体的说，对于一个样本，按照概率\(\alpha_k\)选择第k个分布，然后从第k分布\(\mathcal{N}({x}\mid{\mu}_k,{\sigma}_k)\)生成数据x。
在实际中，我们已知的是样本点，要来估计上面提到的参数\(\alpha_k, \ \mu_k, \ \bf \sigma_k\), 假设我们采用极大似然法，可以写出对数似然函数如下： \[ l = \sum_{i=1}^N \log P(x_i) = \sum_{i=1}^N \log \left(\sum_{k=1}^K \alpha_k \mathcal{N}({x}_i\mid{\mu}_k,{\sigma}_k) \right) \tag{2-2} \] 这是一个复杂的非凸函数。对于每个观测数据点来说，事先并不知道它是属于哪个子分布的（hidden variable）。此外，\(\log P\)中的\(P\)里面还有求和项， K 个高斯模型的和不是一个高斯模型，对于每个子模型都有未知的, \(\alpha_k, \ \mu_k, \ \bf \sigma_k\)直接求导无法计算，需要通过迭代的方法求解。
因此，求解高斯混合模型采用的是EM算法：
- E-step： 根据当前的参数，计算每个点由某个分模型生成的概率
- M-step：使用E步估计出的概率，来改进每个分模型的\(\alpha_k, \ \mu_k, \ \bf \sigma_k\)
高斯分布和K-means都需要指定K值，都是用EM算法来求解，同时往往也只能收敛到局部最优。它比K-means更优的地方在于：
- 可以给出一个样本属于某个类的概率是多少
- 除了聚类外，还可以用于概率密度的估计
- 可以用于生成新的样本点
求解推导
本小节主要是推导EM算法中的E-step和M-step
引入隐变量
首先对式2-1做一下变换，引入一个服从多项分布的随机变量\(\bf z\)，简单的说\(\bf z\)就是一个One-hot的向量: \[ z_{k} = \begin{cases} 1，&\text{样本来自第k个分模型} \\ 0, &\text{否则} \end{cases} \]
假设\(z_k\)之间是独立同分布的，则容易写出联合概率： \[ p({\bf z}) = p(z_1) p(z_2)...p(z_K) = \prod_{k=1}^K \alpha_k^{z_k} \tag{2-2} \] 则给定\(\bf z\)之后，可以求出\(P(x \mid z)\)(相当于选择第k个分布)： \[ p(x\mid {\bf z}) = \prod_{k=1}^K \mathcal{N}({\bf x}\mid{\mu}_k,{\bf \sigma}_k)^{z_k} \tag{2-3} \] 于是可以根据联合概率求出： \[ \begin{align*} p({\bf x}) &= \sum_z p(z)p(x\mid z)\\ &= \sum_{k=1}^K\alpha_k \mathcal{N}({\bf x}\mid{\mu}_k,{\bf \Sigma}_k) \tag{2-4} \end{align*} \] 2-4是忽略了\(z_k=0\)的项，因为值为1。可以发现和2-1等价。
这样，我们引入了隐变量\(\bf z\)，用来描述样本是由第j个高斯分布产生的。
E-step
在EM算法中，我们提到E-step求的Q函数形式为\(P(z \mid x)\)，对于GMM模型来说，我们可以用贝叶斯公式求出\(p(z_k = 1 \mid x)\)： \[ \begin{align*} \gamma(z_k) = P(z_k=1\mid x) &= \frac{P(x, z_k=1)}{P(x)} \\ &= \frac{P(x\mid z_k =1)P(z_k =1)}{P(x)}\\ &= \frac{\alpha_k \mathcal{N}({\bf x}\mid{\mu}_k,{\bf \Sigma}_k)}{\sum_{j=1}^K\alpha_j \mathcal{N}({\bf x}\mid{\mu}_j,{\bf \Sigma}_j)} \tag{2-5} \end{align*} \]
M-Step
现在我们知道了2-5的形式，就可以求导来优化\(\mu\)等参数了。回顾一下对数似然函数2-2： \[ l = \sum_{i=1}^N \log P(x_i) = \sum_{i=1}^N \log \left(\sum_{k=1}^K \alpha_k \mathcal{N}({x}_i\mid{\mu}_k,{\sigma}_k) \right) \tag{2-2} \] 2-2对\(\mu_k\)求导，并令导数为0，有 \[ \begin{align*} \frac{\partial l}{\partial \mu_k} &= \sum_{i=1}^N \frac{\alpha_k \frac{\partial \mathcal{N}({x}_k\mid{\mu}_k,{\sigma}_k)}{\partial \mu_k}}{\sum_{j=1}^K \alpha_j \mathcal{N}({x}_i\mid{\mu}_j,{\sigma}_j)}\\ &=\sum_{i=1}^N \frac{\alpha_k \frac{1}{\sigma_k\sqrt{2\pi}} \exp\left(-\frac{(x_i-\mu_k)^2}{2\sigma_k^2}\right) (\frac{1}{\sigma_k^2}(x_i-\mu_k)(-1)) }{\sum_{j=1}^K \alpha_j \mathcal{N}({x}_i\mid{\mu}_j,{\sigma}_j)}\\ &= -\sum_{i=1}^N \frac{\alpha_k \mathcal{N}({x}_i\mid{\mu}_k,{\sigma}_k) }{\sum_{j=1}^K \alpha_j \mathcal{N}({x}_i\mid{\mu}_j,{\sigma}_j)} \frac{x_i-\mu_k}{\sigma_k^2}\\ &= -\sum_{i=1}^N \gamma(z_{ik})\frac{x_i-\mu_k}{\sigma_k^2} = 0 \end{align*} \] 得： \[ \mu_k = \frac{\sum_{i=1}^N \gamma(z_{ik})\ x_i}{\sum_{i=1}^N \gamma(z_{ik})} \] 2-2对\(\sigma_k\)求导有： \[ \begin{align*} \frac{\partial l}{\partial \sigma_k} &= \sum_{i=1}^N \frac{\alpha_k \frac{\partial \mathcal{N}({x}_k\mid{\mu}_k,{\sigma}_k)}{\partial \sigma_k}}{\sum_{j=1}^K \alpha_j \mathcal{N}({x}_i\mid{\mu}_j,{\sigma}_j)}\\ &= \sum_{i=1}^N \frac{\alpha_k \left( -\frac{1}{\sqrt{2\pi}\sigma_k^2}\exp(-\frac{(x_i-\mu_k)^2}{2\sigma_k^2}) + \frac{1}{\sqrt{2\pi} \sigma_k} \exp(-\frac{(x_i-\mu_k)^2}{2\sigma_k^2})\frac{(x_i-\mu_k)^2}{2\sigma_k^4} 2\sigma_k\right)}{\sum_{j=1}^K \alpha_j \mathcal{N}({x}_i\mid{\mu}_j,{\sigma}_j)}\\ &= \sum_{i=1}^N \frac{\alpha_k \mathcal{N}({x}_i\mid{\mu}_k,{\sigma}_k)\left( -\frac{1}{\sigma_k}+ \frac{(x_i-\mu_k)^2}{\sigma_k^3} \right)}{\sum_{j=1}^K \alpha_j \mathcal{N}({x}_i\mid{\mu}_j,{\sigma}_j)}\\ &= \sum_{i=1}^N \gamma(z_{ik})\left( -\frac{1}{\sigma_k}+ \frac{(x_i-\mu_k)^2}{\sigma_k^3} \right) = 0 \end{align*} \] 令其导数为0，有 \[ \sigma_k^2 = \frac{\sum_{i=1}^N \gamma({z_{ik}})(x_i-\mu_k)^2}{\sum_{i=1}^N\gamma({z_{ik}})} \] 由于\(\sum_k \alpha_k = 1\)，这是有约束条件的，所以需要引入拉格朗日算子然后在对\(\alpha_k\)求偏导： \[ \begin{align*} l_2 &= \sum_{i=1}^N \log \left(\sum_{k=1}^K \alpha_k \mathcal{N}({x}_i\mid{\mu}_k,{\sigma}_k) \right) + \sum_{k=1}^K\lambda_k (\alpha_k - 1)\\ \frac{\partial l_2}{\partial \alpha_k} &= \sum_{i=1}^N \frac{ \mathcal{N}({x}_k\mid{\mu}_k,{\sigma}_k)}{\sum_{j=1}^K \alpha_j \mathcal{N}({x}_i\mid{\mu}_j,{\sigma}_j)} + \lambda= 0 \end{align*} \] 两边同乘以\(\alpha_k\)得到： \[ \sum_{i=1}^N \frac{ \alpha_k\mathcal{N}({x}_k\mid{\mu}_k,{\sigma}_k)}{\sum_{j=1}^K \alpha_j \mathcal{N}({x}_i\mid{\mu}_j,{\sigma}_j)}+\lambda \alpha_k = 0\\ \sum_{i=1}^N\gamma(z_{ik}) = -\lambda a_k \] 这里有个技巧就是对k求和，有 \[ \begin{align*}\sum_{k=1}^ K\sum_{i=1}^N \gamma(z_{ik}) &= \sum_{i=1}^N \frac{ \sum_{k=1}^{K} \alpha_k\mathcal{N}({x}_k\mid{\mu}_k,{\sigma}_k)}{\sum_{j=1}^K \alpha_j \mathcal{N}({x}_i\mid{\mu}_j,{\sigma}_j)} \\ &= N \\ &= -\lambda \sum_{k=1}^K a_k\\ &= -\lambda \end{align*} \] 因此，有 \[ \alpha_k = \frac{\sum_{i=1}^N \gamma(z_{ik})}{N} \]
小结
GMM模型采用EM算法进行求解，步骤如下：
- 设定K的值，随机初始化K个高斯分布的参数\(\alpha_k, \mu_k, \sigma_k\)
- E-step，根据当前的\(\alpha_k\)计算后验概率\(\hat\gamma(z_{ik})\)
\[ \begin{align*} \hat\gamma(z_{ik}) = P(z_k=1\mid x_i) &= \frac{\alpha_k \mathcal{N}({\bf x_i}\mid{\mu}_k,{\bf \Sigma}_k)}{\sum_{j=1}^K\alpha_j \mathcal{N}({\bf x_i}\mid{\mu}_j,{\bf \Sigma}_j)} \end{align*} \]
- M-step, 根据E-step中计算的\(\gamma(z_{ik})\)再计算新的\(\hat\alpha_k, \hat\mu_k, \hat\sigma_k\)
\[ \begin{align*} \hat \mu_k &= \frac{\sum_{i=1}^N \hat\gamma(z_{ik})\ x_i}{\sum_{i=1}^N \gamma(z_{ik})}\\ \hat \sigma_k^2 &= \frac{\sum_{i=1}^N \hat\gamma({z_{ik}})(x_i-\mu_k)^2}{\sum_{i=1}^N\hat\gamma({z_{ik}})}\\ \hat\alpha_k &= \frac{\sum_{i=1}^N \hat\gamma(z_{ik})}{N} \end{align*} \]
- 计算对数似然函数\(\sum_{i=1}^N \log \left(\sum_{k=1}^K \alpha_k \mathcal{N}({x}_i\mid{\mu}_k,{\sigma}_k) \right)\)，看是否收敛，若不收敛，返回2步
参考资料
[1] 《百面机器学习》 [2]高斯混合模型（GMM）及其EM算法的理解
https://brilliant.org/wiki/gaussian-mixture-model/
https://www.geeksforgeeks.org/gaussian-mixture-model/
https://blog.csdn.net/lpsl1882/article/details/53072907
