（概率论与数理统计专业论文）基于分组数据回归系数的估计.pdf

复旦大学硕士学位论文 摘要 随机变量的分组观察值是指在随机试验中，我们只知道随机变量x 是否落入某一已 知区问 t j 一1 ，t j ) ，而不知道随机变量x 的具体观察值在医学研究和经济指标的分析中 常常会遇到这类数据，因此研究这类数据具有重要的理论意义和实际意义近年来，对 于这类数据的分析研究得到了越来越多的关注，一部分问题己经得到解决，但这一领域 还有很多问题有待研究本文研究基于分组数据线性回归模型中回归系数的估计问题。 全文共分为五章第一章首先介绍基于分组数据的线性回归模型提出的背景以及利 用经典的统计方法m l e 对回归系数进行估计时所遇到的困难；第二章将介绍一种处理不 完全数据时的常用算法“e m 算法”以及其扩展算法“s e m 算法”和“e c m 算法”；第三章 讨论了分组数据下经典线性回归模型参数m l e 的存在唯一性，并利用e m 算法s e m 算 法分别获得其近似解和渐近协方差，通过模拟表明此方法的可行性与有效性第四章 利用e c m 算法、s e m 算法分别获得了分组数据下分组异方差线性回归模型参数m l e 的 近似解和渐近协方差，通过模拟表明此方法的可行性与有效性。第五章利用e c m 算法 和s e m 算法分别获得了分组数据下多元线性回归模型参数m l e 的近似解和渐近协方差， 通过模拟表明此方法的可行性与有效性 关键词：m l e ，分组数据，e m 算法，s e m 算法，e c m 算法 复旦大学硕士学位论文 i i a b s t r a c t ag r o u p e do b s e r v a t i o nm e a n st h a txi sk n o w ne i t h e rl i e si n s i d ea ni n t e r v a l 阢一1 ，乃) ， o rn o t h o w e v e r ，w ea r en o ta b l et oo l x * r v et h ee x a c tv a l u eo fx i ti sv e r yi m p o r t a n t b o t hi np r a c t i c ea n di nt h e o r yt os t u d yt h i sk i n do fd a t a ，f o ri to f t e no c c u r si nm e d i c a l r e s e m x c ha n de c o n o m i ca n a l y s i s r e c e n t l yt h i sk i n do fq u e s t i o ni su n d e rd i s c u s s i o n ，a n d s o m eo ft h e mh a v eb e e ns o l v e d i ns p i t eo ft h a t ，t h e r ea r em a n yu n r e s o l v e dq u e s t i o n s i nt h i sf i e l d i nt h i sp a p e r ，w et r yt os t u d yt h ep r o b l e mo fe s t i m a t i n gt h er e g r e s s i o n c o e f f i c i e n t si nt h el i n e a rr e g r e s s i o nm o d e l sb a s e do ng r o u p e dd a t a t h ep a p e rc o n s i s t so f5c h a p t e r s i nc h a p t e r1 w ei n t r o d u c eb a c k g r o u n da tw h i c h t h el i n e a rr e g r e s s i o nm o d e l sb a s e do ng r o u p e dd a t aa r ep r o p o s e da n dt h ep r o b l e me l l - c o u n t e r e dd u r i n gt h ep r o c e s so fc a r r y i n go u tt h ec l a s s i c a l 岫l l l nl i k e l i h o o dm e t h o d i nc h a p t e r2 ，w ei n t r o d u c ea nu s u a la l g o r i t h mf o rh a n d l i n gi n c o m p l e t ed a t a ，e ma l g o - r i t h m ，a n di t se x t e n d e da l g o r i t h m s ，s e ma l g o r i t h ma n de c ma l g o r i t h m i nc h a p t e r3 ， w ed i s c u s st h ee x i s t b e n c ea n du n i q u e n e s so fm l ef o rc l a s s i c a ll i n e a rr e g r e s s i o nm o d e l s b a s e d0 1 1g r o u p e dd a t a a n du s ee ma l g o r i t h ma n ds e ma l g o r i t h mt oo b t a i nt h ea p - p r o x i m a t es o l u t i o no ft h em l e a n di t sa s y m p t o t i cc o v a r i a n c e - v a r i a n c er e s p e c t i v e l y i n c h a p t e r4 ，w eu s ee c ma l g o r i t h mt oo b t a i nt h ea p p r o x i m a t es o l u t i o no ft h em l ef o r g r o u p e dh e t e r o s c e d a s t i ea n dg r o u p e dd a t al i n e a rr e g r e s s i o nm o d e la n ds e ma l g o r i t h m t oe s t i m a t ei t sa s y m p t o t i cc o v a r i a n c e - v a r i a n c e i nc h a p t e r5 w el l 舱e c ma l g o r i t h mt o o b t a i nt h ea p p r o x i m a t es o l u t i o no ft h em l ef o rm u l t i p l el i n e a rr e g r e s s i o nm o d e lw i t h g r o u p e dd a t aa n ds e ma l g o r i t h mt oe s t i m a t ei t sa s y m p t o t i cc o v a r i a n c e - v a r i a n c e b y s i m u l a t i o n sw es h o wt h a tt h e s em e t h o d sa r ea v a i l a b l ea n de f f e c t i v e k e yw o r d s ：m l e ，g r o u p e dd a t a ，e ma l g o r i t h m ，s e ma l g o r i t h m ，e c ma l g o r i t h m 论文独创性声明 本论文足我个人在导师指导下进行的研究【作及取得的研究成果。论文中除 了特别加以标注和致谢的地方外，不包含其他人或其它机构已经发表或撰写过的 研究成果。其他同志对本研究的启发和所傲的贡献均已在论文中作了明确的声明 并表示了谢意。 作者签名：衄日期：丝1 2 ：主：型 论文使用授权声明 本人完全了解复旦大学有关保留、使用学位论文的规定，即：学校有权保留 送交论文的复印件，允许论文被查阅和借阅；学校可以公布论文的全部或部分内 容，可以采用影印、缩印或其它复制手段保存论文。保密的论文在解密后遵守此 规定。 作者签名：虹导师签名： 獬 复旦大学硕士学位论文1 第一章引言 在现实世界中，我们经常要研究变量y 与y 的相关关系。由于这种关系在很多时候 可以被线性回归模型很好的描述，所以线性回归模型成为现代统计学中应用最为广泛的 模型之一 对n 组独立观测k ，玑) ，经典线性回归模型假定因变量鼽与自变量墨具有如下关系： 弘= 卢- i - q ，其中e f 独立同分布于( o ，矿) ， = 1 ，n ( 1 1 ) 为估计模型( 1 1 ) 中的参数( 卢，矿) ，通常采用最小二乘法或极大似然法尽管这两种方法 的思想有所不同，但它们所获得的估计量却是一致的。要获得此种形式的估计量，除了 要满足线性性、独立性，正态性这三个条件外。还要求所观测到的y l 是一个确切值。 然而，在实践中出于对个人隐私、节省成本等因素的考虑，经常无法获得执的确切 值，而只是知道玑落于事先确定的某一区问i 乃一l ，乃) 。我们称以此种形式出现的数据为 分组数据。它显然不同于传统的数据结构。因此，即使在经典模型的线性性、独立性、 正态性三条件都继续满足的条件下，也不能直接套用已有结果来实施统计分析。那么， 对于分组数据下线性回归模型的参数估计问题该如何解决昵。 本文将主要研究分组数据下经典线性回归模型、分组异方差线性回归模型和多元 线性回归模型的参数估计问题。由于此问题仍属于参数型的统计问题，我们自然地想 到利用极大似然估计( m l e ) 去估计模型中的各参数。但是，数据结构的改变使得我们 无法获得m l e 的显式解，此时通常是通过采用一种适当的统计算法去获得m l e 的近似 解。 近年来，随着对不完全数据分析研究的加深，出现了许多新颖的统计算法。e m 算 法就是这些算法中发展很快且应用很广的一种，它不直接对复杂的后验分布进行极大 化或进行模拟，而是在观测数据的基础上人为添加一些数据，从而简化计算并完成一系 列简单的极大化和模拟另外为了弥补e m 算法的某些缺陷，又出现某些e m 算法的变 例，如e c m 算法、m c e m 算法等。这些算法为我们进行分组数据的参数估计提供了新 的方法和思路。本文则要利用这些算法来获取分组数据下以上三种回归模型m l e 的近 似解。 为了做区间估计与假设检验，我们还有必要算出m l e 的大样本渐近协方差阵。为 l l p 。j v i e n g 和r u b i n 提出了e m 的补充算法s e m 算法。在s e m 方法下，只需要完全数据下渐 近协方差阵的程序代码、e m 算法的程序代码以及矩阵规化的标准程序代码就可获得估 计的渐近协方差阵。相较于b o o t s t r a p 和j a c k k n i f e 算法几千次迭代后的重抽样，s e m 算 复旦大学硕士学位论文 2 法显然更具有吸引力，故本文采用此种算法来算出分组数据下回归模型参数的m l e 的 渐近协方差阵。 复旦大学硕士学位论文 3 第二章e m 算法及其变例 2 1e m 算法的e 步和m 步 e m 算法【3 ，6 ，4 ，1 1 是一种迭代方法，最初由d e m p s t e r 、l a i r d 和r u b i n 在1 9 7 7 年提出， 主要用来求后验分布的众数( 即最大似然估计) ，它的每一次迭代由两步组成：e 步( 期 望步) 和m 步( 极大化步) 。记口为我们所关心的未知参数。一般地，以f ( o l 赢) 表示p 的 基于观测数据的后验分布密度函数，称为观测后验分布；f ( o l l k ，。) 表示添加数 据k 胁后得到的关于0 的后验分布密度函数，称为添加后验分布；，( ，磊妇限赢) 表示在 给定日和观测数据i ，如下潜在数据。的条件分布密度函数。我们的目的是计算观测后 验分布，( 口i y k ) 的众数。于是e m 算法如下进行。记鲫为第z + 1 次迭代开始时后验众数 的估计值，则第z + 1 次迭代的两步为： e 步：将，( 口i ，k ，。) 或l o g l ( o i i 幺，v 胁) 关于。的条件分布求期望，从而把k 协积掉， 即 q ( 0 1 0 ( o , ) = e 1 0 9 f ( 口i y 咖，) 即，y 咖】 = 1 0 9 f ( pj ，k ，) 】，( ，施1 ，l k ) d y 赢 ( 2 1 ) j m 步：将q ( 酬口( o ，l k ) 极大化，即找一个点0 0 ) ，使 q ( p + 1 1 0 ( o , l 乞。) = m qp p “， ，咖) ( 2 2 ) 由此形成了一次f 1 3 0 ( 0 一卵+ 1 ) 的迭代，按上述e 步和m 步进行迭代，直到| i 卯+ 1 ) 一卵 1 1 或 者l i q ( 口( 1 + 1 ) 眇) ，l k ) 一q ( o ( o l o ( “， ，咖) i i 充分小时停止。 2 2e m 算法的收敛性 刚刚所描述的e m 算法内在地定义了一个由0 的参数空间e 到它自己的映射：0 一 m ( 口) ， p ( + 1 ) = m ( 卵) ，z = 0 ，l ， 如果口( 】收敛于某个点矿且m ( 口) 是连续的，则0 + 必然满足 矿= m ( 0 + 、 复旦大学硕士学位论文 4 因此，在矿的邻域内进行泰勒展开，我们有 p ( + 1 ) 一矿( 0 ( o 一矿) m 7 p ) ，( 2 3 ) 其中， r ( 口) = ( o u 溉i ( o ) ) 是m ( 口) = ( m ( 口) ， 毛( 口) ) 的d d 维j a b i 缸矩阵。因此， 当d = 1 时，e m 算法是线性收敛的。事实上，对于d 1 ，如果观测信息阵一l ，( p i ，k ) 是 正定的，则e m 算法仍是线性收敛的。 e m 局部收敛速率如下定义： 口：h m 嵝：! 二哑 “o 。i i o ( o o l i 可以证明，当观测信息阵一l ，( 6 i l ，幽) 正定时，p 等于 ( a ) 的最大特征根。在第2 3 一节中， 我们将说明m ( 百) 反映了缺失信息的比例。因此，p 有效地度量了缺失信息的总体比例。 概念上来说，缺失信息的比例等于1 减去完全信息中观测信息的比例。于是，当缺失信 息比例较大时，e m 算法的收敛速度会变得比较慢。e m 算法的线性收敛相对于牛顿的 二阶收敛是极慢的，但是它的简单易行与对似然函数的稳定提升使得它在很多情况下 是非常有吸引力的。d e m p s t e r 、l a i r d 和r u b i n 还进一步地说明了在下述情况下e m 算法 是可靠地收敛的。 1 如果f ( 口i i k ) 有上界，则j ( 删l ，赢) 收敛到某个p 。 2 如果i ( y i o ) 是广义指数族，b - i ( o i y , ) 有界，那么z ( 口( i ，k ) 收敛到一个稳定值f 。 3 如果f ( y l o ) 是一个正则指数族，且z ( 口i ，k ) 有界，那么删收敛到一个稳定点0 2 3s e m 算法 e m 算法除了在收敛速度方面受到批评外，另一个被学者们经常批评的地方就是它 不像牛顿算法一样在获得m l e 近似解的同时获得了观测信息阵的估计。此估计的获得 等同于获得了m l e 的渐近协方差，它为区间估计和假设检验打下了基础。为了改进这 个缺陷，m e n g 和r u b i n 提出了补充的e m 算法，即s e m 算法【7 ，6 ，4 】来实现对渐近协差阵 的估计。下面我们简要介绍一下s e m 算法的想法和实施步骤。完全数据y 的分布可因子 化为 f ( y l o ) = ，( y 咖i o ) f ( 一y 咖，口) ， 对应的似然分解为 t ( o l ) = l ( o w ) 一l o g ( - i y 幽，口) ， ( 2 4 ) 复旦大学硕士学位论文5 其中，z ( 驯l k ) 是要极大化的已观测到的似然，l ( o i y ) 是完全数据似然，l o g ，( k 协l l k ，口) 是 完全数据似然的缺失部分。此外，我们还假定l ( o i y ) 是很容易极大化的。 对( 2 4 ) 式微分两次，对任意i k 有 j ( 卅，k ) = i ( 口l y ) 一a 21 0 9 ，( l l k ，o ) a o a p ， 其中，i ( 0 1 y ) 是基于y 已观测到的信息，最后一项的负值是来自 赢。的缺失信息。 对给定l k 和口的“的分布取期望值，得 = 一，( 2 5 ) 其中， 厶妇= e 【一o i l 。g s 瓦( y f m 刁i , 歹l r o h , 一o ) i v - 幽，o l i o s ， = s ( o l ) ， j c 。= b i s ( o l y ) i y o h , ，o l o ：p 这有一个明显的解释： 已观测到的信息= 完全信息- 缺失信息， 重排式( 2 5 ) ，且在两侧乘以仨，得 d m = 焉= i 一焉，( 2 6 ) 这几矩阵d m 代表缺失信息的比例，它是e m 映射的斜率，控制着e m 收敛的速度。缺失 信息比例越大，则收敛速率越低。 ( 2 6 ) 式蕴含着瓦：= 坛k ( j r d m ) 一，即 = 。( ，一d m ) ， 其中，v k = 瑶，。= 焉是已观测到的数据及完全数据的协方差矩阵。因此， k b = ( ，一d m + d m ) ( i d m ) 一1 = + k 上式中 a v = d m ( i d m ) 一1 ， 这是在方差中由缺失数据带来的增量。s e m 的关键想法是：虽然m 没有一个明确的数 学形式，但它的微商d m 可以从实施e m 步的输出得到估计。这些步有效地数值地微分 复旦大学硕士学位论文6 了m ( 口) 。 详细地说，首先要求p 的m l 估计a ，然后进行一系列s e m 迭代，第( f + 1 ) 次迭代定义如 下： 输入：a 和口( 1 ) 步1 ，进行通常的e 和m 步，求出口u + 1 ) 。 步2 ，安排i = 1 ，计算 口( ( 1 ) = ( 岛，鼠一，碰“，反+ 。，以) 它是口中令第i 个分量等于彭“。 步3 ，将们( 1 ) 作为p 的现在的估计，进行e m j g i 代，求得和+ 1 ) ( ) 。 步4 ，求出此值 r ；) = 警半一，d 步5 ，重复步2 至步4 ，对i = 2 ，d 输出：0 ( 1 + 1 ) 和 r ，! ，9 ：i ，j = 1 ，d ) d m ：是一o o 时的极限矩阵 ) ，当序列苔，髫+ 1 1 ，对某一f 是稳定的时候可获 得元素，这一过程中对不同的元素能使用不同的p 值。当d m 的第i 行的全部元素都 已求出时，在后续的迭代中对i 无需再重复上面的步2 至4 步。 从数值精度的角度来说，选择e m 算法的初始值作为s e m 算法的初始值对于d m 的 计算总是安全的。但是，由于e m 的初始值可能会离m l e 比较远，这就会增加一些不必 要的迭代过程。因此，一般建议采用原始e m 算法中某一迭代值作为初始值，比如，第2 次 迭代值。 由于对d m 的计算本质上是对函数的数值微分，它当然没有对函数本身的估值来的 准确，所以s e m 算法的停止准则应比e m 算法松。通常，我们采用e m 停止准则的平方 根。 s e m 算法有高度吸引力的一个特征是最后获得的。正常情况数值上是很稳定的。 另外一个有吸引力的特征是对程序和数值错误进行内在的诊断。由于y 理论上应该是 对称的，但却可能由于计算精度的不充足或是程序的错误，导致它在数值上是不对称 的。因此，通过s e m 获得协差阵的估计的不对称性是一个程序错误的指示。更进一步， 不管是否对称，v k 可能不是半正定的或者收敛到一个鞍点。 复旦大学硕士学位论文 7 2 4 e c m 算法 e m 算法的吸引点之- - ：是q ( o i o ( o ) 的求导与最大化比观测数据下的极大似然函数来 得简单。然而，在许多重要应用中m 步却是很难实施的。由于这个原因，m e n g 和r l l b i n 提 出了e c m 算法l 8 ，4 1 。e c m 算法用一系列计算上更为简单的条件极大化步( c m ) 来代替e m 的每一m 步，通过限定p 于一个特定的子空间，使得每一次的条件极大化可以获得解析 解或是比较简单的数值解。 我们称第f 次e 步迭代后的一系列c m 步为一个c m 循环。于是，e c m 的第f 次迭代由 第f 次e 步与第z 次c m 循环构成。令s 表示每一个c m 循环中包含c m 步的数目。对于8 = 1 ，s ，第z 次g m 循环的第s 个c m 步是在 吼( = 如( + ( ”1 ) 8 )( 2 7 ) 的限定下对q ( o l o c o ) ) 最大化，其中口( 件( 。一1 ) 司是在当前c m 循环的第s 一1 步找到的最大 值。当c m 循环的整个s 步都完成的时候，我们令0 ( 1 + 1 ) = 卵+ s s ) ，然后继续第l + 1 次迭 代的e 步。 显然，设计一个有效的e c m 算法的关键是选取适当的限定条件。通常，我们很自 然地将0 分割为s 个子向量，即0 = 慨，以) 。然后，在第s 个c m 步中，固定0 。外的全 部子向量对q 关于以极大化。这等同于限制函数乳( 口) = ( 0 l ，伊，l ，伊件1 ，以) 。相 反地，我们也可以固定以对q 关于以外的全部子向量同时极大化。此时，乳( 口) = 以。另 外，在每个c m 步中，吼可能对应不同的0 的分割。其他一些限制机制可根据具体问题去 构造 m 朗g 和r 曲i n 还进一步说明了e c m 算法的收敛性。当允许无限制地在0 的参数空间 上关于它极大化，函数9 的集合是空间充满的时候，在保证e m 收敛的基本同样的条件 下，e c m 也会收敛到一个稳定点 复旦大学硕士学位论文 8 第三章分组数据下的经典线性回归模型 3 1 模型假定 考虑线性模型 挑= 7 + 盯e ，i = 1 ，n ，( 3 1 ) 其中，日独立同分布且密度函数为，盯 0 ，x i $ 1 7 为p 维向量。假设( 一o o ，c o ) 被事先 确定的分割点 t a 分为个区间阮一l ，t j ) ，1 js 膏，七2 ，其中t oe o o ，t k 三 + o o ，m ，五) = ( 一o 。，五) 。瓠的确切值不能被观测到，只知道啦玑 玩，其中啦和以为 事先给定的常数 乃 。进一步假定，当1 i n l 时，一o o = 啦 圾 o o ，即玑落入 区间( 一o 。，噩) ；当n 1 + 1 i n 2 时，一0 0 a i 玩= 0 0 ，即鲥客入区间( t k l ，+ o o ) ； 当几2 + 1 i s n 时，- - 0 0 m o ) ； ( i i i ) 矩阵x = ( x l ，) 7 满秩。 3 2m l e 的存在唯一性 采用基于似然函数的估计方法m l e 来获得回归系数的估计之前，我们首先来研 究m l e 的存在唯一性问题。由于，y ) 与0 = ( 卢，h ) 之间为一一映射，故7 ) 的m l e 存 在唯一条件等价于口的m l e 存在唯一。于是，我们只需要讨论口的m l e 是否存在唯一。 励 一碱p 一搬 p 。一 p 一慨 f 矗 所 一 阮 缸一f m ：l 复旦大学硕士学位论文9 对m l e 的研究本质上是对似然函数极值点的研究。如果似然函数本身具有良好的 性质，将有助于简化问题。函数的凸性能起到这样的作用，故首先研究l ( 口) 在什么条件 下具有凸性。根据假定( i ) 以及下面的引理3 2 1 ，可知我们所要研究的似然函数是具有 凸性的。 引理3 2 1 如果z 凹，具有上凸性，则一f ( 口) 具有下凸性。( p r a t t ( 1 9 8 1 ) 9 1 ) 注1 对于几个常用的分布它们的密度函数具有上凸性。比如： 正态分布，( ) = 了1 互；e 一譬； 极值分布，f ( x ) = e 印( 一e “) ； g a m m a ：分布，f ( x ) = 凹。一1 e - zz 0 ，a l ； w e i b u u 分布，f ( = 1 一e 印( 一如o ) ，z 20 ，0 0 ，a 1 ； p a r e t o 分布，f ( z ) = 1 一x - 4 ，z l ，n 1 ： b e t a 分布，( 。) = 凹4 1 ( 1 一z ) 6 1 ，0 z 1 ，口1 ，b 1 等等。 注2 若一2 ( p ) 严凸且口的m l e 存在，则口的m l e 存在唯一。当t o g ( z ) 严凸时，l o g f ( a ，b ) ( 一o o 口 n 2 时，p 的m l e 存在的充要条件为不存在0 o 使得下列两条件 同时成立： ( i ) 当l s i n 1 以及砌+ 1 i n 时， ，溉一卢o ；( 3 4 ) ( i i ) 当，1 1 + 1 i n ， 7 池一p 0 ( 3 5 ) 证明 - t ( o ) 的下凸性以及，的连续性可推m - t ( o ) 的下半连续。规定h o 时，- t ( o ) = o o ，则一f ( 口) 是尼件1 上的真闭下凸函数。 ( 必要性) 令b 为一非空有界集。假设存在o o b 使得一? ( 口) 达最小，且存在如0 ，使 得( i ) 、( i i ) 成立，又因f 的单调性，故对任意的k 0 ，有 - t ( o o + k o ) = 一j ( 岛+ 筇，h o + k h ) n 1 = 一l o g f ( - c o ，h o b i 一岛+ 蠡( 地一p ) ) i = 1 n 2 一l o g f ( h o a i - e , 届o + 七( 碱所，o o ) = n 1 + 1 复旦大学硕士学位论文1 0 n 一x o g f ( h o 啦一z 岛+ 七( 7 池一卢) ，h o b i 一风+ 七( 碱一卢) ) i ：n 2 + l - _ f ( o o ) ， 因此，- z ( 0 0 + k o ) 关于蚓e 增。如果i ( o o + k o ) = f ( p ) 。则如+ 枷b ，令k 一，得 出b 无界。如果z ( o o + k o ) t c o ) ，则0 0gb 。两种情况都将导出矛盾，故m l e 不存在 ( 充分性) 参考s i l 、砷l l l 】_ e 和b u r 耐g e ( 1 9 8 6 ) 【1 0 】的证明。 接下来讨论礼= n 2 时，口的m l e 的存在情况。 定理3 2 2 当n l = n 或砌一n l = n 时，口的m l e 不存在或有无穷多个。 证明 当扎1 = n 时， n l f ( p ) 一l o g f ( - o o ，她一励 i = 1 若b = s u p z ：，( z ) o = + ，则z ( o ) 0 ，但是1 i m 南一。t ( o ) = 0 ，故f ( p ) 的最大值点 不存在。 若b 6 ，i = 1 ，n ，即郦 0 ，线性方程组x 卢= ( ，皿- b - e ) l 有解。因此，满足x p ( 晒一b ) l 的0 就 有无穷多个，故f ( 口) 的最大值点有无穷多个。 当他一n 1 = n 时，同理可得 定理3 2 3 当0 n i 2 时，口的m l e 不存在。 证明 当0 n l z 时，乃死一t ，c s q 式与c s 劫式不能同时为。，于是方程组 薹三：无解， 故“p ) 的最大值点不存在 定理3 2 4当0 n l n 2 = n ，分组数k = 2 ，p 的维数p = 1 时，口的m l e 有无穷 多个。 证明 当n = n 2 ，k = 2 ，k = 1 时， l ( 0 ) = n l l o g f ( 一，蝎一所+ ( n 2 一n 1 ) l o g f ( 7 灯1 一p ，o o ) ， 令z = f ( 一o o ，l 冗一卢) ，则 f ( 口) = n l l o g z + ( 耽一n 1 ) t o g 当z = 磬时，z ( 口) 达到最大只要口满足f ( 一，晒一所= 嚣，则口就是z ( 口) 的最大值点， 于是2 ( 口) 的最大值点有无穷多个。 定理3 2 5 记观测值落入区间阮一l ，乃) 的样本点数目为，兰l 唧= n 。n l ， 讯中不为零的个数为m 。k ，6 ) 和嗡，呸) 分别是头两个观测值个数不为0 的区间。当p 的 维数p = l 时，口的m l e 存在唯一的充要条件为下列条件之一满足： ( i ) m 3 - ( i i ) f i r , = 2 ，嵋 一o o ； ( i i i ) m = 2 ，6 n 耋且嵫 他；( i i ) 存在1 i l ，如n ， 使得一 晚。 啦。 0 0 。不妨假设存在0 0 使得( 3 4 ) 和( 3 5 ) 成立，则有m 。一岛 0 ， 一岛s0 净硒f ，2 风胁虹，产生矛盾，于是由定理3 2 1 ，口的m l e 存在。另外， 王静( 2 0 0 3 ) 【2 】证明了当，( 功= 了1 磊e - 譬时，- l o g l ( o ) 严t & 。因此，m l e 不仅存在而且 睢一 ( 必要性) 上述条件都不满足= 争f l = n 2 。由定理3 2 2 和定理3 2 4 ，m l e 不存在或不唯 3 3 利用e m 算法获得m l e 的近似解 定理3 3 1分组数据下的回归模型，由e m 算法得到参数的m l e 近似解的迭代公 式为： 争e筹端一h(of(h(ot-xo)确9 复旦大学硕士学位论文 1 2 南= 一：喜序筹筠莉翟筋器厕氓 ( 3 1 0 ) 其中，( 卢( “， ( ) ) 为第f 次迭代解。 证明 使用e m 算法的关键就是要确定适当的潜在数据，以便简化计算。在分组 数据的条件下，只知道y l 落入区间k ，以) ，i = 1 ，n ，将此观测结果记为l k 。如果 将，名妇= ( y l ，) 视为潜在数据，则。包含了，k 的全部信息。确定了潜在数据 后，我们只需要去完成e 步与m 步 首先完成e 步。 由于知道k 胁，也就是挑为一个确定值，所以添加后验分布为， 对( 3 1 1 ) 式取对数后，得 n ，伊，危1 “，) = h f ( h y t 一卢) ， ( 3 1 1 ) b l l o g f ( 3 ，叫x 腑，y 咖) = 1 0 9 f ( h y 一z ：卢) + l o g h ， ( 3 1 2 ) i = l 将上式代入( 2 1 ) 式，得 n q ( 卢， i 卢“，九，1 ；m ) = e 1 0 9 f ( h y i 一励i 卢( 1 ) ，h ( 0 ，】；胁】+ n l o g h ， ( 3 1 3 ) b 1 记条件密厦函数为 舭) - ，慨圳“川= 丽粤高岛， 将丘( t ) 代入( 3 1 3 ) 5 ，得 q ， l o ，l “，k b ) = l o g ，( 船一伪厶( t ) d 亡+ 竹1 0 9 ， ( 3 1 4 ) 于是就完成了期望步。 接下来完成m 步。 首先将q 关于p 求导，得 筹，喜筑r 嬲剐啦， 复旦大学硕士学位论文 令鬻= 0 ，再将( p ，_ 1 1 ) 换成( m ) ，h q + 1 ) 即得( 3 9 ) 式 然后将q 关于h 求导，得 塑o h = 窑r 黼f ( h t 堍疵+ ! h ， 白。一卢) 。” 令鬻= 0 ，再将( 反 ) 换成( 卢( f + ，h ( t + 1 ) 即得( 3 1 0 ) 式。 推论3 3 1 当，( z ) = 去e 印( 一譬) 时，由e m 算法得到参数的m l e 近似解的迭代 公式为： 1 件1 = ( x x ) 一1 z e p 二由j i 蚶，y “，口“l ， ( 3 1 5 ) 盯2 = 去e 【( 一7 ) 2 i y 幽计仰，口】 ( 3 1 6 ) 其中，7 ( 0 = 箭，盯( 1 ) = 南，山= 与萨，风= 号铲，圣( z ) = f ：= o o f ( t ) d t ， e 一棚= 驾锊鬟鑫剑锄卿， ( 3 1 7 ) 琊嘞) 2 i ，) 】_ 盟笔誊云筹趔 2 口( ( z ：7 ( j + 1 ) 一，y ( o ) 【，( 一鸽) 一，( 一璐) 】 + 土l 、落鬣布l _ + ( ，y ( f + 1 ) 一，y ) 2 + 一( 0 2 ( 3 1 8 ) 利用算法所得的具体模拟结果见下一节 3 4 模拟 对回归模型的误差项服从正态分布的情况进行模拟。分别通过对大样本与小样本 下的二元回归模型以及大样本下的四元回归模型的模拟，考察利用e m 算法、s e m 算法 得到的估计的有效性及稳健性。模拟如下进行： 步l ，主观选定p ，盯，n ，一o 。= t o 乃 死 死一l 靠= + o 。，j = l ，k 。 步2 ，由均匀分布u ( o ，1 ) 产生协变量x ；由正态分布( o ，矿) 产生随机误差；由p + 白确定观测值雏。 步3 ，以旧( i + 1 ) 一讲o l i 1 0 8 为停止准则，实施e m 迭代算得m l e 的近似解( p ，子) 。 步4 ，以| | 瑙+ ”一，r i ( f 0 i i 1 0 4 为停止准则，从第5 次e m 开始实施s e m 迭代算得d m 。 复旦大学硕士学位论文 1 4 步5 ，由于此模型下俨j c 。a 凰卵= 0 - 2 x x 。故j 未= c r 2 ( x 僻) ，再将口用子替换，利 用公式。= 。( j d m ) - 1 即得渐近协方差阵。 我们下面依据以上各步骤分别对二元、四元回归模型进行模拟。 对二元回归模型，取p = ( 3 ，2 ，”。，口1 , t o，孔= 2 ，t 2 = 5 ，t 3 = 8 ，t 4 = 1 1 t s = + 先取初始值岛= ( 0 5 ，1 0 ，1 2 ) ，印= 1 0 。分别对n = 2 0 0 0 ，8 0 0 ，3 0 0 ，1 0 0 的情况进行模 拟，模拟结果见3 1 3 8 。 为了考察估计的正态性，我们对n = 8 0 0 的情况进行2 0 0 次模拟，分别画出各估计 的q q 图。 为了考察初始值的选取对估计的影响，再取初始值岛= ( 2 0 ，3 0 ，4 0 ) 7 ，t t o = 1 5 ，对7 1 , = 8 0 0 的情况进行模拟，模拟结果见3 9 - 3 1 0 。 为了考察缺失信息比例对收敛速度的影响，再取分组区间为t o = 一o o ，丑= 3 ，t 2 = 6 ，t s = 9 ，t 4 = + o o 的情况进行模拟，模拟结果见3 i i 3 1 2 。 从表图中可以看出： i 随着样本数的减少，各估计的相对偏差平均值大部分增大，相对偏差标准差全部增 大，但即使在样本较小的情况下，所得估计也是比较有效的。 2 样本量由2 0 0 0 减少到8 0 0 ，再从8 0 0 减到3 0 0 ，e m 迭代算法的平均迭代次数有所增加， 但不并显著。然而，当样本量减少到1 0 0 时，迭代次数却显著的增加了，而且波动 也比较大，最少的迭代次数只有3 9 次，与n = 2 0 0 0 时相近；但最大的迭代次数则高 达7 4 次，远远高于n = 2 0 0 0 ，8 0 0 ，3 0 0 的情况。 3 随着样本数的减少，由s e m 算法计算得到的渐近协差阵v 对角线上的元素全部增大， 非对角线上元素的绝对值全部增大当n = 2 0 0 0 时，v 在1 0 一l o 一的精度下是对 称的；当n = 8 0 0 时，v 在1 0 4 的精度下是对称的；当n = 3 0 0 时，v 在1 0 3 的精度下 是对称的；当n = i 0 0 时，v 在1 0 一1 0 - 2 的精度下是对称的。因此，样本数的减少 使得v 的对称性下降。但是，同v 本身的增大相比，此精度下的对称应是被容许的。 4 对于每个样本量下的各次模拟，尽管随机误差的波动使我们由e m 算法获得的点估 计各不相同，但是由s e m 算法获得的渐近协差阵却是相差不多的。 5 样本数的变化对s e m 算法的迭代次数影响不大 6 初始值离真值较大时，e m 算法迭代次数有所增多却并不显著，各估计的相对偏差无 论是平均值还是标准差都没有显著差异。因此，e m 算法对于初始值的选取并不敏 感，具有一定的稳健性。但是，由s e m 算法获得的渐近协差阵的估计只有在1 0 - 3 的 复旦大学硕士学位论文 精度下才是对称的，也就是说，初始值与真值的较大偏离导致了y 对称性的显著下 降。尽管如此，我们所获得的m l e 估计还是比较准确的。因此，当s e m 迭代所获得 的估计y 对称性不够好时，我们可以根据m l e 估计选择与其接近的初始值重新进行 一系列的迭代。 7 缺失信息比例增大，e m 算法迭代次数显著增多，即收敛速度减慢；s e m 算法的迭 代次数亦有所增多，但可以通过改变s e m 迭代中实施e m 迭代的起始位置来得到改 善。各参数估计相对偏差无论是平均值还是标准差都有增，但与真值相比，仍是比 较有效的估计。由s e m 算法获得的矿相较于分组数为5 组的情况对称性有所下降，但 相较于y 本身的增大，此精度下的对称性也是可以被接受的。 8 q q 图上的各点基本落在一条直线附近，这说明了此估计是具有正态性倾向的。这 与我们算得的y 在较高精度下的对称是一致的 对四元回归模型，取p ；( 3 ，2 ，7 ，- 6 ，一8 ) ，盯= 1 ，t o = - - o o ，噩= 一1 1 ，t 2 = - 7 ，t 3 = 一3 ，t 4 = 1 ，t 5 = 5 ，t 6 = 9 ，乃= + o o 取初始值岛= ( o 5 ，1 0 ，1 2 ，4 ，8 ) ，a o = 1 0 ，n = 8 0 0 。 模拟结果见3 1 3 - 3 1 4 。与3 5 - 3 6 比较可知： 1 随着维数的增多，迭代次数与相对偏差有所增大，但是估计的效果仍是不错的。 2 随着维数的增多，s e m 算法得到的渐近协方差在精度稍微减弱的条件下基本上仍是 对称的。 复旦大学硕士学位论文1 6 表3 1 ：n - - 2 0 0 0 初始值岛= ( 0 5 ，1 0 ，1 2 ) ，印= l o 的1 0 次e m 模拟结果 表3 2 ：n = 2 0 0 0 初始值岛= ( 0 5 ，1 0 ，1 2 ) ，o o = 1 0 的1 0 次s e m 模拟结果 n 0 8 岛 240 * 4 2 5 3 5 5 - 0 0 0 4 2 8 40 0 0 2 4 0 60 0 0 x ； 0 2 - 0 0 0 9 4 - 0 0 0 彭2 3 s0 6 蛳- 0 0 0 5 4 9 2 - 0 0 0 5 2 6 2 3 22 - 0 0 0 5 3 8 00 , 4 2 9 2 5 20 0 0 4 7 9 7 - 0 0 18 0 0 4 s “0 0 0 0 1 9 6 - 0 0 0 5 4 9 00 0 1 0 8 0 6 0 0 0 0 3 7 3 如4220 0 0 0 0 2 10 0 0 1 7 9 40 4 2 3 1 1 0 - 0 0 0 2 2 4 20 0 0 0 2 0 50 0 0 4 2 3 7 - 00 0 5 2 6 70 ，0 0 0 3 8 300 1 0 0 4 0 # l t 晶d0 0 0 0 00 6 7 4 911 3 5 30 ，0 1 3 5 0 5 几1 1 0 0 500 0 0 00 s 4 3 30 0 0 2 7 3 2 如09 7 1 80 0 0 0 010 7 5 0o 0 0 3 s 0 5 0 0 1 1 1 4 00 0 0 6 7 8 70 0 0 0 0 5 20 0 0 0 0 5 2 0 0 0 e 0 0 60 0 0 2 7 0 20 0 0 0 0 1 5 00 0 0 a ) 7 5 0 0 0 6 6 7 200 1 0 0 8 600 0 0 0 5 20 0 0 0 0 4 6 00 0 0 0 5 30 0 0 0