理解EM算法

---

---

EM（Expectation Maximization 期望最大化）算法是一种迭代算法，用于含有隐变量的概率模型参数的极大似然估计，或极大后验概率估计。其每次迭代由E、M两步构成。下面首先给出一般EM算法的求解过程（怎么做），然后结合一个例子来理解，然后讲为什么这么求解，即推导，最后讲述EM算法在高斯混合模型中的应用及小结。

---

一般用Y表示观测随机变量，Z表示隐随机变量，Y和Z在一起称为完全数据，Y称为不完全数据。由于含有隐变量，我们不能直接由P(Y|θ)=∑zP(Z|θ)P(Y|Z,θ)的最大似然估计来得到模型参数θ。EM算法就是在给定Y和Z的联合概率分布为P(Y,Z|θ)的情况下，通过迭代求解L(θ)=lnP(Y|θ)的极大似然估计来估算模型参数的算法。求解步骤如下：

1 选择一个初始的参数θold。

2 E Step 估计P(Z|Y,θold)。

3 M Step 估计θnew

θnew=argmaxθQ(θ,θold)

Q(θ,θold)=∑zP(Z|Y,θold)logP(Y,Z|θ)

4 检查是否到达停止迭代条件，一般是对较小的正数ε1,ε2，若满足：

‖θnew−θold‖<ε1或‖Q(θnew−θold)−Q(θold−θold)‖<ε2

则停止迭代，否则θold←θnew转到步骤2继续迭代。

---

下面结合一个《统计学习方法》中的例子来加深下理解：

例：假设有3枚硬币，分别记做A，B，C。这些硬币正面出现的概率分别是π,p和q。进行如下掷硬币实验：先掷硬币A，根据其结果选出硬币B或C，正面选B，反面选硬币C；然后投掷选重中的硬币，出现正面记作1，反面记作0；独立地重复n次（n=10)，结果为:

1111110000

假设只能观察到投掷硬币的结果，而不知其过程，问如何估计三硬币正面出现的概率，即三硬币的模型参数π,p和q。

解答：我们现在只可以看到硬币最终的结果1111110000，即观测变量Y，至于结果来自于B还是C无从得知，我们设隐变量Z来表示来自于哪个变量，令θ={π,p,q}。观测数据的似然函数可以表示为：

P(Y|θ)=∑zP(Z|θ)P(Y|Z,θ)=n∏j=1[πpyj(1−p)1−yj+(1−π)qyi(1−q)1−yj]

则模型参数θ的最大log似然估计为：

ˆθ=argmaxθlogP(Y|θ)

由于隐变量的存在，使得观测数据的最大似然函数里log里有带有加和(π..+(1−π)..)，导致上式是没有解析解的，只能通过迭代的方法求解，EM算法就是用于求解此类问题的一种迭代算法。

根据第二部分EM算法求解过程，假设第i次迭代参数的估计值为θ(i)=(π(i),p(i),q(i))，EM算法的第i+1次迭代如下：

E步：估计P(Z|Y,θ(i)),即在模型参数θ(i)下观测数据yj来自硬币B的概率

μ(i+1)=P(Z|Y,θ(i))=π(i)(p(i))yj(1−p(i))1−yjπ(i)(p(i))yj(1−p(i))1−yj+(1−π(i))(q(i))yj(1−q(i))1−yj

M步：估计θ(i+1),即：

θ(i+1)=argmaxθQ(θ,θi)

Q(θ,θ(i))=Ez[logP(Y,Z|θ)|Y,θ(i)]=∑zP(Z|Y,θ(i))logP(Y,Z|θ)=n∑j=1[μ(i+1)logπpyj(1−p)1−yj+(1−μ(i+1))log(1−π)qyi(1−q)1−yj]

μ(i+1)是E步得到常数，可通过分别对参数π,p和q求偏导使其为零来最大化上式，获得π,p和q的新的估计值：

π(i+1)=1nn∑j=1μ(i+1)j

p(i+1)=∑nj=1μ(i+1)jyj∑nj=1μ(i+1)j

q(i+1)=∑nj=1(1−μ(i+1)j)yj∑nj=1(1−μ(i+1)j)

在选定参数初始值θ(0)后，根据E步，M步循环迭代，直至满足迭代停止条件，即可得到参数θ的极大似然估计。

由上述EM计算过程和结合例子的应用，相信大家都会用EM算法解决问题了，即知道怎么做了，下面一节主要来讲述为什么这样做，即为什么这样就可以解决此类含有隐变量的最大似然估计。

EM算法的导出

---

面对一个含有隐变量的概率模型，目标是极大化观测数据（不完全数据），即极大化

L(θ)=logP(Y|θ)=log∑zP(Y,Z|θ)=log[∑zP(Z|θ)P(Y|Z,θ)]

极大化的主要困难是上式中含有未观测数据log里有包含和，EM算法则是通过迭代逐步近似极大化L(θ)的。假设在第i次迭代后θ的估计值是θ(i)，我们希望新估计的值θ能使L(θ)增加，即L(θ)>L(θ(i))，并逐步达到极大值，考虑两者的差：

