EM算法应用：高斯混合模型

---

高斯混合模型概要

假设我们有若干个高斯分布模型，每个模型会以一定的概率被选择，然后根据被选择的模型的概率分布产生一个结果。我们希望可以在只有输出序列的情况下，对模型的参数进行估计。

形式化

我们把观察序列记为$Y=y_1,y_2,…,y_j$ ，将模型总数记为 $K$ ，第 $k$ 个高斯分布模型 $\phi\left(y | \theta_{k}\right)$ 会以 $\alpha_k$ 的概率被选择，

$$P(y | \theta)=\sum_{k=1}^{K} \alpha_{k} \phi\left(y | \theta_{k}\right)$$

其中$\phi\left(y | \theta_{k}\right)=\frac{1}{\sqrt{2 \pi} \sigma_{k}} \exp \left(-\frac{\left(y-\mu_{k}\right)^{2}}{2 \sigma_{k}^{2}}\right)$,$\sigma_k$和$\mu_k$是第$k$个模型的参数。在这个问题中，每次生成观察结果$y_j$时使用了哪一个模型是隐变量，用$\gamma_{j k}$表示，定义如下:

$$\gamma_{j k}=\left\{\begin{array}{l}{1,y_j是由模型k生成的} \\ {0,否则}\end{array}\right.$$ $$j=1,2, \cdots, N ; \quad k=1,2, \cdots, K$$

构建模型

根据上面的定义，我们可以得到完全数据的似然函数（联合概率）:

$$P(y, \gamma | \theta)=\prod_{j=1}^{N} P\left(y_{j}, \gamma_{j 1}, \gamma_{j 2}, \cdots, \gamma_{j k} | \theta\right)\\ =\prod_{k=1}^{K} \prod_{j=1}^{N}\left[\alpha_{k} \phi\left(y_{j} | \theta_{k}\right)\right]^{\gamma_{j k}}\\ =\prod_{k=1}^{K} \alpha_{k}^{n_{k}} \prod_{j=1}^{N}\left[\phi\left(y_{j} | \theta_{k}\right)\right]^{\gamma_{k}}\\ =\prod_{k=1}^{K} \alpha_{k}^{n_{k}} \prod_{j=1}^{N}\left[\frac{1}{\sqrt{2 \pi} \sigma_{k}} \exp \left(-\frac{\left(y_{j}-\mu_{k}\right)^{2}}{2 \sigma_{k}^{2}}\right)\right]^{\gamma_{k}}$$

其中，$n_{k}=\sum_{j=1}^{N} \gamma_{j k}, \sum_{k=1}^{K} n_{k}=N$，$n_k$代表第$k$个模型被选择的次数

接着，我们可以写出对数似然函数：

$$\log P(y, \gamma | \theta)=\sum_{k=1}^{K} \left\{ n_{k} \log \alpha_{k}+\sum_{j=1}^{N} \gamma_{j k}\left[\log \left(\frac{1}{\sqrt{2 \pi}}\right)-\log \sigma_{k}-\frac{1}{2 \sigma_{k}^{2}}\left(y_{j}-\mu_{k}\right)^{2}\right]\right\}$$

根据EM算法的框架，我们先来求解Q函数

$$Q\left(\theta, \theta^{(i)}\right)=E_{Z}\left[\log P(Y, Z | \theta) | Y, \theta^{(i)}\right]\\ =\sum_{Z} P(Z | Y, \theta^{(j)}) \log P(Y, Z | \theta)$$

我们可以直接对对数似然函数求期望得到Q函数（第一行），也可以使用第二行的表达式，利用概率分布来直接求得Q函数。下面两种方法都会列出。

求解Q函数（E步）

通过期望求解

我们直接求解对数似然函数的期望，注意此时我们的随机变量是$\gamma$

$$Q\left(\theta, \theta^{(i)}\right)=E\left[\log P(y, \gamma | \theta) | y, \theta^{(i)}\right]\\ =E\left\{\sum_{k=1}^{K} \left\{ n_{k} \log \alpha_{k}+\sum_{j=1}^{N} \gamma_{j k}\left[\log \left(\frac{1}{\sqrt{2 \pi}}\right)-\log \sigma_{k}-\frac{1}{2 \sigma_{k}^{2}}\left(y_{j}-\mu_{k}\right)^{2}\right]\right\}\right\}\\ =E\left\{\sum_{k=1}^{K} n_{k} \log \alpha_{k} + \sum_{k=1}^{K} \sum_{j=1}^{N} \gamma_{j k}\left[\log \left(\frac{1}{\sqrt{2 \pi}}\right)-\log \sigma_{k}-\frac{1}{2 \sigma_{k}^{2}}\left(y_{j}-\mu_{k}\right)^{2}\right]\right\}\\ =E\left\{\sum_{k=1}^{K} \sum_{j=1}^{N} \gamma_{jk} \log \alpha_{k} + \sum_{k=1}^{K} \sum_{j=1}^{N} \gamma_{j k}\left[\log \left(\frac{1}{\sqrt{2 \pi}}\right)-\log \sigma_{k}-\frac{1}{2 \sigma_{k}^{2}}\left(y_{j}-\mu_{k}\right)^{2}\right]\right\}\\ =\sum_{k=1}^{K}\left\{\sum_{j=1}^{N}\left(E \gamma_{j k}\right) \log \alpha_{k}+\sum_{j=1}^{N}\left(E \gamma_{j k}\right)\left[\log \left(\frac{1}{\sqrt{2 \pi}}\right)-\log \sigma_{k}-\frac{1}{2 \sigma_{k}^{2}}\left(y_{j}-\mu_{k}\right)^{2}\right]\right\}$$

在这个式子中，$E\gamma_{jk}$是直接根据当前参数计算的，$\alpha_k,\sigma_k,\mu_k$是第$k$个模型待更新的参数。我们来计算$E\gamma_{jk}$,记为$\hat{\gamma}_{j k}$:

$$\hat{\gamma}_{j k}=E\left(\gamma_{j k} | y, \theta\right)=P\left(\gamma_{j k}=1 | y, \theta\right)\\ =\frac{P\left(\gamma_{j k}=1, y_{j} | \theta\right)}{ P\left( y_{j} | \theta\right)}\\ =\frac{P\left(\gamma_{j k}=1, y_{j} | \theta\right)}{\sum_{k=1}^{K} P\left(\gamma_{j k}=1, y_{j} | \theta\right)}\\ =\frac{ P\left(\gamma_{j k}=1 | \theta\right)P\left(y_{j} | \gamma_{j k}=1, \theta\right)}{\sum_{k=1}^{K} P\left(\gamma_{j k}=1 | \theta\right)P\left(y_{j} | \gamma_{j k}=1, \theta\right) }\\ =\frac{\alpha_{k} \phi\left(y_{j} | \theta_{k}\right)}{\sum_{k=1}^{K} \alpha_{k} \phi\left(y_{j} | \theta_{k}\right)}, \quad j=1,2, \cdots, N ; \quad k=1,2, \cdots, K$$

其中，第二行到第三行使用了全概率公式，第三行到第四行是条件概率公式。我们由此可以得到Q函数：

$$Q\left(\theta, \theta^{(i)}\right)=\sum_{k=1}^{K}\sum_{k=1}^{N} \left\{ \hat{\gamma}_{j k} \log \alpha_{k}+ \hat{\gamma}_{j k}\left[\log \left(\frac{1}{\sqrt{2 \pi}}\right)-\log \sigma_{k}-\frac{1}{2 \sigma_{k}^{2}}\left(y_{j}-\mu_{k}\right)^{2}\right]\right\}$$

通过概率分布求解

我们由

$$Q\left(\theta, \theta^{(i)}\right)=\sum_{Z} P(Z | Y, \theta^{(j)}) \log P(Y, Z | \theta)$$

得到我们需要求解的Q函数为

$$Q\left(\theta, \theta^{(i)}\right)=\sum_{j=1}^{N}\left[\sum_{k=1}^{K} P(Z=k | Y, \theta^{(j)}) \log P(Y, Z=k | \theta)\right]$$

其中$Z$代表选择的模型编号。
我们先来求解根据当前参数估算的条件概率部分

$$P(Z=k | Y, \theta^{(j)})=P(Z=k | Y=y_j, \theta^{(j)})\\ =\frac{P(Z=k,Y=y_j|\theta^{(i)})} {P(Y=y_j|\theta^{(i)})}$$

代入上面的概率分布定义得到

$$P(Z=k,Y=y_j|\theta^{(i)})= \alpha_{k} \phi\left(y | \theta_{k}\right)\\ P(Y=y_j|\theta^{(i)})=\sum_{k=1}^{K} \alpha_{k} \phi\left(y | \theta_{k}\right)\\ \frac{P(Z=k,Y=y_j|\theta^{(i)})} {P(Y=y_j|\theta^{(i)})}=\frac{\alpha_{k} \phi\left(y_{j} | \theta_{k}\right)}{\sum_{k=1}^{K} \alpha_{k} \phi\left(y_{j} | \theta_{k}\right)}$$

可见$P(Z=k | Y, \theta^{(j)})$即为$\gamma_{jk}$。
再次代入联合分布概率到$\log P(Y, Z=k | \theta)$中，我们可以得到和上面直接求期望一样的Q函数：

$$Q\left(\theta, \theta^{(i)}\right)=\sum_{k=1}^{K}\sum_{k=1}^{N} \left\{ \hat{\gamma}_{j k} \log \alpha_{k}+ \hat{\gamma}_{j k}\left[\log \left(\frac{1}{\sqrt{2 \pi}}\right)-\log \sigma_{k}-\frac{1}{2 \sigma_{k}^{2}}\left(y_{j}-\mu_{k}\right)^{2}\right]\right\}$$

极大化Q函数（M步）

我们分别求Q函数对参数的偏导并使偏导数为0可以得到新的参数的表达式。
对于$\mu,\sigma^2$,求导求极值过程无特别难点，结果如下：

$$\hat{\mu}_{k}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k} y_{j}}{\sum_{j=1}^{N} \hat{\gamma}_{j k}}, \quad k=1,2, \cdots, K\\ \hat{\sigma}_{k}^{2}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k}\left(y_{j}-\mu_{k}\right)^{2}}{\sum_{j=1}^{N} \hat{\gamma}_{j k}}, \quad k=1,2, \cdots, K$$

