《EM算法及应用实例.docx》由会员分享，可在线阅读，更多相关《EM算法及应用实例.docx（13页珍藏版）》请在三一文库上搜索。
1、学院：专业：任课 老 师：姓名：统计学模拟及其 R实现论文报告论文 （设计） 题目：EM 算法及应用实例数学与统计学院2014 级统计学 黄恒振 钟思敏EM 算法及应用实例【内容摘要】 EM 算法是一种迭代算法 , 主要用来计算后验分布的众数或 极大似然估计 , 广泛地应用于缺损数据、 截尾数据、成群数据、 带有讨厌参数的 数据等所谓的不完全数据的统计推断问题。在介绍EM 算法的基础上 , 针对 EM算法收敛性定理 ,给出了 EM 的实值实例 , 最后给出 EM算法的医学应用。【关键词】 Bayes；后验分布；极大似然估计；医学实例【 Abstract 】 The EM algorithm i
2、s an iterative algorithm, is used to calculate the posterior distribution of the number or maximumli kelihood estimation, widely used in the defect data, of according to, groups of data, with hate parameters such as the so-called incomplete data statistical inference problem. Based on the introducti
3、on of the EM algorithm, in view of the EM algorithm convergence theorem, gives a real value of EM instance, the medical application of EM algorithm is given.【 Key words 】 Bayes theorem ；Posterior distribution ； Maximuml ikelihood estimation ； Medical instance目录【内容摘要】 2【关键词】 20 引言 41 EM 算法及其理论 42 EM
4、算法的收敛性 62.1 EM 算法的收敛性定理 62.2 EM 算法收敛速度的探讨 83 EM 算法的应用实例 9参考文献 : 110 引言在统计领域里 , 统计计算技术近年来发展很快 , 它使许多统计方法 , 尤其 是 Bayes 统计得到广泛的运用。 Bayes 计算方法有很多 , 大体上可分为两大类 : 一类是直接应用于后验分布以得到后验均值或后验众数的估计 , 以及这种估计 的渐进方差或其近似 ; 另一类算法可以总称为数据添加算法 , 这是近年发展很快 而且应用很广的一种算法 , 它是在观测数据的基础上加一些 “ 潜在数据” , 从 而简化计算并完成一系列简单的极大化或模拟 , 该“
5、潜在数据” 可以是“ 缺 损数据” 或未知参数。其原理可以表述如下 : 设我们能观测到的数据是 Y , 关于 Y 的后验分布 p( |Y) 很复杂, 难以直接进行各种统计计算 , 假如我们能 假定一些没有能观测到的潜在数据 Z 为已知 (譬如, Y 为某变量的截尾观测值 , 则Z 为该变量的真值 ), 则可能得到一个关于 的简单的添加后验分布 p(|Y , Z), 利用 p(|Y , Z) 的简单性我们可以对 Z 的假定作检查和改进 , 如此进行 , 我们就将一个复杂的极大化或抽样问题转变为一系列简单的极大化或抽样。 EM 算法就是一种常用的数据添加算法。1 EM 算法及其理论EM ( Exp
6、ectation-Maximization )算法是一种迭代方法 , 最初由 Dempster 等于 1977 年首次提出 , 主要用来计算后验分布的众数 ( 即极大似然估计 ) 。近十 年来引起了统计学家们的极大兴趣 , 在统计领域得到广泛应用。 这种方法可以广 泛的应用于缺损数据 , 截尾数据 , 成群数据 , 带有讨厌参数的数据等所谓的不 完全数据。它的每一次迭代有两步组成 :E 步(求期望)和 M 步(极大化) 。一般地 , 以 p( |Y) 表示的基于观测数据的后验分布密度函数 , 称为观测后验分布 , p( |Y,Z) 表示添加数据 Z 后得到的关于 的后验分布密度函数 , 称为添
7、加后验分 布, 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
8、) , 使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 = (y 1 , y 2 , y 3 , y 4)=(125, 18 , 20, 34),再设属于各类的概率分布为( + ， (1 - ) ， (1 - ) ， )，其中 (0，1)试估计 。取的先验分布()为(0, 1) 上均匀分布
9、 , 则的观测后验分布为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
10、次迭代中, 假设有估计值 (i), 则可通过 E 步和M 步得到一个新的 估计, 在 E 步中有 :Q(| (i) , Y)=Ez (y 1 -Z +y 4)log +(y2 +y 3)log(1- )| (i), Y= y 1 -Ez(Z| (i) , Y)+y 4 log +(y 2 +y 3)log(1 - ) 因在(i) 和Y 给定下, Z b(y 1,2/( (i) +2), 故Ez(Z| (i) ,Y)= 2y1/( (i) +2) 在 M 步中, 我们将 Q(|(i) , Y) 对求导并令其为 0, 有(i+1) (5)(5) 式给出了有 EM 算法得到的迭代公式。 由此公式进行
11、迭代 , 可得到观测后 验分布的众数。 由(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 =
12、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)+lo
13、gp(Z |Y) p(Z | (i), Y)dZ =Q(| (i) , Y)-H( |(i) , Y)+K( (i) , Y)(8) 其中Q(|(i), Y) 已在(1) 中定义,H(|(i), Y)= logp(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) i 单增, 但只有在严格单增的情况下 , 才能保证 EM 算法求出的 (
14、 ×) 为的 极大似然估计。这是从纯数学的 , 实际应用有一定的困难。 定理 2 (1) 如果 p(|Y)有上界, 则L( (i) |Y)收敛到某个 L*。- H( (i+1)|(i), Y)-H( (i)|(i), Y)由 Jensen 不等式 , 知 Ezlog()|(i) , Y log Ez()|(i) , Y= 0H( | 故(i+1) (i) , Y)-H(i) |(i) ,Y) 0, 又 (i+1) 是使 Q(|(i) , Y) 达到最大的 ,Q(i+1)所以( 6)成立。这说明用(i) , Y)-Q( (i) |(i) , Y) EM 算法求出的(i) 的确使其对数似
15、然 L( (i) |Y) 关于(2) 如果 Q(| )关于和都连续, 则在关于 L 的很一般的条件下 , 由EM 得 到的估计序列 (i) 的收敛值* 是 L 的稳定点。证明 由定理 1 的证明及单调收敛定理易知( 1)成立;(2) 的证明见文献 4 。 定理 2 表明 EM 算法的结果只能保证收敛到后验密度函数的稳定点 , 并不能 保证收敛到极大值点 , 事实上 , 任何一种算法都很难保证其结果为极大值点。 通 常选取几个不同的初值进行迭代 , 然后在诸估计间加以选择 , 以减轻初值选取对 结果的影响。2.2 EM 算法收敛速度的探讨EM 算法是一种求参数极大似然估计的迭代算法 , 在处理不
16、完全数据中有重 要应用。 EM 算法实现简单 , 数值计算稳定 , 存储量小 , 并具有良好的全局收敛 性。但是 , EM 算法收敛速度相当慢 , 只是次线性的收敛速度 , 这个缺点防碍了 EM 算法的应用。现已提出了多种加速 EM 算法收敛的方法 , 其中使用非线性规 划中的 Broyden 对称秩 1 校正公式 , 提出了一种加速 EM 算法收敛的方法。与 其他加速收敛方法相比 , 本方法简便易行 , 不必对不完全数据的似然函数一维 搜索, 在收敛性方面 , 与 EM 算法一样具有全局收敛性。 数值试验结果表明 , 本 方法的收敛速度比 EM 算法快的多。一般地, 在 EM 算法的 M 步
17、求(i+1) 时, 可求解方程组 =0 , 得到 参数的迭代公式 (i+1) = G(i) ) 。这种迭代公式通常使 EM 算法的收敛速度很慢 , 为加速收敛 , 可考虑使用其他方法求解。根据著名的 Fisher 公式 =(i) = g( (i) ), 这里 g( (i) ) 是不完全数据对数似然函数 L()在(i)处的梯度 g(i) )= L( )| =(i) , 于是问题 转化为求(*) , 使 g( (*) )=0 , 从而可以使用非线性规划中的有效方法求解 , 达 到加速收敛的目的。在非线性规划的求解方法中 ,Broyden 对称秩 1 校正公式具有二次终止性 和超线性的收敛速度 ,
18、使用它可 望改善 EM 步长的缓慢变化 , 以加速 EM 算法的收敛。但是 , 如果不对 L( ) 进 行一维搜索 ,Broyden 对称秩 1 校正公式仅具有局部收敛性 , 因此考虑将 Broyden 对称秩 1 校正公式的二次终止性和 EM 算法良好的全局收敛 性相结合 , 形成一种能够加速 EM 算法收敛的新方法 EMB 算法。在迭代 初期, 使用 EM 算法, 而当接近最优解时 , EM 步长的变化极为缓慢 , 这时使用 Broyden 对称秩 1 校正公式进行校正。3 EM 算法的应用实例EM 算法可以应用于医学研究中 , 尤其是临床医学中十分常见的一种数据观 测形式为重复观测 ,
19、其特点是在同一实验单位上进行多次重复观测 , 这个过程 由于各种原因经常导致实验观测数据缺失 , 如动物的意外死亡 , 记录仪器发生 故障 , 被调查者拒绝回答相关调查项目等 , 然而 , 目前绝大多数重复观测数据 统计分析方法都有一个基本假定 , 即每个实验单位具有完全的重复观测数据 ( 通 常假定服从多元正态分布 ) 。因此有必要对缺失数据进行适当的处理 , 否则数据 分析时 , 通常的做法是删去具有缺失的部分观察记录而不考虑记录数据所蕴涵 的信息 , 从而造成信息的损失以及分析结果的偏性。缺失数据的处理方式在很大程度上依赖医学数据观测的实际背景 , 例如 , 常规的实验设计得到的独立结构
20、数据方差分析前的缺项估计的原理 , 都是解误 差离均差平方和对缺项变量的导数等于 0 的方程 , 从而获得缺项的估计值。而 医学重复观测设计得到的数据有别于独立结构数据 , 其数据间存在一定的相 关性 , 因此需要进行专门的缺项估计方法的研究。 为此。我们在软件编程的基础 上, 应用 EM 算法进行了这方面的数据处理。目前 , 可用于进行缺项估计的方法很多 , 但 EM 算法较为优秀 , 它是利用 所有的资料信息来进行缺项估计 , 以“ 补缺” 的方式将含有缺失值的“ 不完 全资料” 转化为“ 完全资料” , 因此, 对 EM 算法得到的 “ 完全资料” 可 以应用所有的重复观测资料分析方法进
21、行处理分析。 EM 方法补缺使估计系数的 方差明显变小 , 从而提高了生长曲线系数的估计精度。当医院科研数据量较大时 , 为方便缺失值的估计 , 就需用 MA TLAB 进行编 程, 另外 , 在实际工作中为了了解多个指标间的关系及变化规律 , 常常需要同 时观测个体的多个反应指标 , 此时 , 虽然数据不是重复观测数据 , 但如果缺失 数据随机发生 , 即缺失数据发生的可能性不被该数据值的大小所影响 , 则同样 可用 EM 算法进行估计 , 因此 EM 算法不但适用于重复观测缺失数据的处理 , 而 且适用于一般情况下多变量缺失数据的处理分析 , 其适用范围更广泛。参考文献 :1 肖枝洪，朱强
22、 统计模拟及 R实现M . 武汉：武汉大学出版社， 2010.42 茆诗松, 王静龙,濮小龙.高等数理统计 M .第2版.北京:高等教育出版社 2006.3 Wu C F J .On th e convergence properties of the EMalgorithm J .The Annals of St ati sti cs , 1983.4 高旅瑞, 陈志, 王家润.一种加速 EM 算法收敛的方法 J . 数理统计与应 用概率, 1998.5 陈长生, 王彤, 徐勇勇, 尚磊. 医学科研中缺失数据的 EM 估计 J . 第四 军医大学学报 , 2002 .6 杨基栋 EM算法理论及其应用 上海：华东师范大学出版社 J 2009.