L(θ)−L(θ(i))=log[∑zP(Z|θ)P(Y|Z,θ)]−logP(Y|θ(i))=log[∑zP(Z|Y,θ(i))P(Z|θ)P(Y|Z,θ)P(Z|Y,θ(i))]−logP(Y|θ(i))≥∑zP(Z|Y,θ(i))logP(Z|θ)P(Y|Z,θ)P(Z|Y,θ(i))−logP(Y|θ(i))=∑zP(Z|Y,θ(i))logP(Z|θ)P(Y|Z,θ)P(Y|θ(i))P(Z|Y,θ(i))

上式中的不等号是由Jensen不等式得到，下面对Jesen不等式做一个简单回顾。

Jensen不等式：如果f是凸函数，X是随机变量，则E[f(X)]≥f(EX)，特别地，如果f是严格凸函数，那么当E[f(X)]=f(EX)当且仅当P(x=E[X])=1，即X为常量。凹函数的性质和凸函数相反。

上式中f(X)为log(X)为凹函数，则E[f(X)]≤f(EX)，是上式中不等号的由来。

继续转回正题。令

B(θ,θ(i))≜L(θ(i))+∑zP(Z|Y,θ(i))logP(Z|θ)P(Y|Z,θ)P(Z|θ(i))P(Z|Y,θ(i))

L(θ)≥B(θ,θ(i))

即B(θ,θ(i))是L(θ)的下界，且L(θ(i))=B(θ(i),θ(i))

因此，可以使B(θ,θ(i))增大的θ也可以使L(θ)增大，为了使L(θ)有尽可能大的增长，选择θ(i+1)使得B(θ,θ(i))达到极大，即

θ(i+1)=argmaxθB(θ,θ(i))

省去对θ极大化而言是常熟的项，有

θ(i+1)=argmaxθ[L(θ(i))+∑zP(Z|Y,θ(i))logP(Z|θ)P(Y|Z,θ)P(Z|θ(i))P(Z|Y,θ(i))]=argmaxθ[∑zP(Z|Y,θ(i))logP(Z|θ)P(Y|Z,θ)]=argmaxθ[∑zP(Z|Y,θ(i))logP(Y,Z|θ)]=argmaxθQ(θ,θ(i))

推导完毕。

EM算法在高斯混合模型中的应用

---

EM算法可以用来估计高斯混合模型中的参数，例如：

假设观测数据y1,y2,...,yN由高斯混合模型生成

P(y|θ)=K∑k=1αkϕ(y|θk)

其中θ=(α1,α2,...,αk;θ1,θ2,...,θk)。

可以设想观测数据yj,j=1,2,...,N是这样产生，先根据概率αk选择第k个高斯分布ϕ(y|θk)，然后根据第k个模型的概率分布ϕ(y|θk)生成观测数据yj。此时观测数据yj是已知的，反映观测数据yj来自哪个分模型是未知的，看作隐变量。可以看出，可以用EM算法估计高斯混合模型（含有隐变量的概率模型参数）的参数θ。

1 首先确定E步，估计P(Z|Y,θ(i))，即在已知第i次迭代参数的情况下，观测数据yj来自计算分模型k的概率。

γ(i+1)jk=P(Z|Y,θ(i))=α(i)kϕ(y|θ(i)k)∑Kk=1α(i)kϕ(y|θ(i)k),j=1,2,...,N;k=1,2,...,K

2 M步，将新估计的Z的概率代进最大似然的公式，对待估计参数分别求偏导，以计算新一轮迭代模型参数（详细的推导不再赘述，感兴趣的可以自行推导）。

μ(i+1)k=∑Nj=1γ(i+1)jkyj∑Nj=1γ(i+1)jk,k=1,2,...,K

(σ2k)(i+1)=∑Nj=1γ(i+1)jk(yj−μ(i+1)k)2∑Nj=1γ(i+1)jk,k=1,2,...,K

α(i+1)k=∑Nj=1γ(i+1)jkN,k=1,2,...,K

重复E步，M步直至收敛。

---

1 EM算法是含有隐变量的概率模型极大似然估计或极大后验概率估计的迭代算法。含有隐变量的概率数据表示为P(Y,Z|θ)，Y表示观测变量，Z是隐变量，θ是模型参数。EM算法通过迭代求解观测数据的对数似然函数L(θ)=log(P|θ)的极大化来实现极大似然估计。每次迭代包含两步：

E步，求解P(Z|Y,θold)；

M步，估计θnew

θnew=argmaxθQ(θ,θold)

2 EM算法应用极其广泛，主要用于含有隐变量的概率模型的学习，但其对参数初始值比较敏感，而且不能保证收敛到全局最优。

---

[1] 统计学习方法.李航

[2] Pattern Recognition And Machine Learning.Christopher M. Bishop

[3] The EM Algorithm,JerryLead's Blog

[4] 三硬币问题