EM算法及其实例分析

EM算法及其应用

0 引言
在统计领域里, 统计计算技术近年来发展很快, 它使许多统计方法, 尤其是Bayes 统计得到广泛的运用。Bayes计算方法有很多, 大体上可分为两大类:一类是直接应用于后验分布以得到后验均值或后验众数的估计, 以及这种估计的渐进方差或其近似;另一类算法可以总称为数据添加算法, 这是近年发展很快而且应用很广的一种算法, 它是在观测数据的基础上加一些“ 潜在数据” , 从而简化计算并完成一系列简单的极大化或模拟, 该“ 潜在数据” 可以是“ 缺损数据” 或未知参数。其原理可以表述如下:设我们能观测到的数据是Y , θ关于Y 的后验分布p(θ|Y)很复杂, 难以直接进行各种统计计算, 假如我们能假定一些没有能观测到的潜在数据Z 为已知(譬如, Y 为某变量的截尾观测值, 则Z 为该变量的真值), 则可能得到一个关于θ的简单的添加后验分布p(θ|Y , Z), 利用p(θ|Y , Z)的简单性我们可以对Z 的假定作检查和改进, 如此进行, 我们就将一个复杂的极大化或抽样问题转变为一系列简单的极大化或抽样。EM 算法就是一种常用的数据添加算法。
1　EM 算法及其理论
EM（Expectation-Maximization）算法是一种迭代方法, 最初由Dempster 等于1977 年首次提出, 主要用来计算后验分布的众数(即极大似然估计)。近十年来引起了统计学家们的极大兴趣, 在统计领域得到广泛应用。这种方法可以广泛的应用于缺损数据, 截尾数据, 成群数据, 带有讨厌参数的数据等所谓的不完全数据。
它的每一次迭代有两步组成:E 步(求期望)和M 步(极大化)。一般地, 以p(θ|Y)表示θ的基于观测数据的后验分布密度函数, 称为观测后验分布, p(θ|Y,Z)表示添加数据Z 后得到的关于θ的后验分布密度函数, 称为添加后验分布, p(Z|θ, Y)表示在给定θ和观测数据Y 下潜在数据Z 的条件分布密度函数。我们的目的是计算观测后验分布p(θ|Y)的众数, 于是, EM 算法按下述步骤进行。
记θ(i)为第i +1 次迭代开始时后验众数的估计值, 则第i +1 次迭代的两步为
E 步:将p(θ|Y,Z)或log p(θ|Y ,Z)后关于Z 的条件分布求期望, 从而把Z 积掉, 即
Q(θ|θ(i), Y) Ez[ logp(θ|Y ,Z)|θ(i),Y]
=∫log[ p(θ|Y , Z)] p(Z|θ(i), Y)dZ
(1)
M 步:将Q(θ|θ(i),Y)极大化, 即找一个点θ(i+1), 使
Q(θ(i+1)|θ(i), Y)= maxθ Q(θ|θ(i),Y)
(2)
如此形成了一次迭代θ(i) → θ(i+1)。将上述E 步和M 步进行迭代直至‖θ(i+1) -θ(i)‖ 或‖Q(θ(i+1)|θ(i), Y)-Q(θ(i)|θ(i), Y)‖ 充分小时停止。

例1　（Rao（1965））设有197种动物服从多项分布，将其分成四类，观测到的数据为y = (y1 , y2 , y3 , y4)=(125, 18 , 20, 34), 再设属于各类的概率分布为
（+ ， (1 -θ)， (1 -θ)，），
其中θ（0，1）试估计θ。
取θ的先验分布Π(θ)为(0, 1)上均匀分布, 则θ的观测后验分布为
p(θ|Y)∝Π(θ)p(Y|θ)=(+ )y1[ (1 -θ)] y2[ (1-θ)]y3()y4
∝(2+θ)y1(1 -θ)y2 +y3 θy4
(3)
假设第一种结果可以分解成两部分, 其发生概率分别为和, 令Z 和y1 -Z 分别表示试验中结果落入这两部分的次数(Z 是不能观测的潜在数据), 则θ的添加后验分布为
p(θ|Y , Z)∝Π(θ)p(Y ,Z|θ)=( )z()y1 -z [ (1-θ)] y2 [ (1-θ)] y3 ()y4 ∝ θy1 -z +y4 (1 -θ)y2 +y 3
(4)
用(3)求θ的后验众数比较麻烦, 而用(4)求后验众数则十分简单。
在第i +1 次迭代中, 假设有估计值θ(i), 则可通过E 步和M 步得到θ一个新的估计, 在E 步中有:
Q(θ|θ(i), Y)=Ez[ (y1 -Z +y4)logθ+(y2 +y3)log(1-θ)|θ(i), Y]
= [ y1 -Ez(Z|θ(i), Y)+y4] logθ+(y2 +y3)log(1 -θ)
因在θ(i) 和Y 给定下, Z ～ b(y1 ,2/(θ(i)+2)), 故Ez(Z|θ(i),Y)= 2y1/(θ(i)+2)
在M 步中, 我们将Q(θ|θ(i), Y)对θ求导并令其为0, 有
θ(i+1) =
=

