[!IMPORTANT]
参与组队学习的同学须知:
本章学习时间:3天
本章配套视频教程:https://www.bilibili.com/video/BV1Mh411e7VU?p=11
第7章 贝叶斯分类器
本章是从概率框架下的贝叶斯视角给出机器学习问题的建模方法,不同于前几章着重于算法具体实现,本章的理论性会更强。朴素贝叶斯算法常用于文本分类,例如用于广告邮件检测,贝叶斯网和EM算法均属于概率图模型的范畴,因此可合并至第14章一起学习。
7.1 贝叶斯决策论
7.1.1 式(7.5)的推导
由式(7.1)和式(7.4)可得
R(ci∣x)=1∗P(c1∣x)+...+1∗P(ci−1∣x)+0∗P(ci∣x)+1∗P(ci+1∣x)+...+1∗P(cN∣x)
又∑j=1NP(cj∣x)=1,则
R(ci∣x)=1−P(ci∣x)
此即式(7.5)。
7.1.2 式(7.6)的推导
将式(7.5)代入式(7.3)即可推得此式
7.1.3 判别式模型与生成式模型
对于判别式模型来说,就是在已知x的条件下判别其类别标记c,即求后验概率P(c∣x),前几章介绍的模型都属于判别式模型的范畴,尤其是对数几率回归最为直接明了,式(3.23)和式(3.24)直接就是后验概率的形式。
对于生成式模型来说,理解起来比较抽象,但是可通过思考以下两个问题来理解。
(1)对于数据集来说,其中的样本是如何生成的?通常假设数据集中的样本服从独立同分布,即每个样本都是按照联合概率分布P(x,c)采样而得,也可以描述为根据P(x,c)生成的。
(2)若已知样本x和联合概率分布P(x,c),如何预测类别呢?若样本x和联合概率分布P(x,c)已知,则可以分别求出x属于各个类别的概率,即P(x,c1),P(x,c2),...,P(x,cN),然后选择概率最大的类别作为样本x的预测结果。
因此,之所以称为"生成式"模型,是因为所求的概率P(x,c)是生成样本x的概率。
7.2 极大似然估计
7.2.1 式(7.12)和(7.13)的推导
根据式(7.11)和式(7.10)可知参数求解式为
θ^c=θcargmaxLL(θc)=θcargmin−LL(θc)=θcargmin−x∈Dc∑logP(x∣θc)
由"西瓜书"上下文可知,此时假设概率密度函数p(x∣c)∼N(μc,σc2),其等价于假设
P(x∣θc)=P(x∣μc,σc2)=(2π)d∣Σc∣1exp(−21(x−μc)TΣc−1(x−μc))
其中,d表示x的维数,Σc=σc2为对称正定协方差矩阵,∣Σc∣表示Σc的行列式。将其代入参数求解式可得
(μ^c,Σ^c)=(μc,Σc)argmin−x∈Dc∑log[(2π)d∣Σc∣1exp(−21(x−μc)TΣc−1(x−μc))]=(μc,Σc)argmin−x∈Dc∑[−2dlog(2π)−21log∣Σc∣−21(x−μc)TΣc−1(x−μc)]=(μc,Σc)argminx∈Dc∑[2dlog(2π)+21log∣Σc∣+21(x−μc)TΣc−1(x−μc)]=(μc,Σc)argminx∈Dc∑[21log∣Σc∣+21(x−μc)TΣc−1(x−μc)]
假设此时数据集Dc中的样本个数为n,即∣Dc∣=n,则上式可以改写为
(μ^c,Σ^c)=(μc,Σc)argmini=1∑n[21log∣Σc∣+21(xi−μc)TΣc−1(xi−μc)]=(μc,Σc)argmin2nlog∣Σc∣+i=1∑n21(xi−μc)TΣc−1(xi−μc)
为了便于分别求解μ^c和Σ^c,在这里我们根据式xTAx=tr(AxxT),xˉ=n1∑i=1nxi将上式中的最后一项作如下恒等变形:
==========i=1∑n21(xi−μc)TΣc−1(xi−μc)21tr[Σc−1i=1∑n(xi−μc)(xi−μc)T]21tr[Σc−1i=1∑n(xixiT−xiμcT−μcxiT+μcμcT)]21tr[Σc−1(i=1∑nxixiT−nxˉμcT−nμcxˉT+nμcμcT)]21tr[Σc−1(i=1∑nxixiT−2nxˉμcT+nμcμcT+2nxˉxˉT−2nxˉxˉT)]21tr[Σc−1((i=1∑nxixiT−2nxˉxˉT+nxˉxˉT)+(nμcμcT−2nxˉμcT+nxˉxˉT))]21tr[Σc−1(i=1∑n(xi−xˉ)(xi−xˉ)T+i=1∑n(μc−xˉ)(μc−xˉ)T)]21tr[Σc−1i=1∑n(xi−xˉ)(xi−xˉ)T]+21tr[Σc−1i=1∑n(μc−xˉ)(μc−xˉ)T]21tr[Σc−1i=1∑n(xi−xˉ)(xi−xˉ)T]+21tr[n⋅Σc−1(μc−xˉ)(μc−xˉ)T]21tr[Σc−1i=1∑n(xi−xˉ)(xi−xˉ)T]+2ntr[Σc−1(μc−xˉ)(μc−xˉ)T]21tr[Σc−1i=1∑n(xi−xˉ)(xi−xˉ)T]+2n(μc−xˉ)TΣc−1(μc−xˉ)
所以
(μ^c,Σ^c)=(μc,Σc)argmin2nlog∣Σc∣+21tr[Σc−1i=1∑n(xi−xˉ)(xi−xˉ)T]+2n(μc−xˉ)TΣc−1(μc−xˉ)
观察上式可知,由于此时Σc−1和Σc一样均为正定矩阵,所以当μc−xˉ=0时,上式最后一项为正定二次型。根据正定二次型的性质可知,此时上式最后一项的取值仅与μc−xˉ相关,并有当且仅当μc−xˉ=0时,上式最后一项取最小值0,此时可以解得
μ^c=xˉ=n1i=1∑nxi
将求解出来的μ^c代回参数求解式可得新的参数求解式,有
Σ^c=Σcargmin2nlog∣Σc∣+21tr[Σc−1i=1∑n(xi−xˉ)(xi−xˉ)T]
此时的参数求解式是仅与Σc相关的函数。
为了求解Σ^c,在这里我们不加证明地给出一个引理:设B为p阶正定矩阵,n>0为实数,在对所有p阶正定矩阵Σ有
2nlog∣Σ∣+21tr[Σ−1B]≥2nlog∣B∣+2pn(1−logn)
当且仅当Σ=n1B时等号成立。(引理的证明可搜索张伟平老师的"多元正态分布参数的估计和数据的清洁与变换"课件)
根据此引理可知,当且仅当Σc=n1∑i=1n(xi−xˉ)(xi−xˉ)T
时,上述参数求解式中argmin后面的式子取到最小值,那么此时的Σc即我们想要求解的Σ^c。
7.3 朴素贝叶斯分类器
7.3.1 式(7.16)和式(7.17)的解释
该式是基于大数定律的频率近似概率的思路,而该思路的本质仍然是极大似然估计,下面举例说明。以掷硬币为例,假设投掷硬币5次,结果依次是正面、正面、反面、正面、反面,试基于此观察结果估计硬币正面朝上的概率。
设硬币正面朝上的概率为θ,其服从伯努利分布,因此反面朝上的概率为1−θ,同时设每次投掷结果相互独立,即独立同分布,则似然为
L(θ)=θ⋅θ⋅(1−θ)⋅θ⋅(1−θ)=θ3(1−θ)2
对数似然为
LL(θ)=lnL(θ)=3lnθ+2ln(1−θ)
易证LL(θ)是关于θ的凹函数,因此对其求一阶导并令导数等于零即可求出最大值点,具体地
∂θ∂LL(θ)=∂θ∂(3lnθ+2ln(1−θ))=θ3−1−θ2=θ(1−θ)3−5θ
令上式等于0可解得θ=53,显然53也是正面出现的频率。
7.3.2 式(7.18)的解释
该式所表示的正态分布并不一定是标准正态分布,因此p(xi∣c)的取值并不一定在(0,1)之间,但是仍然不妨碍其用作"概率",因为根据朴素贝叶斯的算法原理可知,只需p(xi∣c)的值仅仅是用来比大小,因此只关心相对值而不关心绝对值。
7.3.3 贝叶斯估计
贝叶斯学派视角下的一类点估计法称为贝叶斯估计[1],常用的贝叶斯估计有最大后验估计(Maximum
A Posteriori
Estimation,简称MAP)、后验中位数估计和后验期望值估计这3种参数估计方法,下面给出这3种方法的具体定义。
设总体的概率质量函数(若总体的分布为连续型时则改为概率密度函数,此处以离散型为例)
为P(x∣θ),从该总体中抽取出的n个独立同分布的样本构成样本集
D={x1,x2,⋯,xn},则根据贝叶斯式可得,在给定样本集D的条件下,θ的条件概率为
P(θ∣D)=P(D)P(D∣θ)P(θ)=∑θP(D∣θ)P(θ)P(D∣θ)P(θ)
其中P(D∣θ)为似然函数,由于样本集D中的样本是独立同分布的,所以似然函数可以进一步展开,有
P(θ∣D)=∑θP(D∣θ)P(θ)P(D∣θ)P(θ)=∑θ∏i=1nP(xi∣θ)P(θ)∏i=1nP(xi∣θ)P(θ)
根据贝叶斯学派的观点,此条件概率代表了我们在已知样本集D后对θ产生的新的认识,它综合了我们对θ主观预设的先验概率P(θ)和样本集D带来的信息,通常称其为θ的后验概率。
贝叶斯学派认为,在得到P(θ∣D)以后,对参数θ的任何统计推断,都只能基于P(θ∣D)。至于具体如何去使用它,可以结合某种准则一起去进行,统计学家也有一定的自由度。
对于点估计来说,
求使得P(θ∣D)达到最大值的θ^MAP作为θ的估计称为最大后验估计,求P(θ∣D)的中位数θ^Median作为θ的估计称为后验中位数估计,求P(θ∣D)的期望值(均值)θ^Mean作为θ的估计称为后验期望值估计。
7.3.4 Categorical分布
Categorical分布又称为广义伯努利分布,是将伯努利分布中的随机变量可取值个数由两个泛化为多个得到的分布。具体地,设离散型随机变量X共有k种可能的取值{x1,x2,⋯,xk},且X取到每个值的概率分别为P(X=x1)=θ1,P(X=x2)=θ2,⋯,P(X=xk)=θk,则称随机变量X服从参数为θ1,θ2,⋯,θk的Categorical分布,其概率质量函数为
P(X=xi)=p(xi)=θi
7.3.5 Dirichlet分布
类似于Categorical分布是伯努利分布的泛化形式,Dirichlet分布是Beta分布的泛化形式。对于一个k维随机变量x=(x1,x2,⋯,xk)∈Rk,其中xi(i=1,2,⋯,k)满足0⩽xi⩽1,∑i=1kxi=1,若x服从参数为α=(α1,α2,⋯,αk)∈Rk的Dirichlet分布,则其概率密度函数为
p(x;α)=∏i=1kΓ(αi)Γ(∑i=1kαi)i=1∏kxiαi−1
其中Γ(z)=∫0∞xz−1e−xdx为Gamma函数,当α=(1,1,⋯,1)时,Dirichlet分布等价于均匀分布。
7.3.6 式(7.19)和式(7.20)的推导
从贝叶斯估计的角度来说,拉普拉斯修正就等价于先验概率为Dirichlet分布的后验期望值估计。为了接下来的叙述方便,我们重新定义一下相关数学符号。
设有包含m个独立同分布样本的训练集D,D中可能的类别数为k,其类别的具体取值范围为{c1,c2,...,ck}。若令随机变量C表示样本所属的类别,且C取到每个值的概率分别为P(C=c1)=θ1,P(C=c2)=θ2,...,P(C=ck)=θk,那么显然C服从参数为θ=(θ1,θ2,...,θk)∈Rk的Categorical分布,其概率质量函数为
P(C=ci)=P(ci)=θi
其中P(ci)=θi就是式(7.9)所要求解的P^(c),下面我们用贝叶斯估计中的后验期望值估计来估计θi。根据贝叶斯估计的原理可知,在进行参数估计之前,需要先主观预设一个先验概率P(θ),通常为了方便计算后验概率P(θ∣D),我们会用似然函数P(D∣θ)的共轭先验作为我们的先验概率。显然,此时的似然函数P(D∣θ)是一个基于Categorical分布的似然函数,而Categorical分布的共轭先验为Dirichlet分布,所以只需要预设先验概率P(θ)为Dirichlet分布,然后使用后验期望值估计就能估计出θi。
具体地,记D中样本类别取值为ci的样本个数为yi,则似然函数P(D∣θ)可展开为
P(D∣θ)=θ1y1...θkyk=i=1∏kθiyi
则有后验概率
P(θ∣D)=P(D)P(D∣θ)P(θ)=∑θP(D∣θ)P(θ)P(D∣θ)P(θ)=∑θ[∏i=1kθiyi⋅P(θ)]∏i=1kθiyi⋅P(θ)
假设此时先验概率P(θ)是参数为α=(α1,α2,...,αk)∈Rk的Dirichlet分布,则P(θ)可写为
P(θ;α)=∏i=1kΓ(αi)Γ(∑i=1kαi)i=1∏kθiαi−1
将其代入P(D∣θ)可得
P(θ∣D)=∑θ[∏i=1kθiyi⋅P(θ)]∏i=1kθiyi⋅P(θ)=∑θ∏i=1kθiyi⋅∏i=1kΓ(αi)Γ(∑i=1kαi)∏i=1kθiαi−1∏i=1kθiyi⋅∏i=1kΓ(αi)Γ(∑i=1kαi)∏i=1kθiαi−1=∑θ[∏i=1kθiyi⋅∏i=1kθiαi−1]⋅∏i=1kΓ(αi)Γ(∑i=1kαi)∏i=1kθiyi⋅∏i=1kΓ(αi)Γ(∑i=1kαi)∏i=1kθiαi−1=∑θ[∏i=1kθiyi⋅∏i=1kθiαi−1]∏i=1kθiyi⋅∏i=1kθiαi−1=∑θ[∏i=1kθiαi+yi−1]∏i=1kθiαi+yi−1
此时若设α+y=(α1+y1,α2+y2,...,αk+yk)∈Rk,则根据Dirichlet分布的定义可知
P(θ;α+y)θ∑P(θ;α+y)11∑θ[∏i=1kθiαi+yi−1]1=∏i=1kΓ(αi+yi)Γ(∑i=1k(αi+yi))i=1∏kθiαi+yi−1=θ∑∏i=1kΓ(αi+yi)Γ(∑i=1k(αi+yi))i=1∏kθiαi+yi−1=θ∑∏i=1kΓ(αi+yi)Γ(∑i=1k(αi+yi))i=1∏kθiαi+yi−1=∏i=1kΓ(αi+yi)Γ(∑i=1k(αi+yi))θ∑[i=1∏kθiαi+yi−1]=∏i=1kΓ(αi+yi)Γ(∑i=1k(αi+yi))
将此结论代入P(D∣θ)可得
P(θ∣D)=∑θ[∏i=1kθiαi+yi−1]∏i=1kθiαi+yi−1=∏i=1kΓ(αi+yi)Γ(∑i=1k(αi+yi))i=1∏kθiαi+yi−1=P(θ;α+y)
综上可知,对于服从Categorical分布的θ来说,假设其先验概率P(θ)是参数为α的Dirichlet分布时,得到的后验概率P(θ∣D)是参数为α+y的Dirichlet分布,通常我们称这种先验概率分布和后验概率分布形式相同的这对分布为共轭分布。在推得后验概率P(θ∣D)的具体形式以后,根据后验期望值估计可得θi的估计值为
θi=EP(θ∣D)[θi]=EP(θ;α+y)[θi]=∑j=1k(αj+yj)αi+yi=∑j=1kαj+∑j=1kyjαi+yi=∑j=1kαj+mαi+yi
显然,式(7.9)是当α=(1,1,...,1)时推得的具体结果,此时等价于我们主观预设的先验概率P(θ)服从均匀分布,此即拉普拉斯修正。同理,当我们调整α的取值后,即可推得其他数据平滑的公式。
7.4 半朴素贝叶斯分类器
7.4.1 式(7.21)的解释
在朴素贝叶斯中求解P(xi∣c)时,先挑出类别为c的样本,若是离散属性则按大数定律估计P(xi∣c),若是连续属性则求这些样本的均值和方差,接着按正态分布估计P(xi∣c)。现在估计P(xi∣c,pai),则是先挑出类别为c且属性xi所依赖的属性为pai的样本,剩下步骤与估计P(xi∣c)时相同。
7.4.2 式(7.22)的解释
该式写为如下形式可能更容易理解:
I(xi,xj∣y)=n=1∑NP(xi,xj∣cn)logP(xi∣cn)P(xj∣cn)P(xi,xj∣cn)
其中i,j=1,2,...,d且i=j,N为类别个数。该式共可得到2d(d−1)个I(xi,xj∣y),即每对(xi,xj)均有一个条件互信息I(xi,xj∣y)。
7.4.3 式(7.23)的推导
基于贝叶斯定理,式(7.8)将联合概率P(x,c)写为等价形式P(x∣c)P(c),实际上,也可将向量x拆开,把P(x,c)写为P(x1,x2,...,xd,c)形式,然后利用概率公式P(A,B)=P(A∣B)P(B)对其恒等变形
P(x,c)=P(x1,x2,…,xd,c)=P(x1,x2,…,xd∣c)P(c)=P(x1,…,xi−1,xi+1,…,xd∣c,xi)P(c,xi)
类似式(7.14)采用属性条件独立性假设,则
P(x1,...,xi−1,xi+1,...,xd∣c,xi)=j=1j=i∏dP(xj∣c,xi)
根据式(7.25)可知,当j=i时,∣Dc,xi∣=∣Dc,xi,xj∣,若不考虑平滑项,则此时P(xj∣c,xi)=1,因此在上式的连乘项中可放开j=i的约束,即
P(x1,...,xi−1,xi+1,...,xd∣c,xi)=j=1∏dP(xj∣c,xi)
综上可得:
P(c∣x)=P(x)P(x,c)=P(x)P(c,xi)P(x1,…,xi−1,xi+1,…,xd∣c,xi)∝P(c,xi)P(x1,…,xi−1,xi+1,…,xd∣c,xi)=P(c,xi)j=1∏dP(xj∣c,xi)
上式是将属性xi作为超父属性的,AODE尝试将每个属性作为超父来构建SPODE,然后将那些具有足够训练数据支撑的SPODE集成起来作为最终结果。具体来说,对于总共d个属性来说,共有d个不同的上式,集成直接求和即可,因为对于不同的类别标记c均有d个不同的上式,至于如何满足"足够训练数据支撑的SPODE"这个条件,注意式(7.24)和式(7.25)均使用到了∣Dc,xi∣和∣Dc,xi,xj∣,若集合Dxi中样本数量过少,则∣Dc,xi∣和∣Dc,xi,xj∣将会更小,因此在式(7.23)中要求集合Dxi中样本数量不少于m′。
7.4.4 式(7.24)和式(7.25)的推导
类比式(7.19)和式(7.20)的推导。
7.5 贝叶斯网
7.5.1 式(7.27)的解释
在这里补充一下同父结构和顺序结构的推导。同父结构:在给定父节点x1的条件下x3,x4独立
P(x3,x4∣x1)=P(x1)P(x1,x3,x4)=P(x1)P(x1)P(x3∣x1)P(x4∣x1)=P(x3∣x1)P(x4∣x1)
顺序结构:在给定节点x的条件下y,z独立
P(y,z∣x)=P(x)P(x,y,z)=P(x)P(z)P(x∣z)P(y∣x)=P(x)P(z,x)P(y∣x)=P(z∣x)P(y∣x)
7.6 EM算法
"西瓜书"中仅给出了EM算法的运算步骤,其原理并未展开讲解,下面补充EM算法的推导原理,以及所用到的相关数学知识。
7.6.1 Jensen不等式
若f是凸函数,则下式恒成立
f(tx1+(1−t)x2)⩽tf(x1)+(1−t)f(x2)
其中t∈[0,1],若将x推广到n个时同样成立,即
f(t1x1+t2x2+...+tnxn)⩽t1f(x1)+t2f(x2)+...+tnf(tn)
其中t1,t2,...,tn∈[0,1],∑i=1nti=1。此不等式在概率论中通常以如下形式出现
φ(E[X])⩽E[φ(X)]
其中X是随机变量,φ为凸函数,E[X]为随机变量X的期望。显然,若f和φ是凹函数,则上述不等式中的⩽换成⩾也恒成立。
7.6.2 EM算法的推导
假设现有一批独立同分布的样本{x1,x2,...,xm},它们是由某个含有隐变量的概率分布p(x,z;θ)生成,现尝试用极大似然估计法估计此概率分布的参数。为了便于讨论,此处假设z为离散型随机变量,则对数似然函数为
LL(θ)=i=1∑mlnp(xi;θ)=i=1∑mlnzi∑p(xi,zi;θ)
显然,此时LL(θ)里含有未知的隐变量z以及求和项的对数,相比于不含隐变量的对数似然函数,显然该似然函数的极大值点较难求解,而EM算法则给出了一种迭代的方法来完成对LL(θ)的极大化。
下面给出两种推导方法,一个是出自李航老师的《统计学习方法》[2],一个是出自吴恩达老师的CS229,两种推导方式虽然形式上有差异,但最终的Q函数相等,接下来先讲述两种推导方法,最后会给出Q函数是相等的证明。
首先给出《统计学习方法》中的推导方法,设X={x1,x2,...,xm},Z={z1,z2,...,zm},则对数似然函数可以改写为
LL(θ)=lnP(X∣θ)=lnZ∑P(X,Z∣θ)=ln(Z∑P(X∣Z,θ)P(Z∣θ))
EM算法采用的是通过迭代逐步近似极大化L(θ):假设第t次迭代时θ的估计值是θ(t),我们希望第t+1次迭代时的θ能使LL(θ)增大,即LL(θ)>LL(θ(t))。为此,考虑两者的差
LL(θ)−LL(θ(t))=ln(Z∑P(X∣Z,θ)P(Z∣θ))−lnP(X∣θ(t))=ln(Z∑P(Z∣X,θ(t))P(Z∣X,θ(t))P(X∣Z,θ)P(Z∣θ))−lnP(X∣θ(t))
由上述Jensen不等式可得
LL(θ)−LL(θ(t))⩾Z∑P(Z∣X,θ(t))lnP(Z∣X,θ(t))P(X∣Z,θ)P(Z∣θ)−lnP(X∣θ(t))=Z∑P(Z∣X,θ(t))lnP(Z∣X,θ(t))P(X∣Z,θ)P(Z∣θ)−1⋅lnP(X∣θ(t))=Z∑P(Z∣X,θ(t))lnP(Z∣X,θ(t))P(X∣Z,θ)P(Z∣θ)−Z∑P(Z∣X,θ(t))⋅lnP(X∣θ(t))=Z∑P(Z∣X,θ(t))(lnP(Z∣X,θ(t))P(X∣Z,θ)P(Z∣θ)−lnP(X∣θ(t)))=Z∑P(Z∣X,θ(t))lnP(Z∣X,θ(t))P(X∣θ(t))P(X∣Z,θ)P(Z∣θ)
令
B(θ,θ(t))=LL(θ(t))+Z∑P(Z∣X,θ(t))lnP(Z∣X,θ(t))P(X∣θ(t))P(X∣Z,θ)P(Z∣θ)
则
LL(θ)⩾B(θ,θ(t))
即B(θ,θ(t))是LL(θ)的下界,此时若设θ(t+1)能使得B(θ,θ(t))达到极大,即
B(θ(t+1),θ(t))⩾B(θ,θ(t))
由于LL(θ(t))=B(θ(t),θ(t)),那么可以进一步推得
LL(θ(t+1))⩾B(θ(t+1),θ(t))⩾B(θ(t),θ(t))=LL(θ(t))
LL(θ(t+1))⩾LL(θ(t))
因此,任何能使得B(θ,θ(t))增大的θ,也可以使得LL(θ)增大,于是问题就转化为了求解能使得B(θ,θ(t))达到极大的θ(t+1),即
θ(t+1)=argmaxθB(θ,θ(t))=argmaxθ(LL(θ(t))+Z∑P(Z∣X,θ(t))lnP(Z∣X,θ(t))P(X∣θ(t))P(X∣Z,θ)P(Z∣θ))
略去对θ极大化而言是常数的项
θ(t+1)=argmaxθ(Z∑P(Z∣X,θ(t))ln(P(X∣Z,θ)P(Z∣θ)))=argmaxθ(Z∑P(Z∣X,θ(t))lnP(X,Z∣θ))=argmaxθQ(θ,θ(t))
到此即完成了EM算法的一次迭代,求出的θ(t+1)作为下一次迭代的初始θ(t)。综上,EM算法的"E步"和"M步"可总结为以下两步。
E步:计算完全数据的对数似然函数lnP(X,Z∣θ)关于在给定观测数据X和当前参数θ(t)下对未观测数据Z的条件概率分布P(Z∣X,θ(t))的期望Q(θ,θ(t)):
Q(θ,θ(t))=EZ[lnP(X,Z∣θ)∣X,θ(t)]=Z∑P(Z∣X,θ(t))lnP(X,Z∣θ)
M步:求使得Q(θ,θ(t))达到极大的θ(t+1)。
接下来给出CS229中的推导方法,设zi的概率质量函数为Qi(zi),则LL(θ)可以作如下恒等变形
LL(θ)=i=1∑mlnp(xi;θ)=i=1∑mlnzi∑p(xi,zi;θ)=i=1∑mlnzi∑Qi(zi)Qi(zi)p(xi,zi;θ)
其中∑ziQi(zi)Qi(zi)p(xi,zi;θ)可以看做是对Qi(zi)p(xi,zi;θ)关于zi求期望,即
zi∑Qi(zi)Qi(zi)p(xi,zi;θ)=Ezi[Qi(zi)p(xi,zi;θ)]
由Jensen不等式可得
ln(Ezi[Qi(zi)p(xi,zi;θ)])⩾Ezi[ln(Qi(zi)p(xi,zi;θ))]
lnzi∑Qi(zi)Qi(zi)p(xi,zi;θ)⩾zi∑Qi(zi)lnQi(zi)p(xi,zi;θ)
将此式代入LL(θ)可得
LL(θ)=i=1∑mlnzi∑Qi(zi)Qi(zi)p(xi,zi;θ)⩾i=1∑mzi∑Qi(zi)lnQi(zi)p(xi,zi;θ)1◯
若令B(θ)=i=1∑mzi∑Qi(zi)lnQi(zi)p(xi,zi;θ),则此时B(θ)为LL(θ)的下界函数,那么这个下界函数所能构成的最优下界是多少?即B(θ)的最大值是多少?显然,B(θ)是LL(θ)的下界函数,反过来LL(θ)是其上界函数,所以如果能使得B(θ)=LL(θ),则此时的B(θ)就取到了最大值。根据Jensen不等式的性质可知,如果能使得Qi(zi)p(xi,zi;θ)恒等于某个常量c,大于等于号便可以取到等号。因此,只需任意选取满足Qi(zi)p(xi,zi;θ)=c的Qi(zi)就能使得B(θ)达到最大值。由于Qi(zi)是zi的概率质量函数,所以Qi(zi)同时也满足约束0⩽Qi(zi)⩽1,∑ziQi(zi)=1,结合Qi(zi)的所有约束可以推得
Qi(zi)p(xi,zi;θ)=c
p(xi,zi;θ)=c⋅Qi(zi)
zi∑p(xi,zi;θ)=c⋅zi∑Qi(zi)
zi∑p(xi,zi;θ)=c
Qi(zi)p(xi,zi;θ)=zi∑p(xi,zi;θ)
Qi(zi)=zi∑p(xi,zi;θ)p(xi,zi;θ)=p(xi;θ)p(xi,zi;θ)=p(zi∣xi;θ)
所以,当且仅当Qi(zi)=p(zi∣xi;θ)时B(θ)取到最大值,将Qi(zi)=p(zi∣xi;θ)代回LL(θ)和B(θ)可以推得
LL(θ)=i=1∑mlnzi∑Qi(zi)Qi(zi)p(xi,zi;θ)=i=1∑mlnzi∑p(zi∣xi;θ)p(zi∣xi;θ)p(xi,zi;θ)=i=1∑mzi∑p(zi∣xi;θ)lnp(zi∣xi;θ)p(xi,zi;θ)=max{B(θ)}2◯3◯4◯5◯
其中式4◯是式1◯中不等式取等号时的情形。由以上推导可知,此时对数似然函数LL(θ)等价于其下界函数的最大值max{B(θ)},所以要想极大化LL(θ)可以通过极大化max{B(θ)}来间接极大化LL(θ),因此,下面考虑如何极大化max{B(θ)}。假设已知第t次迭代的参数为θ(t),而第t+1次迭代的参数θ(t+1)可通过如下方式求得
θ(t+1)=argθmaxmax{B(θ)}=argθmaxi=1∑mzi∑p(zi∣xi;θ(t))lnp(zi∣xi;θ(t))p(xi,zi;θ)=argθmaxi=1∑mzi∑p(zi∣xi;θ(t))lnp(xi,zi;θ)6◯7◯8◯
此时将θ(t+1)代入LL(θ)可推得
LL(θ(t+1))=max{B(θ(t+1))}=i=1∑mzi∑p(zi∣xi;θ(t+1))lnp(zi∣xi;θ(t+1))p(xi,zi;θ(t+1))⩾i=1∑mzi∑p(zi∣xi;θ(t))lnp(zi∣xi;θ(t))p(xi,zi;θ(t+1))⩾i=1∑mzi∑p(zi∣xi;θ(t))lnp(zi∣xi;θ(t))p(xi,zi;θ(t))=max{B(θ(t))}=LL(θ(t))9◯10◯11◯12◯13◯14◯
其中,式9◯和式10◯分别由式5◯和式4◯推得,式11◯由式1◯推得,式12◯由式7◯推得,式13◯和式14◯由式2◯至式5◯推得。此时若令
Q(θ,θ(t))=i=1∑mzi∑p(zi∣xi;θ(t))lnp(xi,zi;θ)
由式9◯至式14◯可知,凡是能使得Q(θ,θ(t))达到极大的θ(t+1)一定能使得LL(θ(t+1))⩾LL(θ(t))。综上,EM算法的"E步"和"M步"可总结为以下两步。
E步:令Qi(zi)=p(zi∣xi;θ)并写出Q(θ,θ(t));
M步:求使得Q(θ,θ(t))到达极大的θ(t+1)。
以上便是EM算法的两种推导方法,下面证明两种推导方法中的Q函数相等。
Q(θ∣θ(t))=Z∑P(Z∣X,θ(t))lnP(X,Z∣θ)=z1,z2,...,zm∑{i=1∏mP(zi∣xi,θ(t))ln[i=1∏mP(xi,zi∣θ)]}=z1,z2,...,zm∑{i=1∏mP(zi∣xi,θ(t))[i=1∑mlnP(xi,zi∣θ)]}=z1,z2,...,zm∑{i=1∏mP(zi∣xi,θ(t))[lnP(x1,z1∣θ)+lnP(x2,z2∣θ)+...+lnP(xm,zm∣θ)]}=z1,z2,...,zm∑[i=1∏mP(zi∣xi,θ(t))⋅lnP(x1,z1∣θ)]+...+z1,z2,...,zm∑[i=1∏mP(zi∣xi,θ(t))⋅lnP(xm,zm∣θ)]
其中z1,z2,...,zm∑[i=1∏mP(zi∣xi,θ(t))⋅lnP(x1,z1∣θ)]可作如下恒等变形:
=========z1,z2,...,zm∑[i=1∏mP(zi∣xi,θ(t))⋅lnP(x1,z1∣θ)]z1,z2,...,zm∑[i=2∏mP(zi∣xi,θ(t))⋅P(z1∣x1,θ(t))⋅lnP(x1,z1∣θ)]z1∑z2,...,zm∑[i=2∏mP(zi∣xi,θ(t))⋅P(z1∣x1,θ(t))⋅lnP(x1,z1∣θ)]z1∑P(z1∣x1,θ(t))lnP(x1,z1∣θ)z2,...,zm∑[i=2∏mP(zi∣xi,θ(t))]z1∑P(z1∣x1,θ(t))lnP(x1,z1∣θ)z2,...,zm∑[i=3∏mP(zi∣xi,θ(t))⋅P(z2∣x2,θ(t))]z1∑P(z1∣x1,θ(t))lnP(x1,z1∣θ){z2∑z3,...,zm∑[i=3∏mP(zi∣xi,θ(t))⋅P(z2∣x2,θ(t))]}z1∑P(z1∣x1,θ(t))lnP(x1,z1∣θ){z2∑P(z2∣x2,θ(t))z3,...,zm∑[i=3∏mP(zi∣xi,θ(t))]}z1∑P(z1∣x1,θ(t))lnP(x1,z1∣θ){z2∑P(z2∣x2,θ(t))×z3∑P(z3∣x3,θ(t))×...×zm∑P(zm∣xm,θ(t))}z1∑P(z1∣x1,θ(t))lnP(x1,z1∣θ)×{1×1×...×1}z1∑P(z1∣x1,θ(t))lnP(x1,z1∣θ)
所以
z1,z2,...,zm∑[i=1∏mP(zi∣xi,θ(t))⋅lnP(x1,z1∣θ)]=z1∑P(z1∣x1,θ(t))lnP(x1,z1∣θ)
同理可得
z1,z2,...,zm∑[i=1∏mP(zi∣xi,θ(t))⋅lnP(x2,z2∣θ)]z1,z2,...,zm∑[i=1∏mP(zi∣xi,θ(t))⋅lnP(xm,zm∣θ)]=z2∑P(z2∣x2,θ(t))lnP(x2,z2∣θ)⋮=zm∑P(zm∣xm,θ(t))lnP(xm,zm∣θ)
将上式代入Q(θ∣θ(t))可得
Q(θ∣θ(t))=z1,z2,...,zm∑[i=1∏mP(zi∣xi,θ(t))⋅lnP(x1,z1∣θ)]+...+z1,z2,...,zm∑[i=1∏mP(zi∣xi,θ(t))⋅lnP(xm,zm∣θ)]=z1∑P(z1∣x1,θ(t))lnP(x1,z1∣θ)+...+zm∑P(zm∣xm,θ(t))lnP(xm,zm∣θ)=i=1∑mzi∑P(zi∣xi,θ(t))lnP(xi,zi∣θ)
参考文献
[1] 陈希孺. 概率论与数理统计. 中国科学技术大学出版社, 2009.
[2] 李航. 统计学习方法. 清华大学出版社, 2012.