EM算法
这里作为我的博客备份,因为markdown解析各家标准并不能做到完全一致,特别是针对一些数学公式,如有排版问题,请访问原文EM算法 获得更好的排版体验
最近的工作中要借鉴到EM的思想,把EM又重新整理了下,总结此文。
基本原理
我们观察到的样本集,有时候仅仅是最终反映出来结果,即可观测的部分,而内在的一些隐含的因素往往不能被直接观察到,但算法不能忽略这些隐含状态,建模时要包含这些隐状态,而这样的模型有时候会比较复杂,导致直接优化原目标函数,不得不退而求其次,利用琴生不等式的性质,优化下确界,这就是EM的初衷了。
包含隐变量的loglikelihood的估计常常会用到EM,这种情况下,优化目标函数形式是:
L=∑ilog∑zp(x,z;θ)
其中x表示观测到的样本集,z表示隐变量,\(\theta\) 代表模型参数。注意到这个目标函数log里面包含求和符号\(\sum\),这会导致模型优化上的复杂性,为了简化计算,利用琴生不等式有:
L=∑ilog∑zq(z)p(x,z;θ)q(z)
⩾∑i∑zq(z)logp(x,z;θ)q(z)
EM算法不再优化原目标函数,而是优化更容易优化目标函数的下确界,在优化下确界时分为E步和M步,E步中计算目标函数,对隐变量进行分析,M步中利用E步对隐变量分析的结果,重新评估模型参数来达到最大化目标函数的目的,换句话说E步计算隐变量,M步计算模型参数,不停地重复E步和M步直到收敛。
这里还有一个问题没有解决,优化时引入了q(z),但q(z)的公式到底应该怎么表示呢?注意到上述应用琴生不等式的步骤可以表述为:
log(Ez[p(x,z;θ)q(z)])⩾Ez[log(p(x,z;θ)q(z))]
上式如果等号成立,必有
p(x,z;θ)q(z)=c
s.t.∑zq(z)=1
可以得到:
q(z)=p(x,z;θ)∑zp(x,z;θ)
=p(x,z;θ)p(x;θ)
=p(z|x;θ)
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)(zi=A)=p(xi,zi=A)∑zip(xi,zi)=p(t−1)Ahi(1−p(t−1)A)nip(t−1)Ahi(1−p(t−1)A)ni+p(t−1)Bhi(1−p(t−1)B)ni
q(t)(zi=B)=p(xi,zi=B)∑zip(xi,zi)=p(t−1)Bhi(1−p(t−1)B)nip(t−1)Ahi(1−p(t−1)A)ni+p(t−1)Bhi(1−p(t−1)B)ni
第t轮,E步:
\(\quad\quad\quad\quad\)最大化loglikelihood L,求得使L最大化的\({P_A}^{(t)}\)和\({P_B}^{(t)}\)供下轮迭代使用。
L=∑ilog(∑zp(xi,zi))
=∑ilog(∑zq(z)p(xi,zi)q(z))
⩾∑i∑zq(z)log(p(xi,zi)q(z))
=∑iq(zi=A)log(p(xi,zi=A)q(zi=A))+q(zi=B)log(p(xi,zi=B)q(zi=B))
=∑iq(zi=A)log(pAhi(1−pA)niq(zi=A))+q(zi=B)log(pBhi(1−pB)niq(zi=B))
=∑iq(zi=A)[hilog(pA)+nilog(1−pA)]+q(zi=B)[hilog(pB)+nilog(1−pB)]+C
其中\(q(z_i=A)\)和\(q(z_i=B)\)已经在E步求得,所以是常数
为求使得L最大的\(p_A\)和\(p_B\),有:
∂L∂pA=∑iq(zi=A)hipA−q(zi=A)ni1−pA=0
∑iq(zi=A)hi(1−pA)−q(zi=A)nipA=0
pA=∑iq(zi=A)hi∑iq(zi=A)hi+q(zi=A)ni
同理:
pB=∑iq(zi=B)hi∑iq(zi=B)hi+q(zi=B)ni
整个迭代过程如下图所示:
想查看这个代码,可以从这里下载,迭代结果如下:
[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
联系我: