EM算法解析

本文介绍EM算法的相关推导过程，同时会介绍混合高斯模型。

1.EM算法

1.1 EM算法的推导

观测数据（不完全数据）X关于参数$\theta$的对数似然函数为：

$$\begin{align*} L(\theta) &=log P(x|\theta)\\ &=log\sum_{z}P(x,z|\theta)\\ &=log\sum_{z}[q(z)\frac{P(x,z|\theta)}{q(z)}]\\ &\geq\sum_{z}q(z)log\frac{P(x,z|\theta)}{q(z)}\\ &=\sum_{z}q(z)logP(x,z|\theta)+\sum_{z}q(z)logq(z) \end{align*}$$

另一方面：

$$\begin{align*} L(\theta) &= logP(x|\theta) \\ &=\sum_{z}q(z)logP(x|\theta) \\ &=\sum_{z}q(z)log\frac{P(x,z|\theta)}{P(z|x,\theta)}\\ &=\sum_{z}q(z)log\frac{P(x,z|\theta)q(z)}{p(z|x,\theta)q(z)} \\ &=\sum_{z}q(z)logP(x,z|\theta) + \sum_{z}log\frac{q(z)}{p(z|x,\theta)} + \sum_{z}q(z)logq(z)\\ &=\sum_{z}q(z)logP(x,z|\theta) + KL(q(z)||p(z|x,\theta)) + \sum_{z}q(z)logq(z) \end{align*}$$
令:
$$B(\theta,q(z))=\sum_{z}q(z)logP(x,z|\theta)+\sum_{z}q(z)logq(z)$$
根据（公式1）:
$$L(\theta)-B(\theta,q(z)) = KL(q(z)||p(z|x,\theta))$$
由上式可知， $B(\theta,q(z))$是 $L(\theta)$的一个下界。二者的差值即为概率分布 $q(z)$与 $p(z|x,\theta)$的KL距离。

当$\theta$固定时，$L(\theta)$的值确定。但是随着$q(z)$的变化，$B(\theta,q(z))$的值会不断变化。接下来要做的是调整$q(z)$的分布，使$B(\theta,q(z))$的值尽可能的靠近$L(\theta)$的值，即尽可能的提高$L(\theta)$的下界。

根据KL距离的定义，当$q(z)=p(z|x,\theta)$时，$KL(q(z)||p(z|x,\theta)) = 0$，有下式成立：

$$L(\theta)=B(\theta)$$
EM算法：
E步：
在时刻t的参数 $\theta^t$固定的情况下，求能使 $B(\theta^t,q(z);x)$达到最大的下界的期望值：
$$q^t(z) = argmax_{q(z)}B(\theta^t,q(z);x)=p(z|x,\theta^t)$$
M步：
在 $q^t(z)$固定的情况下，求使 $B(\theta^t,q(z);x)$取得极大值的新的 $\theta^{t+1}$，由于 $B(\theta^t,q(z);x)$的第二项与 $\theta$无关，因此只需要求使第一项达到最大值的 $\theta^{t+1}$即可。
$$\theta^{t+1}=argmax_{\theta}E_{q^t(z)}P(x,z|\theta)=argmax_{\theta}\sum_{z}q^t(z)logP(x,z|\theta)$$

1.2 EM算法的收敛性

假定$\theta^t$和$\theta^{t+1}$是第t次和t+1次迭代后的结果，求证：

$$L(\theta^{t+1}) \geq L(\theta^t)$$
证明：
（1）在选定 $\theta^t$后，由E步 $q^t(z) = p(z,x|\theta^t)$，则：
$$L(\theta^t) = E_{q^t(z)}P(x,z|\theta^t) + H(q^t(z))= \sum_{z}q^t(z)logP(x,z|\theta^t) + \sum_{z}q^t(z)logq^t(z)$$
（2）在经过M步的推导得到 $\theta^{t+1}$后：
$$L(\theta^{t+1}) \geq \sum_{z}q^t(z)logP(x,z|\theta^{t+1}) + \sum_{z}q^t(z)logq^t(z)$$
（3）根据M步的定义，将 $\theta^t$调整到 $\theta^{t+1}$后，使下界最大化：
$$E_{q^t(z)}P(x,z|\theta^{t+1}) \geq E_{q^t(z)}P(x,z|\theta^t)$$
（4）综合上述三个公式可得：
$$L(\theta^{t+1}) \geq L(\theta^t)$$
也就是说，经过一次次的迭代，EM算法总会收敛到一个局部最大值（但不一定能收敛到全局最大值）。

2.混合高斯模型（GMM）

假设数据{$x_1,x_2,……,x_n$}由混合高斯模型生成：

$$P(x,\theta)=\sum_{k=1}^K\theta_z\phi(x|\theta_x)=\sum_{k=1}^K\pi_kN(\mu_k,\sigma_k)$$
其中\theta_z由两部分组成：
$$\theta_z = \{\pi_1,\pi_2,……,\pi_k\}$$
$$\theta_x = \{\mu_1,\mu_2,……,\mu_k;\sigma_1,\sigma_2,……,\sigma_k\}$$

定义：

$$\begin{equation} Z = \left[ \begin{array}{c} z_{1}\\ z_{2}\\ …\\ z_{n} \end{array} \right] = %右括号 \left[ \begin{array}{c} 1 & 0 & 0 &…& 0\\ 0 & 1 & 0 &…& 0\\ …&… &…&…&…\\ 0 & 0 & 0 &…& 1 \end{array} \right] %右括号 \end{equation}$$
其中：
$$\begin{equation} Z_{nk}= \begin{cases} 1 & \text{第n个观察值来自第k个模型}\\ 0 & \text{否则} \end{cases} \end{equation}$$

有了观测数据$x_n$及未观测数据$z_{nk}$，那么完全数据是：

$$(x_n,z_{n1},z_{n2,…,z_{nk}}),n=1,2,3,…N$$
所以，完全数据的似然函数为：
$$\begin{align*} p(x,z|\theta) & = P(x_1,x_2,…,x_n,z_1,z_2,…,z_n|\theta)\\ & = \prod_{n=1}^NP(x_n,z_{n1},z_{n2,…,z_{nk}}|\theta)\\ & = \prod_{n=1}^NP(z_n|\theta_z)\prod_{n=1}^NP(x_n|z_n,\theta_x)\\ & = \prod_{n=1}^N\prod_{k=1}^KN(\mu_k,\sigma_k)^{Z_{nk}}\prod_{n=1}^N\prod_{k=1}^K\pi_k^{Z_{nk}} \end{align*}$$

对数似然函数为：

$$L(\theta) = logP(x,z|\theta) = \sum_{n=1}^N\sum_{k=1}^KZ_{nk}[logN(\mu_k,\sigma_k) + log\pi_k]$$