对于参数$\alpha_k$，我们需要利用额外的约束$\sum_{k=1}^{K}\alpha_k=1$来求解。
构造拉格朗日函数

$$L(\alpha,\lambda) = Q + \lambda(1-\sum_{k=1}^{K}\alpha_k)$$

分别对每一个$\alpha_k$和$\lambda$求偏导，得到

$$\frac{\partial L}{\partial \alpha_1} = \sum_{j=1}^{N}\frac{\gamma_{j1}}{\alpha_1}-\lambda\\ \cdots\\ \frac{\partial L}{\partial \alpha_K} = \sum_{j=1}^{N}\frac{\gamma_{jK}}{\alpha_K}-\lambda\\ \frac{\partial L}{\partial \lambda} = 1-\sum_{k=1}^{K}\alpha_k$$

令上述式子等于0，并求解方程组可以得到

$$\lambda = \sum_{k=1}^{K}\sum_{k=1}^{N}\gamma_{jk}\\ \hat\alpha_k = \frac{\sum_{j=1}^N \gamma_{jk}}{\sum_{k=1}^{K}\sum_{k=1}^{N}\gamma_{jk}}=\frac{n_k}{N}$$

至此，高斯混合模型推导完毕。

模型算法总结

输入：观测数据，高斯混合模型
输出：高斯混合模型参数
（1）初始化参数开始迭代
（2）E步：根据当前模型参数计算分模型$k$对观测变量$y_j$的响应度:

$$\hat{\gamma}_{j k}=\frac{\alpha_{k} \phi\left(y_{j} | \theta_{k}\right)}{\sum_{k=1}^{K} \alpha_{k} \phi\left(y_{j} | \theta_{k}\right)}, \quad j=1,2, \cdots, N ; \quad k=1,2, \cdots, K$$

(3) M步：更新参数

$$\hat{\mu}_{k}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k} y_{j}}{\sum_{j=1}^{N} \hat{\gamma}_{j k}}, \quad k=1,2, \cdots, K\\ \hat{\sigma}_{k}^{2}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k}\left(y_{j}-\mu_{k}\right)^{2}}{\sum_{j=1}^{N} \hat{\gamma}_{j k}}, \quad k=1,2, \cdots, K\\ \hat{\alpha}_{k}=\frac{\sum_{j=1}^{N} \hat{\gamma}_{j k}}{N}, \quad k=1,2, \cdots, K$$

(4) 重复(2)(3)直到收敛

参考文献

1. 《统计学习方法》,李航