EM算法及高斯混合模型算法推导

期望最大化(EM)算法

1.前言

​ 概率模型有时候既含有观测变量，又含有隐变量。只含有观测变量的情况下，直接对观测值进行极大似然估计便能够求出参数；比如抛一枚不均匀硬币n次，极大似然估计能够求解出正反面分别出现的概率。在含有隐变量的情况下，无法通过极大似然估计求得；比如手中有三枚不均匀硬币，先从中选取一枚硬币，然后再抛，得到的正反面为观测值；如果直接用极大似然估计，无法体现选择硬币的过程，错误地将三枚硬币视为同一枚硬币。本文要介绍的期望最大化算法便是用来求解含有隐变量的概率模型。

2.算法导出

​ EM算法实质是极大似然估计的变型，是由极大化似然函数导出的。回到三硬币的例子，将选择硬币的事件记为隐变量Z，所求参数记为$\theta$，观测值硬币正反面记为$Y$，于是似然函数为
$$L(\theta )=\log P(Y|\theta )$$
​ 上式中不包含隐变量Z，于是引入隐变量Z可得
$$L(\theta )=\log P(Y|\theta )=\log \sum _{Z}P(Y,Z|\theta )$$
​ 上式极大化无法直接通过偏导求出，可采用迭代法求解，假设当前参数值为${\theta }_i$，迭代法希望每一步$\theta$都能够使得$L(\theta )-L({\theta }_i)$最大化，即
$$\begin{array}{rl}L(\theta )-L({\theta }_i)& =\log \sum _{Z}P(Y,Z|\theta )-\log P(Y|{\theta }_i)\\ & =\log \sum _Z\frac{P(Y,Z|\theta )}{P(Y|{\theta }_i)}\\ & =\log \sum _{Z}P(Z|Y,{\theta }_i)\frac{P(Y,Z|\theta )}{P(Y|{\theta }_i)P(Z|Y,{\theta }_i)}\end{array}$$
​ 又因为log函数是凹函数并且${\sum }_{Z}P(Z|Y,{\theta }_i)=1$，由琴生不等式可得
$$\begin{array}{rl}L(\theta )-L({\theta }_i)& =\log \sum _{Z}P(Z|Y,{\theta }_i)\frac{P(Y,Z|\theta )}{P(Y|{\theta }_i)P(Z|Y,{\theta }_i)}\\ & \ge \sum _{Z}P(Z|Y,{\theta }_i)\log \frac{P(Y,Z|\theta )}{P(Y|{\theta }_i)P(Z|Y,{\theta }_i)}\\ & =\sum _{Z}P(Z|Y,{\theta }_i)\log \frac{P(Y,Z|\theta )}{P(Y,Z|{\theta }_i)}\\ & =\sum _{Z}P(Z|Y,{\theta }_i)\log P(Y,Z|\theta )-\sum _{Z}P(Z|Y,{\theta }_i)\log P(Y,Z|{\theta }_i)\\ & \stackrel{△}{=}B(\theta ，{\theta }_i)\end{array}$$
​ 所以，$B(\theta ，{\theta }_i)$是$L(\theta )-L({\theta }_i)$的下界，即最大化下界等价于最大化原函数
$$\begin{array}{rl}\theta & =\arg \underset{\theta }{\max }L(\theta )-L({\theta }_i)\\ & =\arg \underset{\theta }{\max }B(\theta ,{\theta }_i)\\ & =\arg \underset{\theta }{\max }\sum _{Z}P(Z|Y,{\theta }_i)\log P(Y,Z|\theta )-\sum _{Z}P(Z|Y,{\theta }_i)\log P(Y,Z|{\theta }_i)\end{array}$$
​ 又因为上式中，后项与因变量$\theta$无关，可舍去，即
$$\begin{array}{rl}\theta & =\arg \underset{\theta }{\max }\sum _{Z}P(Z|Y,{\theta }_i)\log P(Y,Z|\theta )\\ & =\arg \underset{\theta }{\max }{E}_Z[\log P(Y,Z|\theta )|Y,{\theta }_i]\\ & \stackrel{△}{=}\arg \underset{\theta }{\max }Q(\theta ,{\theta }_i)\end{array}$$
​ 于是极大似然估计等价于最大化$Q(\theta ,{\theta }_i)$函数。

​ 综上所述，导出期望最大化算法:

• 第一步(E步)：求解期望Q函数
  $$Q(\theta ,{\theta }_i)=E_Z[\log P(Y,Z|\theta )|Y,{\theta }_i]$$

• 第二步(M步)：最大化Q函数

3.高斯混合模型

​ 接下来举例介绍下EM算法应该如何使用，高斯分布是常用的概率分布，但是在实际问题中有些概率模型并不是单一高斯分布能够表示。比如说身高问题，我们可以假设身高满足单一高斯分布状态，两头少中间多，但身高跟性别有一定关系，更好地假设是，男生女生身高各满足一个不同的高斯分布。假设男女生出现的概率是${\alpha }_i,i=0,1$,并且身高分别满足高斯分布$\phi ({\mu }_i,{\sigma }_i),i=0,1$，那么一个人身高为x的概率是$P(x)={\sum }_{i=0}^1{\alpha }_i\phi ({\mu }_i,{\sigma }_i)$。这样的概率分布就是混合高斯分布。

​ 混合高斯分布的一般形式是
$$P(Y|\alpha ,\mu ,\sigma )=\sum _{k=1}^K{\alpha }_k\phi ({\mu }_k,{\sigma }_k)$$
​ 假设观测值集合为Y，使用EM算法求解高斯混合模型:

​ 1.求解Q函数(E步)

​ 混合高斯分布可以视为两步，第一步以一定概率选择某一高斯分布，第二步根据该高斯分布确定概率。观测值记为Y，隐变量为选择哪一种高斯分布，记为$\gamma$，且定义
$${\gamma }_{jk}=\{\begin{array}{ll}1,& 第j个样本选择第k个高斯分布\\ 0,& 否则\end{array}$$
​ 完整似然函数为
$$\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[{\alpha }_k\phi ({\mu }_k,{\sigma }_k){]}^{{\gamma }_{jk}}\\ & =\prod _{j=1}^N\prod _{k=1}^K[{\alpha }_k\frac{1}{\sqrt{2\pi }{\delta }_k}{e}^{-\frac{(y_j-{\mu }_k{)}^2}{2{\delta }_k^2}}{]}^{{\gamma }_{jk}}\end{array}$$
​

​ 取对数得
$$\begin{array}{rl}\log P(Y,\Gamma |\theta )& =\sum _{j=1}^N\sum _{k=1}^K{\gamma }_{jk}\log [{\alpha }_k\frac{1}{\sqrt{2\pi }{\delta }_k}{e}^{-\frac{(y_j-{\mu }_k{)}^2}{2{\delta }_k^2}}]\\ & =\sum _{j=1}^N\sum _{k=1}^K{\gamma }_{jk}[\log {\alpha }_k-\log (\sqrt{2\pi })-\log {\delta }_k-\frac{(y_j-{\mu }_k{)}^2}{2{\delta }_k^2}]\end{array}$$
​ 于是，
$$\begin{array}{rl}Q(\theta ,{\theta }_i)& =E_{\Gamma }[\log P(Y,\Gamma |\theta )|Y,{\theta }_i]\\ & =\sum _{j=1}^N\sum _{k=1}^{K}E[{\gamma }_{jk}|Y,{\theta }_i][\log {\alpha }_k-\log (\sqrt{2\pi })-\log {\delta }_k-\frac{(y_i-{\mu }_k{)}^2}{2{\delta }_k^2}]\end{array}$$
​ 上式中需要计算$E[{\gamma }_{jk}|Y,{\theta }_i]$，记为${\hat{\gamma }}_{jk}$，即有
$$\begin{array}{rl}{\hat{\gamma }}_{jk}& =E[{\gamma }_{jk}|Y,{\theta }_i]\\ & =P({\gamma }_{jk}=1|Y,{\theta }_i)\end{array}$$
​ 由贝叶斯公式得
$$\begin{array}{rl}P({\gamma }_{jk}=1|Y,{\theta }_i)& =\frac{P(y_j|{\gamma }_{jk}=1,{\theta }_i)P({\gamma }_{jk}=1|{\theta }_i)}{{\sum }_{k=1}^{K}P(y_j|{\gamma }_{jk}=1,{\theta }_i)P({\gamma }_{jk}=1|{\theta }_i)}\\ & =\frac{\phi (y_j|{\theta }_k){\alpha }_k}{{\sum }_{k=1}^K\phi (y_j|{\theta }_k){\alpha }_k}\end{array}$$
​ 于是$\widehat{{\gamma }_{jk}}$表示响应度，即对于第j个样本数据，属于第k个高斯分布的概率，可以由上一步参数来进行计算。

​ 于是有，
$$\begin{array}{rl}Q(\theta ,{\theta }_i)& =\sum _{j=1}^N\sum _{k=1}^K\widehat{{\gamma }_{jk}}[\log {\alpha }_k-\log (\sqrt{2\pi })-\log {\delta }_k-\frac{(y_j-{\mu }_k{)}^2}{2{\delta }_k^2}]\end{array}$$
​ 2.求最大化操作(M步)

​ Q函数对$\mu ,\delta$分别求偏导可得
$$\frac{\partial Q}{\partial {\mu }_k}=\sum _{j=1}^N{\hat{\gamma }}_{jk}\frac{y_j-{\mu }_k}{{\delta }_k^2}$$

$$\frac{\partial Q}{\partial {\delta }_k}=\sum _{j=1}^N{\hat{\gamma }}_{jk}[-\frac{1}{{\delta }_k}+\frac{(y_j-{\mu }_k{)}^2}{{\delta }_k^3}]$$

​ 将以上两式得零可以推导出
$${\hat{\mu }}_k=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}\cdot y_j}{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}}，k=1,2,\dots ,K$$

