EM算法

EM算法

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

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

基本原理

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

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

L=ilogzp(x,z;θ)

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

L=ilogzq(z)p(x,z;θ)q(z)

izq(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组实验,每组实验开始前先选定一个硬币,然后用这个选定的硬币进行十次投掷,对结果进行记录,如下图所示:

image

可以看到,第一组,四组使用了硬币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(t1)Ahi(1p(t1)A)nip(t1)Ahi(1p(t1)A)ni+p(t1)Bhi(1p(t1)B)ni

q(t)(zi=B)=p(xi,zi=B)zip(xi,zi)=p(t1)Bhi(1p(t1)B)nip(t1)Ahi(1p(t1)A)ni+p(t1)Bhi(1p(t1)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))

izq(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(1pA)niq(zi=A))+q(zi=B)log(pBhi(1pB)niq(zi=B))

=iq(zi=A)[hilog(pA)+nilog(1pA)]+q(zi=B)[hilog(pB)+nilog(1pB)]+C

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

为求使得L最大的\(p_A\)和\(p_B\),有:

LpA=iq(zi=A)hipAq(zi=A)ni1pA=0

iq(zi=A)hi(1pA)q(zi=A)nipA=0

pA=iq(zi=A)hiiq(zi=A)hi+q(zi=A)ni

同理:

pB=iq(zi=B)hiiq(zi=B)hi+q(zi=B)ni

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

image

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

[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

联系我:image

阅读更多
想对作者说点什么?

博主推荐

换一批

没有更多推荐了,返回首页