EM算法
EM算法基本思想
最大期望算法(Expectation-Maximization algorithm, EM),是一类通过迭代进行极大似然估计的优化算法,通常作为牛顿迭代法的替代,用于对包含隐变量或缺失数据的概率模型进行参数估计。
最大期望算法基本思想是经过两个步骤交替进行计算:
第一步是计算期望(E),利用对隐藏变量的现有估计值,计算其最大似然估计值
第二步是最大化(M),最大化在E步上求得的最大似然值来计算参数的值。
M步上找到的参数估计值被用于下一个E步计算中,这个过程不断交替进行。
EM算法推导
对于mmm个样本观察数据x=(x1,x2,...,xm)x=(x^{1},x^{2},...,x^{m})x=(x1,x2,...,xm),现在想找出样本的模型参数θ\thetaθ,其极大化模型分布的对数似然函数为:
θ=argmaxθ∑i=1mlogP(x(i);θ)
\theta = \mathop{\arg\max}_\theta\sum\limits_{i=1}^m logP(x^{(i)};\theta)
θ=argmaxθi=1∑mlogP(x(i);θ)
如果得到的观察数据有未观察到的隐含数据z=(z(1),z(2),...z(m))z=(z^{(1)},z^{(2)},...z^{(m)})z=(z(1),z(2),...z(m)),极大化模型分布的对数似然函数则为:
(a)θ=argmaxθ∑i=1mlogP(x(i);θ)=argmaxθ∑i=1mlog∑z(i)P(x(i),z(i);θ)
\theta =\mathop{\arg\max}_\theta\sum\limits_{i=1}^m logP(x^{(i)};\theta) = \mathop{\arg\max}_\theta\sum\limits_{i=1}^m log\sum\limits_{z^{(i)}}P(x^{(i)}, z^{(i)};\theta) \tag{a}
θ=argmaxθi=1∑mlogP(x(i);θ)=argmaxθi=1∑mlogz(i)∑P(x(i),z(i);θ)(a)
由于上式不能直接求出θ\thetaθ,采用缩放技巧:
∑i=1mlog∑z(i)P(x(i),z(i);θ)=∑i=1mlog∑z(i)Qi(z(i))P(x(i),z(i);θ)Qi(z(i))⩾∑i=1m∑z(i)Qi(z(i))logP(x(i),z(i);θ)Qi(z(i))
\sum\limits_{i=1}^m log\sum\limits_{z^{(i)}}P(x^{(i)}, z^{(i)};\theta) = \sum\limits_{i=1}^m log\sum\limits_{z^{(i)}}Q_i(z^{(i)})\frac{P(x^{(i)}, z^{(i)};\theta)}{Q_i(z^{(i)})} \\ \geqslant \sum\limits_{i=1}^m \sum\limits_{z^{(i)}}Q_i(z^{(i)})log\frac{P(x^{(i)}, z^{(i)};\theta)}{Q_i(z^{(i)})}
i=1∑mlogz(i)∑P(x(i),z(i);θ)=i=1∑mlogz(i)∑Qi(z(i))Qi(z(i))P(x(i),z(i);θ)⩾i=1∑mz(i)∑Qi(z(i))logQi(z(i))P(x(i),z(i);θ)
上式用到了Jensen不等式:
log∑jλjyj⩾∑jλjlogyj    ,λj⩾0,∑jλj=1
log\sum\limits_j\lambda_jy_j \geqslant \sum\limits_j\lambda_jlogy_j\;\;, \lambda_j \geqslant 0, \sum\limits_j\lambda_j =1
logj∑λjyj⩾j∑λjlogyj,λj⩾0,j∑λj=1
并且引入了一个未知的新分布Qi(z(i))Q_i(z^{(i)})Qi(z(i))。
此时,如果需要满足Jensen不等式中的等号,所以有:
P(x(i),z(i);θ)Qi(z(i))=c,c为常数
\frac{P(x^{(i)}, z^{(i)};\theta)}{Q_i(z^{(i)})} =c, c为常数
Qi(z(i))P(x(i),z(i);θ)=c,c为常数
由于Qi(z(i))Q_i(z^{(i)})Qi(z(i))是一个分布,所以满足
∑zQi(z(i))=1
\sum\limits_{z}Q_i(z^{(i)}) =1
z∑Qi(z(i))=1
综上,可得:
Qi(z(i))=P(x(i),z(i);θ)∑zP(x(i),z(i);θ)=P(x(i),z(i);θ)P(x(i);θ)=P(z(i)∣x(i);θ)
Q_i(z^{(i)}) = \frac{P(x^{(i)}, z^{(i)};\theta)}{\sum\limits_{z}P(x^{(i)}, z^{(i)};\theta)} = \frac{P(x^{(i)}, z^{(i)};\theta)}{P(x^{(i)};\theta)} = P( z^{(i)}|x^{(i)};\theta)
Qi(z(i))=z∑P(x(i),z(i);θ)P(x(i),z(i);θ)=P(x(i);θ)P(x(i),z(i);θ)=P(z(i)∣x(i);θ)
如果Qi(z(i))=P(z(i)∣x(i);θ)Q_i(z^{(i)}) = P( z^{(i)}|x^{(i)};\theta)Qi(z(i))=P(z(i)∣x(i);θ) ,则第(1)式是我们的包含隐藏数据的对数似然的一个下界。如果我们能极大化这个下界,则也在尝试极大化我们的对数似然。即我们需要最大化下式:
argmaxθ∑i=1m∑z(i)Qi(z(i))logP(x(i),z(i);θ)Qi(z(i))
\mathop{\arg\max}_\theta \sum\limits_{i=1}^m \sum\limits_{z^{(i)}}Q_i(z^{(i)})log\frac{P(x^{(i)}, z^{(i)};\theta)}{Q_i(z^{(i)})}
argmaxθi=1∑mz(i)∑Qi(z(i))logQi(z(i))P(x(i),z(i);θ)
简化得:
argmaxθ∑i=1m∑z(i)Qi(z(i))logP(x(i),z(i);θ)
\mathop{\arg\max}_\theta \sum\limits_{i=1}^m \sum\limits_{z^{(i)}}Q_i(z^{(i)})log{P(x^{(i)}, z^{(i)};\theta)}
argmaxθi=1∑mz(i)∑Qi(z(i))logP(x(i),z(i);θ)
以上即为EM算法的M步,∑z(i)Qi(z(i))logP(x(i),z(i);θ)\sum\limits_{z^{(i)}}Q_i(z^{(i)})log{P(x^{(i)}, z^{(i)};\theta)}z(i)∑Qi(z(i))logP(x(i),z(i);θ)可理解为logP(x(i),z(i);θ)logP(x^{(i)}, z^{(i)};\theta)logP(x(i),z(i);θ)基于条件概率分布Qi(z(i))Q_i(z^{(i)})Qi(z(i))的期望。以上即为EM算法中E步和M步的具体数学含义。
图解EM算法
考虑上一节中的(a)式,表达式中存在隐变量,直接找到参数估计比较困难,通过EM算法迭代求解下界的最大值到收敛为止。
图片中的紫色部分是我们的目标模型p(x∣θ)p(x|\theta)p(x∣θ),该模型复杂,难以求解析解,为了消除隐变量z(i)z^{(i)}z(i)的影响,我们可以选择一个不包含z(i)z^{(i)}z(i)的模型r(x∣θ)r(x|\theta)r(x∣θ),使其满足条件r(x∣θ)⩽p(x∣θ)r(x|\theta) \leqslant p(x|\theta)r(x∣θ)⩽p(x∣θ)。
求解步骤如下:
(1)选取θ1\theta_1θ1,使得r(x∣θ1)=p(x∣θ1)r(x|\theta_1) = p(x|\theta_1)r(x∣θ1)=p(x∣θ1),然后对此时的rrr求取最大值,得到极值点θ2\theta_2θ2,实现参数的更新。
(2)重复以上过程到收敛为止,在更新过程中始终满足r⩽pr \leqslant pr⩽p.
EM算法流程
输入:观察数据x=(x(1),x(2),...x(m))x=(x^{(1)},x^{(2)},...x^{(m)})x=(x(1),x(2),...x(m)),联合分布p(x,z;θ)p(x,z ;\theta)p(x,z;θ),条件分布p(z∣x;θ)p(z|x; \theta)p(z∣x;θ),最大迭代次数JJJ
1)随机初始化模型参数θ\thetaθ的初值θ0\theta^0θ0。
2)for j from 1 to Jfor \ j \ from \ 1 \ to \ Jfor j from 1 to J:
a) E步。计算联合分布的条件概率期望:
Qi(z(i))=P(z(i)∣x(i),θj)
Q_i(z^{(i)}) = P( z^{(i)}|x^{(i)}, \theta^{j})
Qi(z(i))=P(z(i)∣x(i),θj)
L(θ,θj)=∑i=1m∑z(i)P(z(i)∣x(i),θj)logP(x(i),z(i);θ) L(\theta, \theta^{j}) = \sum\limits_{i=1}^m\sum\limits_{z^{(i)}}P( z^{(i)}|x^{(i)}, \theta^{j})log{P(x^{(i)}, z^{(i)};\theta)} L(θ,θj)=i=1∑mz(i)∑P(z(i)∣x(i),θj)logP(x(i),z(i);θ)
b) M步。极大化L(θ,θj)L(\theta, \theta^{j})L(θ,θj),得到θj+1\theta^{j+1}θj+1:
θj+1=argmaxθL(θ,θj)
\theta^{j+1} = \mathop{\arg\max}_\theta L(\theta, \theta^{j})
θj+1=argmaxθL(θ,θj)
c) 如果θj+1\theta^{j+1}θj+1收敛,则算法结束。否则继续回到步骤a)进行E步迭代。
输出:模型参数θ\thetaθ。