图解 Expectation Maximization 期望最大化 与应用例子

前言

网上对EM算法介绍已经很详尽，但是没看到比较详细的案例，理解起来有一些抽象。本文对EM的算法做一些总结，重点是介绍EM的案例，使得对该算法有一个直观的理解。

介绍

EM算法主要是针对存在隐变量的问题，即数据不完整的条件下去做参数估计。与之相反，当数据完整的时候，我们采用最大似然法就能解决问题。

似然函数

$$\begin{gathered} L(\theta )=\prod p(x_1,\cdots ,x_n;\theta )) \\ {\theta }^{\ast }=\underset{\theta \in \Theta }{argmax}L(\theta ) \\ ln(L(\theta ))=\sum _{i=1}^{n}ln(p(x_i;\theta )) \\ \frac{\partial ln(L(\theta ))}{\partial \theta }=0 \end{gathered}$$
求似然函数的极大意义在于，找出参数的最佳估计值。直观的感觉是，给定的样本和参数$\theta$越是符合，似然函数函数值自然越大。

Jensen不等式

这是一个容易理解的不等式。对于函数f，如果f满足，在其定义域内的所有实数x，满足$f^{{}^{\prime \prime }}(x)\ge 0$,则f为凸函数，且有$E(f(X))\ge f(E(X))$
由下图可以看到，f满足$(f(a)+f(b))/2\ge f((a+b)/2)$

用归纳法证明一下$E(f(X))\ge f(E(X))$
显然，当k=1，2时，命题成立
假定k时，命题成立
$$E(X)=\sum _{i=1}^k{p}_i{x}_{i}s.t.\sum _{i=1}^k{p}_i=1$$
现在求k+1时，命题是否成立
$$E(X)=\sum _{i=1}^{k+1}{p}_i{x}_{i}s.t.\sum _{i=1}^{k+1}{p}_i=1$$
令$s{p}_k={\sum }_{i=1}^k{p}_k\Rightarrow E(X)=p_{k+1}{x}_{k+1}+s{p}_k{\sum }_{i=1}^k\frac{p_i}{s{p}_k}{x}_i$
令$X_k={\sum }_{i=1}^k\frac{p_i}{s{p}_k}{x}_i$,代入上式得到$E(X)=p_{k+1}{x}_{k+1}+s{p}_k{X}_k$
由于$p_{k+1}+s{p}_k=1$，则有$E(f(X))=p_{k+1}f(x_{k+1})+s{p}_{k}f(X_k)\ge f(p_{k+1}{x}_{k+1}+s{p}_k{X}_k)=f(E(X))$
在f为严格凸函数时，上式只有在E(X)为常数时候才成立，即x是一个常数。
当f为凹函数时，上述不等式需要反号。我们知道最大似然函数ln函数是凹函数，所以上述需要取反。

EM 算法

假设完整数据为X ,Z
$X=[x_1,\cdots ,x_n]$为观测数据
$Z=[z_1,\cdots ,z_n]$为未观测的数据，即隐变量
则有
$$L(\theta )=\sum _{i=1}^{k}lnq(x_i;\theta )=\sum _{i=1}^{k}ln[\sum _{z}q(x_i,z_i;\theta )]$$
假设$Q_i$表示隐含变量Z的某种分布，$Q_i$满足
$$\sum _z{Q}_i(z)=1s.t.Q_i\ge 0$$
也就是说在$x=x_i$,条件下的Q_i(z)概率密度$q(z_i|x_i;\theta )$
$$L(\theta )=\sum _{i=1}^{k}ln[\sum _{z}q(x_i,z_i;\theta )]=\sum _{i=1}^{k}ln[\sum _z{Q}_i\frac{q(x_i,z_i;\theta )}{Q_i}]$$
这里的函数f 为 ln， 凹函数。令有变量$XX=\frac{q(x_i,z_i;\theta )}{Q_i}$
则有$E(XX)={\sum }_i{Q}_{i}XX$
由Jensen不等式可以得到
$$L(\theta )=\sum _{i=1}^{k}ln[\sum _z{Q}_{i}XX]\ge \sum _{i=1}^k\sum _z{Q}_{i}ln(XX)=J(Z,Q)$$
由于缺失了数据Z，$L(\theta )$本身并不存在一个解析解，直接用最大似然是做不到的。一切都只是为了得到一个比较合理的估计解。所以，利用Jensen不等式只是想求得一个更易于求解的下界表达式，再利用这个下界$J(Z,Q)$去逼近$L(\theta )$，直到得到一个最优解。
要想上式成立，唯有XX为常数才能做到，因此，为了能求出解，需要假设$XX=c$,则有
$$XX=\frac{q(x_i,z_i;\theta )}{Q_i}=C\Rightarrow \sum _{z}q(x_i,z_i;\theta )=q(x_i;\theta )=\sum _{i=1}^k{Q}_{i}C=C$$

EM算法是一种迭代逼近的算法，分为两步：

1. E步，固定$\theta$,计算$Q_i(z)$,建立$L(\theta )$的下界。我们知道z是无法观测的，但是有了给定的$\theta$就可以计算$Q_i(z)$
2. M步,得到$Q_i(z)$后，求${\theta }^{\ast }=\underset{\theta \in \Theta }{argmax}{\sum }_{i=1}^k{\sum }_z{Q}_{i}ln(XX)$
   不断重复1，2，直到$\theta$收敛。根据不同初始值，结果并不相同，因此，起始点是比较重要的。

EM例子

