EM算法

EM算法

这里作为我的博客备份，因为markdown解析各家标准并不能做到完全一致，特别是针对一些数学公式，如有排版问题，请访问原文EM算法 获得更好的排版体验

最近的工作中要借鉴到EM的思想，把EM又重新整理了下，总结此文。

基本原理

我们观察到的样本集，有时候仅仅是最终反映出来结果，即可观测的部分，而内在的一些隐含的因素往往不能被直接观察到，但算法不能忽略这些隐含状态，建模时要包含这些隐状态，而这样的模型有时候会比较复杂，导致直接优化原目标函数，不得不退而求其次，利用琴生不等式的性质，优化下确界，这就是EM的初衷了。

包含隐变量的loglikelihood的估计常常会用到EM，这种情况下，优化目标函数形式是：

$L = \sum_i {log{\sum_z {p(x,z;\theta)}}}$

其中x表示观测到的样本集，z表示隐变量，\(\theta\) 代表模型参数。注意到这个目标函数log里面包含求和符号\(\sum\)，这会导致模型优化上的复杂性，为了简化计算，利用琴生不等式有：

$L = \sum_i {log{\sum_z {q(z) \frac{p(x,z;\theta)}{q(z)}}}}$

$\quad \geqslant \sum_i \sum_z {q(z)} {log{ \frac{p(x,z;\theta)}{q(z)}}}$

EM算法不再优化原目标函数，而是优化更容易优化目标函数的下确界，在优化下确界时分为E步和M步，E步中计算目标函数，对隐变量进行分析，M步中利用E步对隐变量分析的结果，重新评估模型参数来达到最大化目标函数的目的，换句话说E步计算隐变量，M步计算模型参数，不停地重复E步和M步直到收敛。

这里还有一个问题没有解决，优化时引入了q(z)，但q(z)的公式到底应该怎么表示呢？注意到上述应用琴生不等式的步骤可以表述为：

$log({E_z{[\frac{p(x,z;\theta)}{q(z)}]})} \geqslant E_z{[log({\frac{p(x,z;\theta)}{q(z)}})] }$

上式如果等号成立，必有
$\frac{p(x,z;\theta)}{q(z)} = c$

$s.t. \sum_z q(z) = 1$
可以得到：

$q(z) = \frac{p(x,z;\theta)}{\sum_z p(x,z;\theta)}$

$\quad = \frac{p(x,z;\theta)}{p(x;\theta)}$

$\quad = {p(z|x;\theta)}$

EM实战

参考文献1举了一个很好掷硬币例子，可惜没有针对这个例子进行完整的理论分析，也没有相应的代码提供，这里我进行了一些完善，希望可以更加容易的理解EM算法。

考虑两个非标准的硬币，这里的非标准是指掷硬币时正反面概率不是相等的。对这两个硬币共进行了5组实验，每组实验开始前先选定一个硬币，然后用这个选定的硬币进行十次投掷，对结果进行记录，如下图所示：

可以看到，第一组，四组使用了硬币B，第二、三、五组使用硬币A，使用最大似然估计这两个硬币扔出正面的概率，A硬币共进行3组实验共投掷30次，其中24次正面，所以投掷正面的概率为0.8；B硬币共进行了两组实验共投掷20次，其中9次正面，所以投掷正面的概率为0.45。

上面的小例子应该是大部分人学习概率的第一个例子，下面对这个小例子进行一点点小小的变化，引入EM算法中重要的概念—–隐变量。仍然是上面的五组实验数据，现在假设不知道每组实验是由哪个硬币投掷的，只观测到每次投掷的正反面，如何利用最大似然求每个硬币投掷时正面朝上的概率呢？

在新的状况下，我们只观测到了一部分数据，即只有每次投掷结果是正面还是反面被观察到，另外一部分数据没有被观察到，即每组实验是采用的哪个硬币。这里可以发现，现在的情况和EM算法提出的动机一样了，所以也就意味着可以用EM算法来解决这个问题。为了使用EM算法来解决问题，首先引入一些变量：

\(i \in { 1,2,3,4,5 }\)，表示第i组实验

\(h_i\) 表示第i组实验中正面朝上的次数

\(n_i\) 表示第i组试验中反面朝上的次数

x_i 表示第i组实验的观测结果

z_i 表示第i组实验是有硬币\(z_i\)掷出的

\({p_A}^{(t)}\) 表示第t轮迭代A硬币正面朝上的概率

\({p_B}^{(t)}\) 表示第t轮迭代B硬币正面朝上的概率

有了这些变量，按照EM算法的流程，首先随机初始化\({p_A}^{t=0}\)和\({p_A}^{t=0}\)，然后E步、M步交替执行，知道收敛。

