EM算法指的是最大期望算法(Expectation Maximization Algorithm),是一种迭代算法,用于含有隐变量的概率参数模型的最大似然或极大后验概率估计。
EM算法的引入
例(三硬币模型)
假设有三枚硬币A,B,C,每个硬币正面出现的概率分别为\(\pi,p,q\)。进行如下的掷硬币实验,先掷硬币A,如果正面朝上则选硬币B,反面朝上选C;然后投掷选择的硬币,出现的结果正面为1,反面为0;独立地重复n次实验,这里n=10,结果如下:
假设智能观测到掷硬币的结果,不能观测掷硬币的过程。问如果估计三硬币正面出现的概率,即三硬币模型的参数?
- 解答:
针对某个输出值y,它在参数\(\theta=(\pi,p,q)\)下的概率分布为:
将观测数据表示为\(Y=(Y_1,Y_2,\cdots,Y_n)^T\),未观测数据表示为\(Z=(Z_1,Z_2,\cdots,Z_n)^T\),则,观测数据的对数似然函数为:
或者可以写成:
直接连乘的似然函数求导太复杂,所以一般用极大似然估计都会转化成对数似然函数。
本题的目标就变成了求解参数\(\theta\)的对数极大似然估计,即:
这个问题没有解析解,只有通过迭代的方式求解。
- 接下来的部分将介绍如何求解
预备知识
凹凸函数
定义:如果函数总是位于任何一条弦的下方,则该函数是凸的;如果函数总是位于任何一条弦的上方,则函数是凹的。
用数学语言表述如下:
如果函数f在某个区间上存在非负(正)的二阶导数,即\(f''(x)\ge 0\),则f为该区间的凸函数。当x是向量时,如果其Hessian矩阵H是半正定的,那么f是凸函数。如果\(f''(x)\ge 0\)或者\(H>0\),那么f是严格凸函数。
凸函数的例子有\(x^2,|x|,e^x,x \log x(x\ge 0)\);凹函数的例子有log x,\(\sqrt{x}(x\ge 0)\)等。
Jensen不等式:
如果f是凸函数,X是随机变量,那么
特别地,如果f是严格凸的,等式成立当且仅当P(X=E[X])=1,即X是常量。
Jensen不等式应用于凹函数时,不等号方向反向。
\(p(x;\theta)\)和\(p(x|\theta)\)
- 定义:
假定一个关于参数\(\theta\),具有离散型概率分布P的随机变量X,则在给定X的输出x时,参数的似然函数可表示为:
其中,\(p(x)\)表示X取x时的概率。上式常常写为\(P(X=x|\theta)\)或者\(P(X=x;\theta)\).需要注意的是,此处并非条件概率,因为\(\theta\)并不(总)是随机变量。
\(p(x;\theta)\)写分号时,表示待估参数是固定的(只是当前未知),加分号带参数是为了说明这里有个\(\theta\)的参数,\(p(x;\theta)\)的意思就是随机变量X=x在参数\(\theta\)下的概率。本节中,这两种记号有着相同的意思。
EM算法
- 一般地,用Y表示观测随机变量的数据,Z表示隐随机变量的数据,Y和Z连在一起称为完全数据,观测数据Y又称为不完全数据。不完全数据Y的似然函数是\(\log P(Y|\theta)\),完全数据的对数似然函数为\(\log P(Y,Z|\theta)\)。
- EM算法通过迭代求\(L(\theta)=\log P(Y|\theta)\)的极大似然估计。
给定训练样本为\({\{y_1,y_2,\cdots,y_n}\}\),样例之间相互独立,我们想找到每个样例的隐含类别z,能使得p(y,z)最大。(这里\(z_i=j\)可以看成样本\(y_i\)被划分成j类) 记模型参数为\(\theta\)。
训练样本的对数极大似然函数如下:
第一步是对极大似然函数取对数,第二步是对每个样例的每个可能的类别j求联合分布概率和,但是直接求\(\theta\)一般比较困难,因为有隐变量z的存在,但是一般确定了z之后,求解就很容易了。
思想
EM是一种解决存在隐含变量优化问题的有效方法。既然不能直接最大化\(L(\theta)\),我们可以不断建立\(L(\theta)\)的下界(E步),然后不断优化下界(M步)。
分析
事实上,\(z_i\)也是个随机变量,对于每一个\(y_i\),都有多种分类的情况。设第i个样本\(y_i\)在z上的概率分布为\(Q_i(z_i)\),即\(Q_i(z_i=j)\)表示\(y_i\)被划分到类j的概率,因此有\(\sum_{z_i}Q_i(z_i)=1\)。
(如果z是连续性的,那么\(Q_i(z_i)\)是概率密度函数,需要将求和符号换做积分符号.)
比如要将班上的同学聚类,假设隐变量z是身高,那么就是连续的高斯分布。如果隐变量是男女,那么就是伯努利分布了。
进行变换得到以下公式:
这里,由(1)到(2)式比较直接,分子分母同乘以一个相等的;(2)到(3)利用了Jensen不等式,log(x)是凹函数,这里x=\(\dfrac{p(y_i,z_i;\theta)}{Q_i(z_i)}\).
这个过程可以看作是对\(L(\theta)\)求了下界。假设\(\theta\)已经给定,那么\(L(\theta)\)的值就取决于\(Q_i(z_i)\)和\(p(y_i,z_i;\theta)\)了。我们通过调整这两个概率使下界不断上升,逼近\(L(\theta)\)的值。
由Jensen不等式可知,要想等式成立,随机变量x=\(\dfrac{p(y_i,z_i;\theta)}{Q_i(z_i)}\)应为常数。故我们有:
\(\dfrac{p(y_i,z_i;\theta)}{Q_i(z_i)}=c,\) c为常数
已知\(\sum_{z_i}Q_i(z_i)=1\),我们对上式进行变换并关于\(z_i\)求和得:
代回原式解得:
至此,我们推出了在固定\(\theta\)之后,\(Q_i(z_i)\)的计算公式就是后验概率,解决了\(Q_i(z_i)\)的选择问题,建立了\(L(\theta)\)的下界。接下来的M步,就是在给定\(Q_i(z_i)\)后,调整\(\theta\),去极大化\(L(\theta)\)的下界。
小结
一般的EM算法步骤如下:
循环重复直至收敛{
(E步)对于每一个i,在给定\(\theta\)的条件下计算
也等价于最大化$$E_z[\log\dfrac{p(y_i,z_i;\theta)}{Q_i(z_i)}]$$
(M步)计算
这里,我们抹掉了对于\(Q_i(z_i)\)而言是常数项的值
}
- 如果我们定义
从上面的推导我们知道,\(L(\theta)\ge J(Q,\theta)\),EM可以看作是J的坐标上升法,E步固定\(\theta\),优化Q,M步固定Q优化\(\theta\).
另一种推导方法
我们面对一个含有隐变量的概率模型,目标是极大化观测数据(不完全数据)Y关于参数\(\theta\)的对数似然函数,即极大化:
这一极大化的主要困难是上式中包含有未观测数据并包含和的对数。
事实上,EM算法是通过迭代逐步近似极大化\(L(\theta)\)的。假设在第t次迭代后\(\theta\)的估计值是\(\theta^{(t)}\),我们希望新的估计值能使\(L(\theta)\)增加,即
\(L(\theta)>L(\theta^{(t)})\),并逐步达到最大值,为此,考虑两者的差:
利用Jensen不等式,得到其下界:
令$$B(\theta,\theta^{(t)}) \hat{=} L(\theta{(t)})+\sum_ZP(Z|Y,\theta) \log\dfrac{P(Z|\theta)P(Y|Z,\theta))}{P(Z|Y,\theta{(t)})P(Y|\theta)}$$
则 $$L(\theta)\ge B(\theta,\theta^{(t)})$$
即函数\(B(\theta,\theta^{(t)})\)是\(L(\theta)\)的一个下界,且由上可知:
因此,任何可以使\(B(\theta,\theta^{(t)})\)增大的\(\theta\),也可以使\(L(\theta)\)增大。
为了使\(L(\theta)\)有尽可能大的增长,选择\(\theta^{(t+1)}\)使得\(B(\theta,\theta^{(t)})\)达到极大,即
现在求\(\theta^{(t+1)}\)的表达式,省去对\(\theta\)的极大化而言是常数的项,过程如下:
上式等价于EM算法的一次迭代,即求Q函数及其最大化。EM算法是通过不断求解下界的极大化逼近求解对数似然函数极大化的算法。
算法步骤
输入: 观测便利数据Y,隐变量数据Z,联合分布\(P(Y,Z|\theta)\),条件分布\(P(Z|Y,\theta)\);
输出:模型参数\(\theta\)
(1)选择参数的初值\(\theta^{(0)}\),开始迭代;
(2)E步:记\(\theta^{(t)}\)为第t次迭代参数\(\theta\)的估计值,在第t+1次迭代的E步,计算:
这里\(P(Z|Y,\theta^{(t)})\)是在给定观测数据Y和当前的参数估计\(\theta^{(t)}\)下隐变量数据Z的条件概率分布;
(3)M步: 求使得\(Q(\theta,\theta^{(t)})\)极大化的\(\theta\),确定第t+1次迭代的参数的估计值\(\theta^{(t+1)}\):
定义(Q函数)
完全数据的对数似然函数\(P(Y,Z|\theta)\)关于在给定观测数据和当前参数\(\theta^{(t)}\)下对未观测数据Z的条件概率分布\(P(Z|Y,\theta^{(t)}\)的期望称为Q函数,即
EM收敛性证明
假定\(\theta^{(t)}\)和\(\theta^{(t+1)}\)是EM第t次和第t+1次迭代后的结果。如果我们证明了\(l(\theta^{(t)})\le l(\theta^{(t+1)})\),也就是说极大似然估计单调增加,那么我们最终会到达极大似然估计的最大值。
证明如下:
选定\(\theta^{(t)}\)后,我们得到E步
这一步保证了在给定\(\theta^{(t)}\)时,Jensen不等式成立,也就是
然后进行M步,固定\(Q_i^{(t)}(z_i)\),并将\(\theta^{(t)}\)视作变量,对上面的\(L(\theta^{(t)})\)求导后,得到\(\theta^{(t+1)}\),这样经过一些推导会有以下式子成立:
上面(4)式由(3)式可得,当固定了\(\theta\)之后,不等式成立;
(5)式利用了M步的定义,
(6)式是之前等式的结果。
综上,我们证明了\(L(\theta)\)会单调增加,又因为\(P(Y|\theta)\)是有界的,因此EM算法是收敛的。
另一种证明方法
- 定理:
设\(P(Y|\theta)\)为观测数据的似然函数,\(\theta^{(t)}(t=1,2,\cdots,n)\)为EM算法得到的参数估计序列,\(P(Y|\theta^{(t)})\)为对应的似然函数序列,则\(P(Y|\theta^{(t)})\)是单调增的,即:
去对数有:
已知$$Q(\theta,\theta^{(t)})=\sum_Z P(Z|Y,\theta^{(t)})\log P(Y,Z|\theta)$$
令$$H(\theta,\theta^{(t)})=\sum_Z P(Z|Y,\theta^{(t)})\log P(Z|Y,\theta)$$
则$$\log P(Y|\theta)=Q(\theta,\theta{(t)})-H(\theta,\theta)$$
在上式中分别令\(\theta\)为\(\theta^{(t)}\)和\(\theta^{(t+1)}\)并相减,则有:
这里的不等号由Jensen不等式得到。故等式的右端非负,定理得证。
综上,\(L(\theta)\)单调递增,且\(P(Y|\theta)\)有界,故\(P(Y|\theta)\)有上界,则
\(L(\theta^{(t)})=\log P(Y|\theta^{(t)})\)收敛到某一值L*,EM算法收敛性得证。
下证明上式右端大于等于0:
已知\(\theta^{(t+1)}\)使得\(Q(\theta,\theta^{(t)})\)达到最大,因此有:
续(三硬币模型)
由第一部分推导可知:
\(P(y_i;\theta)=\pi p^{y_i}(1-p)^{1-y_i}+(1-\pi)q^{y_i}(1-q)^{1-y_i}\)
计算\(Q_i(z_i)\):
同理,\(1-\mu_i=Q_i(z_i=0)=\dfrac{(1-\pi)q^{y_i}(1-q)^{1-y_i}}{\pi p^{y_i}(1-p)^{1-y_i}+(1-\pi)q^{y_i}(1-q)^{1-y_i}}\)
这里\(z_i=1\)表示第一次掷A硬币得到正面的结果;\(z_i=0\)表示第一次掷A硬币得到反面的结果。
由EM算法步骤:
- (1)初始化: 选取模型的初值:\(\theta^{(0)}=(\pi^{(0)},p^{(0)},q^{(0)})\).
- (2)E步:我们计算在当前迭代的模型参数下(也就是第t次迭代后),观测数据来自硬币B的概率:
观测数据来自硬币C的概率:
此时Q函数的形式为:
- (3)M步: 计算$$\theta^{(t+1)}=\arg\max\limits_\theta Q(\theta,\theta^{(t)})$$
对$ Q(\theta,\theta^{(t)})\(关于三个参数求偏导数,获得参数\)\pi,p,q$的新的估计值。
令上面三个偏导函数等于0,我们得到参数的估计值如下:
- (4)重复第(2)和(3)步,直到收敛。
注意事项
步骤(1) 参数的初值可以任意选择,但需注意EM算法对初值是敏感的;
步骤(2) E步求\(Q(\theta,\theta^{(t)})\).Q函数中Z是未观测数据,Y是观测数据。函数中第一个变量表示要极大化的参数,第二个变量表示参数的当前估计值。每次迭代实际在求Q函数及其极大;
步骤(3)M步求\(Q(\theta,\theta^{(t)})\)的极大化,得到\(\theta^{(t+1)}\),完成一次迭代\(\theta^{(t)}\rightarrow\theta^{(t+1)}\).每次迭代使似然函数增大或达到局部极值;
步骤(4)给出停止迭代的条件,一般是对较小的正数\(\varepsilon_1,\varepsilon_2\),若满足:
\(\qquad \lVert \theta^{(t+1)}-\theta^{(t)}\rVert\le\varepsilon_1\),或,\(\lVert Q(\theta,\theta^{(t+1)})-Q(\theta,\theta^{(t)})\rVert\le\varepsilon_2\)
则停止迭代。