EM算法的全称是Expectation Maximization algorithm。它主要是为了解决存在隐变量Z的情况下如何去估计$P(Y|\theta)$。
在实际中可能是对模型做极大思然估计，也可能是做极大后验概率的估计。资料来自李航的小蓝书的第九章，EM算法及其推广，因为在VAE模型中，如何去估计参数，也是参照了EM算法的思想。去优化对数似然函数的一个下界。

三硬币模型

首先先看Bernoulli Distribution的一些定义，比如一个硬币正面计作Y=1反面计作Y=0，抛一枚硬币，是正面的概率为$\Phi$。那么就有

$P(Y=y) = \Phi^y(1-\Phi)^{(1-y)}$

在三硬币模型中

$P(Y|\theta) = \sum_{Z} P(Y,Z|\theta)P(Z|\theta) = \prod_{j}^{n} [\pi p^y_j (1-p)^{1-y_j} + (1-\pi)q^{y_j}(1-q)^{1-y_j}]$

隐变量指的是针对分布$P(Z|\theta)$采样出来的变量，这里记作$\mu_i$，其中$i$表示第$i$的样本，表示这个是由硬币B还是硬币C得到的。

$\mu_j^{(i+1)} = \frac{\pi_{(i)} (p{(i)})^{y_j} (1-p^{(i)})^{(1-y_j)}}{\pi_{(i)} (p{(i)})^{y_j} (1-p^{(i)})^{(1-y_j)} + (1-\pi_{(i)}) (q{(i)})^{y_j} (1-q^{(i)})^{(1-y_j)}}$

这一步是对隐变量的期望的估计，这就是EM算法的E步。接下去是M步，这一步是对Q函数的最大化

$$\pi^{(i+1)} = \frac{1}{n} \sum_{j=1}^{n} \mu_j^{(i+1)} \\ p^{(i+1)} = \frac{\sum_{j=1}^n \mu_j^{(i+1)} y_j}{\sum_{j=1}^n \mu_j^{(i+1)}} \\ q^{(i+1)} = \frac{\sum_{j=1}^n (1-\mu_j^{(i+1)}) y_j}{\sum_{j=1}^n (1-\mu_j^{(i+1)})}$$

EM算法原理

这里设置我们的对数似然函数

$$L(\theta) = log P(Y|\theta) \\ = log \sum_Z P(Y,Z|\theta) \\ = log (\sum_Z P(Y|Z,\theta)P(Z|\theta))$$

我们选择使用迭代的方式来更新的话，就说明我们需要让$L(\theta) > L(\theta^{(i)})$

$$L(\theta) - L(\theta^{(i)}) = log (\sum_Z P(Y|Z,\theta)P(Z|\theta)) - log P(Y|\theta^{(i)}) \\ = log (\sum_Z P(Y|Z,\theta^{(i)}) \frac{\sum_Z P(Y|Z,\theta)P(Z|\theta)}{\sum_Z P(Y|Z,\theta^{(i)})}) - log P(Y|\theta^{(i)}) \\ \geq \sum_Z P(Y|Z,\theta^{(i)}) log\frac{\sum_Z P(Y|Z,\theta)P(Z|\theta)}{\sum_Z P(Y|Z,\theta^{(i)})} - logP(Y|\theta^{(i)}) \\ = \sum_Z P(Y|Z,\theta^{(i)}) log\frac{\sum_Z P(Y|Z,\theta)P(Z|\theta)}{\sum_Z P(Y|Z,\theta^{(i)})P(Y|\theta^{(i)})}$$

由此我们可以得到一个$L(\theta)$的一个下界

$B(\theta, \theta^{(i)}) = L(\theta^{(i)}) + \sum_Z P(Y|Z,\theta^{(i)}) log\frac{\sum_Z P(Y|Z,\theta)P(Z|\theta)}{\sum_Z P(Y|Z,\theta^{(i)})P(Y|\theta^{(i)})}$

同时

$B(\theta^{(i)}, \theta^{(i)}) = L(\theta^{(i)})$

所以在$i+1$步迭代的时候

$$\begin{align} \theta^{(i+1)} & = argmax_{\theta} B(\theta, \theta^{(i)}) \\ & = argmax_{\theta} (\sum_Z P(Z|Y, \theta^{(i)})log(P(Y|Z, \theta)P(Z|\theta))) \\ & = argmax_{\theta} (\sum_Z P(Z|Y, \theta^{(i)}) logP(Y,Z|\theta)) \\ & = argmax_{\theta} Q(\theta, \theta^{(i)}) \end{align}$$

因此Q函数实际上就是去掉了B函数在求极大值时不含参数的项。

# Gaussian Mixture Model
比如我们有三个一维的高斯模型

$$\theta_1 = (\alpha_1, \mu_1, \sigma_1^2) \\ \theta_2 = (\alpha_2, \mu_2, \sigma_2^2) \\ \theta_3 = (\alpha_3, \mu_3, \sigma_3^2)$$

对于一个观测序列$y_1, y_2,…, y_N$，从高斯混合中观测到$y_i$的概率是

$P(y_i|\theta_1, \theta_2, \theta_3) = \sum_{k} \alpha_k \Phi(y_i|\mu_k, \sigma_k)$