第t轮，E步：

\(\quad\quad\quad\quad\)计算loglikelihood L，并对隐变量进行分析：

$q^{(t)}(z_i=A) = \frac{p(x_i,z_i=A)}{\sum_{z_i} p(x_i,z_i)} = \frac{{p_A^{(t-1)}}^{h_i}{(1 - p_A^{(t-1)})}^{n_i}}{{p_A^{(t-1)}}^{h_i}{(1 - p_A^{(t-1)})}^{n_i} + {p_B^{(t-1)}}^{h_i}{(1 - p_B^{(t-1)})}^{n_i}}$

$q^{(t)}(z_i=B) = \frac{p(x_i,z_i=B)}{\sum_{z_i} p(x_i,z_i)} = \frac{{p_B^{(t-1)}}^{h_i}{(1 - p_B^{(t-1)})}^{n_i}}{{p_A^{(t-1)}}^{h_i}{(1 - p_A^{(t-1)})}^{n_i} + {p_B^{(t-1)}}^{h_i}{(1 - p_B^{(t-1)})}^{n_i}}$

第t轮，E步：

\(\quad\quad\quad\quad\)最大化loglikelihood L，求得使L最大化的\({P_A}^{(t)}\)和\({P_B}^{(t)}\)供下轮迭代使用。

$L = \sum_i log(\sum_z p(x_i,z_i))$

$\quad = \sum_i log(\sum_z q(z)\frac{p(x_i,z_i)}{q(z)})$

$\quad \geqslant \sum_i \sum_z q(z) log(\frac{p(x_i,z_i)}{q(z)})$

$\quad \quad\quad \quad\quad \quad\quad \quad\quad \quad \quad \quad\quad \quad\quad \quad\quad = \sum_i { q(z_i=A) log(\frac{p(x_i,z_i=A)}{q(z_i=A)}) + q(z_i=B) log(\frac{p(x_i,z_i=B)}{q(z_i=B)}) }$

$\quad \quad\quad \quad\quad \quad\quad \quad\quad \quad \quad \quad\quad \quad\quad \quad\quad = \sum_i { q(z_i=A) log(\frac{{p_A}^{h_i}{(1-p_A)}^{n_i}}{q(z_i=A)}) + q(z_i=B) log(\frac{{p_B}^{h_i}{(1-p_B)}^{n_i}}{q(z_i=B)}) }$

$\quad \quad\quad \quad\quad \quad\quad \quad\quad \quad \quad \quad\quad \quad\quad \quad\quad = \sum_i { q(z_i=A)[ h_ilog(p_A) + n_i log(1-p_A)] + q(z_i=B)[ h_ilog(p_B) + n_i log(1-p_B)] + C }$

其中\(q(z_i=A)\)和\(q(z_i=B)\)已经在E步求得，所以是常数

为求使得L最大的\(p_A\)和\(p_B\)，有：

$\frac{\partial L}{\partial p_A} = \sum_i \frac{q(z_i = A)h_i}{p_A} - \frac{q(z_i = A)n_i}{1-p_A} = 0$

$\sum_i {q(z_i = A) h_i (1-p_A)} - {q(z_i = A)n_ip_A} = 0$

$p_A = \frac{\sum_i q(z_i=A)h_i}{\sum_i q(z_i=A)h_i + q(z_i=A)n_i}$

同理：

$p_B = \frac{\sum_i q(z_i=B)h_i}{\sum_i q(z_i=B)h_i + q(z_i=B)n_i}$

整个迭代过程如下图所示：

想查看这个代码，可以从这里下载，迭代结果如下：

[epoch  0] pa=0.600000 pb=0.500000
[epoch  1] pa=0.713012 pb=0.581339
[epoch  2] pa=0.745292 pb=0.569256
[epoch  3] pa=0.768099 pb=0.549536
[epoch  4] pa=0.783165 pb=0.534617
[epoch  5] pa=0.791055 pb=0.526281
[epoch  6] pa=0.794532 pb=0.522390
[epoch  7] pa=0.795929 pb=0.520730
[epoch  8] pa=0.796466 pb=0.520047
[epoch  9] pa=0.796668 pb=0.519770
[epoch 10] pa=0.796744 pb=0.519659
ealy stop delta1=0.000028 delta2=0.000045

Ref

1）《What is the expectation maximization algorithm?》 Chuong B Do & Serafim Batzoglou

2）《The Expectation Maximization Algorithm》 frank dellaert

3）《CS229 Lecture notes-8 The EM algorithm》 Andrew Ng

联系我：