(5)
(5)式给出了有EM 算法得到的迭代公式。由此公式进行迭代, 可得到观测后验分布的众数。由(5) 式不用迭代也可求出θ的估计。事实上, 它是一个迭代公式, 若收敛到 θ=(159θ(i)+68)/( 197θ(i) +144), 解
之有 θ= 0.626 821 497。

2　EM 算法的收敛性
2.1　EM 算法的收敛性定理
EM 算法的最大优点是简单和稳定。EM 算法的最大目的是提供一个简单的迭代算法来计算后验众数（或MLE极大似然估计）, 我们不禁问,由EM 算法得到的估计序列是否收敛？如果收敛, 其结果是否为p(θ|Y)的最大值或局部最大值。为此我们给出以下两个定理。记EM 算法得到的估计序列为θ(i) , i = 1, 2, …，L(θ|Y)= logp(θ|Y)。
定理1　EM 算法在每一次迭代后均提高(观测)后验密度函数值, 即
p(θ(i+1)|Y)≥ p(θ(i)|Y)
(6)
证明　由全概率公式（乘数公式）有
p(θ, Z |Y)= p(Z |θ, Y)p(θ|Y)= p(θ|Y , Z)p(Z |Y)
将上式后面两项取对数有
logp(θ|Y)= logp(θ|Y , Z)-logp(Z |θ, Y)+logp(Z |Y)
(7)
设现有估计θ(i), 将上式对Z 关于p(Z |θ(i), Y)求期望, 有
log p(θ|Y)
=∫[ logp(θ|Y , Z)-logp(Z |θ, Y)+logp(Z |Y)] p(Z |θ(i), Y)dZ
==Q(θ|θ(i), Y)-H(θ|θ(i), Y)+K(θ(i) , Y)
(8)
其中Q(θ|θ(i), Y)已在(1)中定义,
H(θ|θ(i), Y)=∫log[p(Z|θ,Y)] p(Z |θ(i), Y)dZ ,
K(θ(i), Y)=∫log[ p(Z |Y)] p(Z |θ(i), Y)dZ
分别在(8)中取θ为θ(i) 和θ(i+1)并相减, 有
logp(θ(i+1)|Y)-logp(θ(i)|Y)= [ Q(θ(i+1)|θ(i), Y)-Q(θ(i)|θ(i), Y)] -[ H(θ(i+1)|θ(i), Y)-H(θ(i)|θ(i), Y)]
由Jensen 不等式, 知
Ez[log()|θ(i), Y] log{ Ez()|θ(i), Y]}= 0
故
H(θ(i+1)|θ(i), Y)-H(θ(i)|θ(i),Y)≤ 0,又θ(i+1) 是使Q(θ|θ(i), Y)达到最大的, 故
Q(θ(i+1)|θ(i), Y)-Q(θ(i)|θ(i), Y)≥ 0 .
所以（6）成立。这说明用EM 算法求出的θ(i)的确使其对数似然L(θ(i)|Y)关于i 单增, 但只有在严格单增的情况下, 才能保证EM 算法求出的θ(×) 为θ的极大似然估计。这是从纯数学的, 实际应用有一定的困难。
定理2　(1)如果p(θ|Y)有上界, 则L(θ(i)|Y)收敛到某个L*。
(2)如果Q(θ|φ)关于θ和φ都连续, 则在关于L 的很一般的条件下, 由EM 得到的估计序列θ(i) 的收敛值θ* 是L 的稳定点。
证明　由定理1的证明及单调收敛定理易知（1）成立;(2)的证明见文献[ 4]。
定理2表明EM 算法的结果只能保证收敛到后验密度函数的稳定点, 并不能保证收敛到极大值点, 事实上, 任何一种算法都很难保证其结果为极大值点。通常选取几个不同的初值进行迭代, 然后在诸估计间加以选择,以减轻初值选取对结果的影响。
2.2　EM 算法收敛速度的探讨
EM 算法是一种求参数极大似然估计的迭代算法, 在处理不完全数据中有重要应用。EM 算法实现简单, 数值计算稳定, 存储量小, 并具有良好的全局收敛性。但是, EM 算法收敛速度相当慢, 只是次线性的收敛速度, 这个缺点防碍了EM 算法的应用。现已提出了多种加速EM 算法收敛的方法, 其中使用非线性规划中的Broyden 对称秩1 校正公式, 提出了一种加速EM 算法收敛的方法。与其他加速收敛方法相比, 本方法简便易行, 不必对不完全数据的似然函数一维搜索, 在收敛性方面, 与EM 算法一样具有全局收敛性。数值试验结果表明, 本方法的收敛速度比EM 算法快的多。
一般地, 在EM 算法的M 步求θ(i+1)时, 可求解方程组 =0 , 得到参数的迭代公式θ(i+1) = G(θ(i))。这种迭代公式通常使EM 算法的收敛速度很慢, 为加速收敛, 可考虑使用其他方法求解。
根据著名的Fisher 公式 θ=θ(i) = g(θ(i)), 这里g(θ(i))是不完全数据对数似然函数L(θ)在θ(i)处的梯度g(θ(i))=▽ L(θ)|θ=θ(i) , 于是问题转化为求θ(), 使g(θ())=0 , 从而可以使用非线性规划中的有效方法求解, 达到加速收敛的目的。
在非线性规划的求解方法中,Broyden 对称秩1 校正公式具有二次终止性和超线性的收敛速度, 使用它可
望改善EM 步长的缓慢变化, 以加速EM 算法的收敛。但是, 如果不对L(θ)进行一维搜索,Broyden 对称秩1 校正公式仅具有局部收敛性, 因此考虑将Broyden 对称秩1 校正公式的二次终止性和EM 算法良好的全局收敛
性相结合, 形成一种能够加速EM 算法收敛的新方法——— EMB 算法。在迭代初期, 使用EM 算法, 而当接近最优解时, EM 步长的变化极为缓慢, 这时使用Broyden 对称秩1 校正公式进行校正。
3　EM 算法的应用实例
EM 算法可以应用于医学研究中, 尤其是临床医学中十分常见的一种数据观测形式为重复观测, 其特点是在同一实验单位上进行多次重复观测, 这个过程由于各种原因经常导致实验观测数据缺失, 如动物的意外死亡, 记录仪器发生故障, 被调查者拒绝回答相关调查项目等, 然而, 目前绝大多数重复观测数据统计分析方法都有一个基本假定, 即每个实验单位具有完全的重复观测数据(通常假定服从多元正态分布)。因此有必要对缺失数据进行适当的处理, 否则数据分析时, 通常的做法是删去具有缺失的部分观察记录而不考虑记录数据所蕴涵的信息, 从而造成信息的损失以及分析结果的偏性。
缺失数据的处理方式在很大程度上依赖医学数据观测的实际背景, 例如, 常规的实验设计得到的独立结构数据方差分析前的缺项估计的原理, 都是解误差离均差平方和对缺项变量的导数等于0 的方程, 从而获得缺项的估计值。而医学重复观测设计得到的数据有别于独立结构数据, 其数据间存在一定的相
关性, 因此需要进行专门的缺项估计方法的研究。为此。我们在软件编程的基础上, 应用EM 算法进行了这方面的数据处理。
目前, 可用于进行缺项估计的方法很多, 但EM 算法较为优秀, 它是利用所有的资料信息来进行缺项估计, 以“ 补缺” 的方式将含有缺失值的“ 不完全资料” 转化为“ 完全资料” , 因此, 对EM 算法得到的 “ 完全资料” 可以应用所有的重复观测资料分析方法进行处理分析。EM 方法补缺使估计系数的方差明显变小, 从而提高了生长曲线系数的估计精度。
当医院科研数据量较大时, 为方便缺失值的估计, 就需用MA TLAB 进行编程, 另外, 在实际工作中为了了解多个指标间的关系及变化规律, 常常需要同时观测个体的多个反应指标, 此时, 虽然数据不是重复观测数据, 但如果缺失数据随机发生, 即缺失数据发生的可能性不被该数据值的大小所影响, 则同样可用EM 算法进行估计, 因此EM 算法不但适用于重复观测缺失数据的处理, 而且适用于一般情况下多变量缺失数据的处理分析, 其适用范围更广泛。

参考文献:
[1] 肖枝洪，朱强 统计模拟及R实现[M] .武汉：武汉大学出版社，2010.4
[2] 茆诗松, 王静龙,濮小龙.高等数理统计[ M] .第2版.北京:高等教育出版社, 2006.
[3] Wu C F J .On th e convergence properties of the EMalgorithm[ J] .Th e Annals of St ati sti cs , 1983.
[4] 高旅瑞, 陈志, 王家润.一种加速EM 算法收敛的方法[ J] .数理统计与应用概率, 1998.
[5] 陈长生, 王彤, 徐勇勇, 尚磊.医学科研中缺失数据的EM 估计[ J] .第四军医大学学报, 2002 .
[6]杨基栋 EM算法理论及其应用 上海：华东师范大学出版社[J] 2009.
你好！ 这是你第一次使用 Markdown编辑器 所展示的欢迎页。如果你想学习如何使用Markdown编辑器, 可以仔细阅读这篇文章，了解一下Markdown的基本语法知识。