这个例子用的是李航的统计学习方法中的三个硬币例子。文中只给出了结果，并未对过程做详细推导，大概是作者觉得过于简单不需要去赘述吧。

令 θ = { π , p , q } \theta = \{ \pi ,p,q\} θ={π,p,q}，观测结果为y(同前面的x变量，为了和书本一致，这里使用y代替)，则有
$p(y|\theta )=\pi p^y{p}^{1-y}+(1-\pi )q^y{q}^{1-y}$
y的观测结果可能是来自硬币B，也可能来自硬币C，但是这一点我们是无法从y中获得的，因此，在这个问题中，我们将y的观测结果来源取决于A的投币结果，这个结果我们无法观测，因此作为隐变量。
$$\begin{gathered} Y={\left(Y_1,Y_2,\cdots ,Y_n\right)}^T,Z={\left(Z_1,Z_2,\cdots ,Z_n\right)}^T \\ P(Y|\theta )=\sum _{z}P(Z|\theta )P(Y|Z,\theta ) \\ P(Y|\theta )=\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] \\ \theta =(\pi ,p,q) \\ \hat{\theta }=\arg \underset{\theta }{\max }\log P(Y|\theta ) \end{gathered}$$

详细计算方法：
给定 ${\theta }^{(i)}=\left({\pi }^{(i)},p^{(i)},q^{(i)}\right).$
E步，计算$Q_i$
$\theta$已知， $Q_i$未知
前面提到过，$Q_i(z_i)=q(z_i|y_i;\theta )=q(y_i,z_i;\theta )/q(y_i;\theta )$
我们知道，前面的变量$XX=\frac{q(y_i,z_i;\theta )}{Q_i}=q(y_i,\theta )$,在给定$\theta$下，显然，该变量是常数项，使得 $L(\theta )=J(Z,Q)$

求后验概率，发生y_i时，正面为A的机率,也就是观测结果来自B的机率
$Q_{i1}(z_i=z_1(A正面))$即书中提到的下面的公式
$${\mu }^{(i+1)}=\frac{{\pi }^{(i)}{\left(p^{(i)}\right)}^{y_j}{\left(1-p^{(i)}\right)}^{1-y_j}}{{\pi }^{(i)}{\left(p^{(i)}\right)}^{y_j}{\left(1-p^{(i)}\right)}^{1-y_j}+\left(1-{\pi }^{(i)}\right){\left(q^{(i)}\right)}^{y_j}{\left(1-q^{(i)}\right)}^{1-y_j}}$$
可以得到
$Q_{i2}(z_i=z_2=A反面)=1-Q_{i1}(z_i=z_1)=1-{\mu }^{(i+1)}$

M步
$\theta$未知， $Q_i$已知
有了$Q_{i1}={\mu }_i,Q_{i2}$我们更新$$\begin{gathered} J(Q,Z)=Score=\sum _{i=1}^k\sum _z{Q}_i\ln (XX)=\sum _{i=1}^k(Q_{i1}ln[(q(y_i,z_1;\theta )/Q_{i1})]+Q_{i2}\ln (q(y_i,z_2;\theta )/Q_{i2})) \\ \frac{\partial (Q_{i1}ln[(q(y_i,z_1;\theta )/Q_{i1})]}{\partial \pi }=Q_{i1}{q}^{\prime }(y_i,z_1;\theta )/q(y_i,z_1;\theta )= \\ {Q}_{i1}{p}^y(1-p{)}^{1-y}/(\pi p^y(1-p{)}^{1-y})=Q_{i1}/\pi \\ \frac{\partial (Q_{i2}ln[(q(y_i,z_2;\theta )/Q_{i2})]}{\partial \pi }=Q_{i2}{q}^{\prime }(y_i,z_2;\theta )/q(y_i,z_2;\theta )= \\ -Q_{i2}{q}^y(1-q{)}^{1-y}/((1-\pi )q^y(1-q{)}^{1-y})=-Q_{i2}/(1-\pi ) \\ \frac{\partial Score}{\partial \pi }=\sum _i{Q}_{i1}/\pi -Q_{i2}/(1-\pi )=0\Rightarrow \\ \sum _i[(1-\pi )Q_{i1}-\pi Q_{i2}]=\sum _i(Q_{i1}-\pi )\Rightarrow \\ {\pi }^{(i+1)}=\frac{1}{n}\sum _{j=1}^n{\mu }_j^{(i)\cdot 1} \end{gathered}$$
同理可以得到
$$\begin{array}{l}{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\left(1-{\mu }_j^{(i+1)}\right)y_j}{{\sum }_{j=1}^n\left(1-{\mu }_j^{(i+1)}\right)}\end{array}$$

总结

用一张图总结EM的迭代是怎么工作的
通过前t-1次迭代，得到了${\theta }^t$,这时，$J(Z,Q)$对应的曲线如绿色曲线所式，这时候，E步构造新的$Q_i$参数，得到新的$J(Z,Q)$曲线（蓝色）,使得在${\theta }^t$位置,$J(Z,Q)=L(\theta )$。M步，对$J(Z,Q)$求导，求该曲线的极值位置作为${\theta }^{t+1}$，直到收敛。从图可以直观得看到，EM过程肯定是收敛得，因为极值保证了，新的J(Z,Q)一定大于等于老的J(Z,Q)。但是初始值若给的不合适，最终的结果未必理想。

参考文献

1. https://wenku.baidu.com/view/3396bb4d6294dd88d0d26bee.html
2. https://www.cnblogs.com/jerrylead/archive/2011/04/06/2006936.html