$${\hat{\delta }}_k^2=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}\cdot (y_j-{\mu }_k{)}^2}{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}}，k=1,2,\dots ,K$$

​ 结合约束条件${\sum }_{k=1}^K{\alpha }_k=1$，求Q函数拉格朗日函数有
$$L(\lambda ,\alpha )=Q(\alpha )+\lambda \cdot (\sum _{k=1}^K{\alpha }_k-1)$$
​ 分别对$\alpha ，\lambda$求偏导得
$$\frac{\partial L}{\partial {\alpha }_k}=\sum _{j=1}^N{\hat{\gamma }}_{jk}\frac{1}{{\alpha }_k}+\lambda ,k=1,2,\dots ,K$$

$$\frac{\partial L}{\partial \lambda }=\sum _{k=1}^K{\alpha }_k-1$$
​ 令上两式为零，可得
$${\hat{\alpha }}_k=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}}{{\sum }_{j=1}^N{\sum }_{k=1}^K{\hat{\gamma }}_{jk}}=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}}{N}，k=1,2,\dots ,K$$
​ 综上所述，高斯混合模型(GMM)模型算法流程如下：

• 1.初始化各参数

• 2.按照上一次参数更新响应度
  $${\hat{\gamma }}_{jk}=\frac{\phi (y_j|{\theta }_k){\alpha }_k}{{\sum }_{k=1}^K\phi (y_j|{\theta }_k){\alpha }_k}，k=1,2,\dots ,K;j=1,2,\dots ,N$$

• 3.更新模型参数

  $${\hat{\mu }}_k=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}\cdot y_j}{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}}，k=1,2,\dots ,K$$

  $${\hat{\delta }}_k^2=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}\cdot (y_j-{\mu }_k{)}^2}{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}}，k=1,2,\dots ,K$$

  $${\hat{\alpha }}_k=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}}{{\sum }_{j=1}^N{\sum }_{k=1}^K{\hat{\gamma }}_{jk}}=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}}{N}，k=1,2,\dots ,K$$

  ​
  ​

• 4.重复上述过程直至收敛

4.参考资料

• 统计学习方法 - 李航