这个其实就是对应公式

$P(Y|\theta) = \sum_Z P(Z|\theta)*P(Y|Z, \theta)$

对应的Q函数是

$Q(\theta, \theta^{(i)}) = \sum_Z P(Z|Y, \theta^{(i)}) log P(Y,Z|\theta) \\$

此处需要区分一下$P(Z|Y, \theta^{(i)})$和$P(Z|\theta^{(i)})$，前者是我观测到$Y$的情况下，判断这个来自哪个分布，而后者则是在未观测到$Y$值的情况下就作出判断，所以后者其实对应的就是$\alpha_k$，前者我们给一个符号$\gamma_{jk}$，表示第$P(Z=k|Y=y_j, \theta^{(i)})$的情况下的估计。

因为对应到Q函数中$log$前的一项，所以需要这个东西。因为我们在做第$i+1$步更新的时候，$\theta^{(i)}$实际上就已经知道了，所以我们对于$P(Z=k|Y=y_j, \theta^{(i)})$也可以估计了。具体如下

$$\begin{align} P(Z=k|Y=y_j, \theta^{(i)}) & = \frac{P(Z=k, Y=y_j|\theta^{(i)})}{P(Y=y_j|\theta^{(i)})} \\ & = \frac{P(Z=k, Y=y_j|\theta^{(i)})}{\sum_Z P(Y=y_j, Z|\theta^{(i)})} \\ & = \frac{P(Z=k, Y=y_j|\theta^{(i)})}{\sum_k P(Y=y_j, Z=k|\theta^{(i)})} \\ & = \frac{\alpha_k \Phi(y_j|\theta_k^{(i)})}{\sum_t^K \alpha_t \Phi(y_j|\theta_t^{(i)})} \\ E(\gamma_{jk}) & = P(Z=k|Y=y_j, \theta^{(i)}) \\ \hat{\gamma_{jk}} & = E(\gamma_{jk}) \end{align}$$

到这里就是E步了，有了这个东西我们就可以写出完整的Q函数了

$$\begin{align} P(Y, Z|\theta) & = \prod_j^N P(Y=y_j, Z|\theta) \\ & = \prod_j^N \prod_k^K P(Y=y_j, Z=\gamma_{jk}) \\ & = \prod_j^N \prod_k^K (\alpha_k \Phi(y_j|\theta_k)) ^{\gamma_{jk}} \\ log P(Y, Z|\theta) & = \sum_k^K \sum_j^N \gamma_{jk} log \alpha_k + \gamma_{jk} log \Phi(y_j|\theta_k) \\ & = \sum_k^K n_k log \alpha_k + \sum_k^K \sum_j^N \gamma_{jk} log \Phi(y_j|\theta_k) \\ Q(\theta, \theta^{(i)}) & = \sum_k^K n_k log \alpha_k + \sum_k^K \sum_j^N E(\gamma_{jk}) log \Phi(y_j|\theta_k) \\ & = \sum_k^K n_k log \alpha_k + \sum_k^K \sum_j^N \hat{\gamma_{jk}} log \Phi(y_j|\theta_k) \end{align}$$

这里有个变量$\hat{\gamma_{jk}}$是需要在E步计算的，因为它是一个跟$\theta$相关的量，所以用$\theta^{(i)}$来代替，计算出来之后

$$\begin{align} \frac{\partial Q(\theta, \theta^{(i)})}{\partial \mu_k} & = \frac{\partial \sum_j^N \frac{(-\hat{\gamma_{jk}}(y_j-\mu_k)^2)}{2\sigma_k^2}}{\partial \mu_k} = \frac{\partial \sum_j^N \frac{\hat{\gamma_{jk}}(\mu_k - y_j)^2}{2\sigma_k^2}}{\partial \mu_k} = 0 \\ \sum_j^N \hat{\gamma_{jk}}y_j & = \sum_j^N \hat{\gamma_{jk}} \mu_k \\ \mu_k & = \frac{\sum_j^N \hat{\gamma_{jk}} y_j}{\sum_j^N \hat{\gamma_{jk}}} \end{align}$$

$$\begin{align} \frac{\partial Q(\theta, \theta^{(i)})}{\partial \sigma_k} & = \sum_j^N \hat{\gamma_{jk}} (-\frac{1}{\sigma_k} + \frac{(y_j - \mu_k)^2}{\sigma^3}) \\ & = \sum_j^N \hat{\gamma_{jk}} ((y_j - \mu_k)^2 - \sigma_k^2) = 0 \\ \sigma_k^2 & = \frac{\hat{\gamma_{jk}}(y_j-\mu_k)^2}{\sum_j^N \hat{\gamma_{jk}}} \end{align}$$

在对$\alpha_k$求导的时候有$\sum_k \n_k = 1$的限制，因此需要加上一个限制

$$\begin{align} \frac{\partial}{\partial \alpha_k} (Q(\theta, \theta^{(i)})) + \lambda(\sum_k \alpha_k - 1) &= \frac{n_k}{\alpha_k} + \lambda = 0 \\ \alpha_k & = \frac{n_k}{-\lambda} \\ \sum_k \alpha_k & = \sum_k \frac{n_k}{-\lambda} \\ \sum_k n_k & = -\lambda \end{align}$$