【西瓜书笔记】8. EM算法（上）

EM算法的引入

引入EM算法的原因：

概率模型有时候既含有观测变量，又含有隐变量或者潜在变量。如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法，或者贝叶斯估计法估计模型参数。但是当模型含有隐变量时，就不能简单地使用这些估计方法。EM算法就是含有隐变量的概率模型参数的极大似然估计法。

EM算法的例子

《统计学习方法》例9.1（三硬币模型）：

假设有3枚硬币，分别记作A，B，C。这些硬币正面出现的概率分别是 $\pi$ ，$p$ 和 $q$ 。进行如下掷硬币试验: 先掷硬币A，根据其结果选出硬币B或硬币C，正面选硬币B，反面选硬币C;然后掷选出的硬币，掷硬币的结果，出现正面记作1，出现反面记作0;独立地重复n次实验(这里，n=10)，观测结果如下
$$1,1,0,1,0,0,1,0,1,1$$
假设只能观测到掷硬币的结果，不能观测掷硬币的过程。问如何估计三硬币正面出现的概率，即三硬币模型的参数。

对于每一次实验可以进行如下建模：
$$P(y\mid \theta )=\sum _{z}P(y,z\mid \theta )=\sum _{z}P(z\mid \theta )P(y\mid z,\theta )$$
随机变量y是观测变量，表示一次试验观测的结果是1或0;随机变量z是隐变量，表示未观测到的掷硬币A的结果。这里其实利用了$P(A)={\sum }_{B}P(A,B)$，以及$P(A,B)=P(A)\cdot P(A|B)$。然后我们有
$$\begin{array}{rl}P(y\mid \theta )& =\sum _{z}P(y,z\mid \theta )=\sum _{z}P(z\mid \theta )P(y\mid z,\theta )\\ & =P(z=1\mid \theta )P(y\mid z=1,\theta )+P(z=0\mid \theta )P(y\mid z=0,\theta )\\ & =\{\begin{array}{ll}\pi p+(1-\pi )q,& \text{ if }y=1\\ \pi (1-p)+(1-\pi )(1-q),& \text{ if }y=0\end{array}\\ & =\pi p^y(1-p{)}^{1-y}+(1-\pi )q^y(1-q{)}^{1-y}\end{array}$$
这里$\theta =(\pi ,p,q)$是模型参数。将观测数据表示为$Y={\left(Y_1,Y_2,\dots ,Y_n\right)}^T$，未观测数据表示为$Z={\left(Z_1,Z_2,\dots ,Z_n\right)}^T$,则观测数据的似然函数为每次实验累乘的结果：
$$\begin{gathered} P(Y\mid \theta )=\sum _{Z}P(Z\mid \theta )P(Y\mid Z,\theta )=\prod _{j=1}^{n}P\left(y_j\mid \theta \right) \\ =\prod _{j=1}^n\left[\pi p^{y_j}(1-p{)}^{1-y_j}+(1-\pi )q^{y_j}(1-q{)}^{1-y_j}\right] \end{gathered}$$
考虑求模型参数$\theta =(\pi ,p,q)$的极大似然估计，也就是使用对数似然函数来进行参数估计可得：
$$\begin{array}{rl}\hat{\theta }& =\arg \underset{\theta }{\max }\ln P(Y\mid \theta )\\ & =\arg \underset{\theta }{\max }\ln \prod _{j=1}^n\left[\pi p^{y_j}(1-p{)}^{1-y_j}+(1-\pi )q^{y_j}(1-q{)}^{1-y_j}\right]\\ & =\arg \underset{\theta }{\max }\sum _{j=1}^n\ln \left[\pi p^{y_j}(1-p{)}^{1-y_j}+(1-\pi )q^{y_j}(1-q{)}^{1-y_j}\right]\end{array}$$
上式没有解析解，也就是没办法直接解出$\pi ,p,q$恰好等于某个常数a, b, c。因此我们只能用迭代的方法进行求解。

EM算法的导出

Jesen(琴生)不等式：

若$f$是凸函数，则：
$$f\left(t{x}_1+(1-t)x_2\right)\le tf\left(x_1\right)+(1-t)f\left(x_2\right)$$
其中，$t\in [0,1]$。同理，如果$f$是凹函数，则只需将上式中的$\le$换成$\ge$即可。

将上式中的$t$推广到$n$个变量，同样成立：
$$f\left(t_1{x}_1+t_2{x}_2+\dots +t_n{x}_n\right)\le t_{1}f\left(x_1\right)+t_{2}f\left(x_2\right)+\dots +t_{n}f\left(x_n\right)$$
其中，$t_1,t_2,\dots ,t_n\in [0,1],t_1+t_2+\dots +t_n=1$. 在概率论中常以以下形式出现
$$\phi (E[X])\le E[\phi (X)]$$
其中，$X$是随机变量，$\phi$是凸函数，$E[X]$表示$X$的期望。

我们面对一个含有隐变量的概率模型，目标是极大化观测数据Y关于参数θ的对数似然 函数，即极大化:
$$L(\theta )=\ln P(Y\mid \theta )=\ln \sum _{Z}P(Y,Z\mid \theta )=\ln \left(\sum _{Z}P(Y\mid Z,\theta )P(Z\mid \theta )\right)$$
注意到这一极大化的主要困难是上式中有未观测数据Z并有包含和(Z为离散型时)或者积分(Z为连续型时)的对数。EM算法采用的是通过迭代逐步近似极大化$L(\theta )$。假设在第i次迭代后$\theta$的估计值是${\theta }^{(i)}$，我们希望新的估计值$\theta$能够使$L(\theta )$增加，即$L(\theta )>L\left({\theta }^{(i)}\right)$并逐步达到极大值。为此，我们考虑两者的差:
$$\begin{array}{rl}L(\theta )-L\left({\theta }^{(i)}\right)& =\ln \left(\sum _{Z}P(Y\mid Z,\theta )P(Z\mid \theta )\right)-\ln P\left(Y\mid {\theta }^{(i)}\right)\\ & =\ln \left(\sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\frac{P(Y\mid Z,\theta )P(Z\mid \theta )}{P\left(Z\mid Y,{\theta }^{(i)}\right)}\right)-\ln P\left(Y\mid {\theta }^{(i)}\right)\end{array}$$
套用琴生不等式可有
$$\begin{array}{rl}& \ln \left(\sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\frac{P(Y\mid Z,\theta )P(Z\mid \theta )}{P\left(Z\mid Y,{\theta }^{(i)}\right)}\right)-\ln P\left(Y\mid {\theta }^{(i)}\right)\\ & \ge \sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\ln \frac{P(Y\mid Z,\theta )P(Z\mid \theta )}{P\left(Z\mid Y,{\theta }^{(i)}\right)}-\ln P\left(Y\mid {\theta }^{(i)}\right)\\ & =\sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\ln \frac{P(Y\mid Z,\theta )P(Z\mid \theta )}{P\left(Z\mid Y,{\theta }^{(i)}\right)}-1\cdot \ln P\left(Y\mid {\theta }^{(i)}\right)\end{array}$$
这里$P\left(Z\mid Y,{\theta }^{(i)}\right)$相当于式(7)中的$t_i$，对数函数相当于$f$。这不过这里是凹函数，所以不等式方向相反。又因为$1={\sum }_{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)$于是
$$\begin{array}{rl}L(\theta )-L\left({\theta }^{(i)}\right)& \ge \sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\ln \frac{P(Y\mid Z,\theta )P(Z\mid \theta )}{P\left(Z\mid Y,{\theta }^{(i)}\right)}-\sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\cdot \ln P\left(Y\mid {\theta }^{(i)}\right)\\ & =\sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\left(\ln \frac{P(Y\mid Z,\theta )P(Z\mid \theta )}{P\left(Z\mid Y,{\theta }^{(i)}\right)}-\ln P\left(Y\mid {\theta }^{(i)}\right)\right)\\ & =\sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\ln \frac{P(Y\mid Z,\theta )P(Z\mid \theta )}{P\left(Z\mid Y,{\theta }^{(i)}\right)P\left(Y\mid {\theta }^{(i)}\right)}\end{array}$$
所以
$$L(\theta )\ge L\left({\theta }^{(i)}\right)+\sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\ln \frac{P(Y\mid Z,\theta )P(Z\mid \theta )}{P\left(Z\mid Y,{\theta }^{(i)}\right)P\left(Y\mid {\theta }^{(i)}\right)}$$
令
$$B\left(\theta ,{\theta }^{(i)}\right)=L\left({\theta }^{(i)}\right)+\sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\ln \frac{P(Y\mid Z,\theta )P(Z\mid \theta )}{P\left(Z\mid Y,{\theta }^{(i)}\right)P\left(Y\mid {\theta }^{(i)}\right)}$$
则有
$$L(\theta )\ge B\left(\theta ,{\theta }^{(i)}\right)$$
在第一次迭代的时候，${\theta }_i$都是随机初始化的，通常初始化为$\pi =p=q=0.5$。现在我们不去极大化$L(\theta )$，因为前面说过这很困难。我们转而去极大化它的下界$B\left(\theta ,{\theta }^{(i)}\right)$，得到一个新的$\theta$，然后把这个新的$\theta$代入到$L(\theta )$，看是不是能使得$L(\theta )$变大。也就说$B\left(\theta ,{\theta }^{(i)}\right)$是$L(\theta )$的一个下界，此时若设${\theta }^{(i+1)}$使得$B\left(\theta ,{\theta }^{(i)}\right)$达到极大(不是最大)，也即
$$B\left({\theta }^{(i+1)},{\theta }^{(i)}\right)\ge B\left({\theta }^{(i)},{\theta }^{(i)}\right)$$
进一步可得
$$\begin{gathered} L\left({\theta }^{(i+1)}\right)\ge B\left({\theta }^{(i+1)},{\theta }^{(i)}\right)\ge B\left({\theta }^{(i)},{\theta }^{(i)}\right)=L\left({\theta }^{(i)}\right) \\ \Rightarrow L\left({\theta }^{(i+1)}\right)\ge L\left({\theta }^{(i)}\right) \end{gathered}$$
这里注意:
$$\begin{gathered} B\left({\theta }^{(i)},{\theta }^{(i)}\right)=L\left({\theta }^{(i)}\right)+\sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\ln \frac{P(Y\mid Z,{\theta }^{(i)})P(Z\mid {\theta }^{(i)})}{P\left(Z\mid Y,{\theta }^{(i)}\right)P\left(Y\mid {\theta }^{(i)}\right)} \\ =L\left({\theta }^{(i)}\right)+\sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\ln \frac{P(Y,Z\mid {\theta }^{(i)})}{P(Y,Z\mid {\theta }^{(i)})} \\ =L({\theta }^{(i)}) \end{gathered}$$
所以，任何可以使$B(\theta ,{\theta }^{(i)})$增大的$\theta$，也可使$L(\theta )$增大，于是问题转化为了求解能使得$B\left(\theta ,{\theta }^{(i)}\right)$达到极大的${\theta }^{(i+1)}$，即
$$\begin{array}{rl}{\theta }^{(i+1)}& =\underset{\theta }{\mathrm{argmax}}B\left(\theta ,{\theta }^{(i)}\right)\\ & =\underset{\theta }{\mathrm{argmax}}\left(L\left({\theta }^{(i)}\right)+\sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\ln \frac{P(Y\mid Z,\theta )P(Z\mid \theta )}{P\left(Z\mid Y,{\theta }^{(i)}\right)P\left(Y\mid {\theta }^{(i)}\right)}\right)\\ & =\underset{\theta }{\mathrm{argmax}}\left(\sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\ln (P(Y\mid Z,\theta )P(Z\mid \theta ))\right)\\ & =\underset{\theta }{\mathrm{argmax}}\left(\sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\ln P(Y,Z\mid \theta )\right)\\ & =\underset{\theta }{\mathrm{argmax}}Q\left(\theta ,{\theta }^{(i)}\right)\end{array}$$
到此即完成了EM算法的一次迭代，求出的${\theta }^{(i+1)}$作为下一次迭代的初始${\theta }^{(i)}$。综上可以总结出EM算法的“E步”和“M步”分别为:

• E步：计算完全数据的对数似然函数$\ln P(Y,Z\mid \theta )$关于在给定观测数据$Y$和当前参数$\theta$下对未观测数据$Z$的条件概率分布$P\left(Z\mid Y,{\theta }^{(i)}\right)$的期望，也就是Q函数
  $$Q\left(\theta ,{\theta }^{(i)}\right)=E_Z\left[\ln P(Y,Z\mid \theta )\mid Y,{\theta }^{(i)}\right]=\sum _{Z}P\left(Z\mid Y,{\theta }^{(i)}\right)\ln P(Y,Z\mid \theta )$$

• M步：求使得Q函数到达极大的${\theta }^{(i+1)}$.

(souce:https://www.borealisai.com/en/blog/tutorial-5-variational-auto-encoders/)