（运筹学与控制论专业论文）异方差截尾线性回归模型参数mle的em算法.pdf

摘要 在工业试验中，由于试验条件的限制，试验响应取截尾数据常常是必须 的在截尾数据线性回归模型中，标准的回归技术不再适用一种很自然 的方法就是把截尾时间当做实际失效时间，然后用标准的回归方法建模， 进行分析但是，忽略截尾信息可能导致错误的判定，因为未观测到的失 效时间可能与截尾时间存在很大的差异因此需要一种合适的方法来进行 分析，极大似然估计( m l e ) 因为其有许多优良的性质而被采用 a i t k i n ( 1 9 8 1 ) 给出了同方差右截尾线性回归模型中基于e m 算法的参数 m l e 的方法，并且证明e m 算法和迭代最小二乘方法在除了误差方差的估 计上不同外，其本质上是相同的，且由e m 算法的性质保证了迭代的收敛 性但是对截尾数据用似然方法的一个问题是m l e 可能不存在s i v a p u l l e 和b u r r i d g e ( 1 9 8 6 ) ，h a m a d a 和t s e ( 1 9 8 8 ) 给出了同方差截尾线性回归模型中参 数m l e 存在的充分必要条件，并且此条件可以转化为一个线性规划的求 解问题，从而可以用线性规划的方法和结果来判断其存在性对于异方差 截尾线性回归模型参数m l e 的存在性和参数估计计算的研究，文献中还很 少 本文考虑异方差截尾线性回归模型中参数m l e 的存在性问题，给出了 正态情形下m l e 存在的充分必要条件，并且同样这些条件可以转化为线 性规划解的存在性问题之后给出了用似然方程的方法求参数p ，的估计 卢，砖的两步程序，接着用e m 算法来给出参数的m l e ，证明了e m 算法与 此似然方程法的一致性，由e m 算法的性质保证了迭代的收敛性同时证 明了e m 算法给出的回归系数的m l e 与广义最小二乘估计是相同的，但对 于方差的估计不同，还考虑了方差估计的偏度问题最后通过模拟试验说 明估计的精度和截尾数据的截尾时间有关 关键词：截尾数据；回归模型；极大似然估计；e m 算法；迭代；广义最 小二乘估计 a b s t r a c t i ni n d u s t r ye x p e r i m e n t s ，b e c a u s eo fr e s t r a i no fe x p e r i m e n tc o n d i t i o n ，e x p e r i m e n t r e s p o n s ei so f t e nc e n s o r e dd a t a i nc e n s o r e dl i n e a rr e g r e s s i o nm o d e l s ，s t a n d a r dr e g r e s s i o nm e t h o d si sn o ta p p l i c a b l e i ti sao r d i n a r ym e t h o dt h a tt r e a t i n gt h ec e n s o r e dd a t a a st h ef a i l u r ed a t a ，t h e ni t i sa n a l y z e db ys t a n d a r dr e g r e s s i o nm e t h o d sb u ti tm a y r e s u l ti nw r o n gd e t e m i n e ，h e c a u s et h e r er n a yb eag r e a td i f f e r e n c eb e t w e e nc e n s o r e dt i m e a n dt r u ef a i l u r et i m es oi ti sn e c e s s a r yf i n d i n gam e t h o dt oa n a l y z ed a t a m a x i m u m l i k e l i h o o de s t i m a t ei sa d o p t e db e c a u s eo fi t sg o o dp r o p e r t i e s m a x i m u ml i k e l i h o o de s t i m a t eo ft h ep a r a m e t e r si n l i n e a rm o d e l sw i t hr i g h tt e n s o r e dd a t ai si n t r o d u c e db ya i t k i n ( 1 9 8 1 ) u s i n gt h ee ma l g o r i t h m a n dp r o v et h a te m a l g o r i t h ma n di t e r a t i v el e a s ts q u a r e sm e t h o d sa r ee s s e n t i a l l yt h es a m ee x c e p tf o rt h e e s t i m a t eo ft h ee r r o rv a r i a n c e 0 - 2 i t e r a t i v ec o n v e r g e n c ei so b t a i n e db ye ma l g o r i t h m p r o p e r t i e s b u tt h e r ei sap r o b l e mu s i n gm a x i m u ml i k e l i h o o de s t i m a t e i t i nc e n s o r e d d a t am o d e l s ，t h a ti sm l em a yn o te x i s t e n c e f o rc e n s o r el i n e a rr e g r e s s i o nm o d e l sw i t h s a i u ev a l i a n c e s ，as u f f i c i e n ta n dn e c e s s a r yc o n d i t i o ni so b t a i n e df o rt h ee x i s t e n c eo f m l eb ys i v a p u l l ea n db u r r i d g e ( 1 9 8 6 ) ，h a m a d aa n dt s e ( 1 9 8 8 ) a n dt h e s ec o n d i t i o n s m a yb et r a n s f o r mt ol i n e a rp r o g r a m m i n gp r o b l e m i ti sf e wt oe x i s t e n c eo fm l ea n d c a l c u l a t i o no fe s t i m a t eo fc h ep a r a m e t e r si r tc e n s o r e dl i n e a rr e g r e s s i o nm o d e l sw i t h u n e q u a lv a r i a n c e si nr e f e r e n c e s i nt h i sp a p e r ，w ec o n s i d e rt h ee x i s t e n c eo fm l ei nc e n s o r e dl i n e a rr e g r e s s i o nm o d e l sw i t hu n e q u a lv a r i a r m e s as u f f i c i e n ta n dn e c e s s a r yc o n d i t i o ni so b t a i n e df o re x i s t e n c e o fm l ei nn o r l l r a l ，a n di tm a yb et r a n s f o r mt ol i n e a rp r o g r a m m i n gp r o b l e mt o ot h e n w ep r o p o s et w os t e p sp r o c e d u r eo fo b t a i n i n gp a r a m e t e r l s ( 卢，口；) e s t i m a t e s ( 3 ，鲆) b y l i k e l i h o o de q u a t i o n ，o b t a i nm l e o fp a r a m e t e rb ye ma l g o r i t h m a n dw ep r o v et h e c o n s i s t e n c yo fe ma l g o r i t h ma n dt h i so r d i n a r yi t e r a t i v ea l g o r i t h m ，t h u se n s u r ei t sc o n v e r g e n c ef r o me ma l g o r i t h mp r o p e r t i e sa tt h e s 8 l n et i m e ，w ep r o v et h a tm l eo fr e g r e s s i o n c o f f i c i e n t sa r ei d e n t i c a lt og e n e r a l i z e dl e a s ts q u a r e se s t i m a t e lb u tt h o s ef o rt h ev a r i a n c e e s t i m a t e s8 ed i f f e r e n tt h eb i a s e so ft h ev a r i a n c ee s t i m a t e sa r ed i s c u s s e da tl a s t ：w e d e s c r i b et h a te s t i m a t ea e c l l r a c yi sr e l a t e dt oc e n s o r e dt i m eo fc e n s o r e dd a t a k e y w o r d s ：c e n s o r e dd a t a ；r e g r e s s i o nm o d e l ；i n a x i l n n l t ll i k e l i h o o de s t i r n a t e ；e ma l g o r i t h m ；i t e r a t i v e ；g e n e r a l i z e dl e a s ts q u a r e se s t i m a t e 引言 引言 问题的背景 在改善可靠性的试验中，由于时间，费用，环境等条件的限制，常常得到的只能 是截尾数据对于当今高度可靠的产品，可能很难在一个合理的时间段内甚至是在加 速强度的条件下观测到试验失效因此也不可能用一种专门的监测装置来记录每一个 试验单元的失效时间，那么可以周期性地检查每个单元是否失效，这样得到的就是所 谓的区间截尾数据，左截尾和右截尾数据都是它的特殊情形 虽然截尾数据容易得到，但是它相对于完全数据来说包含的信息少，很难对其进 行分析在早期的工作中，或者是把截尾数据看作实际失效数据，或者是直接忽略掉 截尾数据这样的方法虽然简单易行，并且在有些情形下也得到了可行的结果，但是 显然它有很大的缺陷因此，就需要有一些方法来分析截尾数据，在截尾数据线性回 归模型中，极大似然估计因为其有很多优良的统计性质常常用来估计其未知参数 s c h m e e 和h a h n ( 1 9 7 9 ) 提出了同方差右截尾线性回归模型中正态分布条件下用 迭代最小二乘来估计未知参数的一种迭代算法，这种方法虽然简单易行，甚至容易被 非统计人员接受和应用，但其明显存在迭代的收敛性问题，且其估计的精度问题也较 极大似然估计差 a i t k i n ( 1 9 8 1 ) 接着给出了这个模型中基于e m 算法的求参数的极 大似然估计的方法。并且证明e m 算法和迭代最小二乘方法在除了方差的估计上不 同外，其本质上是相同的，且由e m 算法的性质保证了迭代的收敛性 s i v a p u l i e 和 b u r r i d g e ( 1 9 8 6 ) ，h m n a d a 和t s e ( 1 9 8 8 ) 给出了同方羞截尾线性回归模型中参数m l e 存在的充分必要条件 对于异方差截尾线性回归模型的研究，文献中还很少，但是在实际中这类模型常 常能遇到比如在线性回归模型中，若观测向量是1 2 件样品的某项指标，它们是在k 台仪器上测试得到的不妨设前n ，个是第一台仪器上测得的，接下来的毗个是从第 二台仪器上测得的，依次类推那么模型显然可以写为： u = z ；卢+ e i j ，i = 1 ，2 ，k ；j = 1 ，2 ，啦 其中玑，表示在第i 台仪器上测量的第j 个样品的指标，z 。表示对应的样品上自变 量的取值，e 。j 是对应的试验测量误差假设所有这些测量都是相互独立的，但不同 的仪器测量精度不同。即对于测量误差e m 下面的假定是合理的： v a r ( e i ) = 口。2 ，i = 1 ，2 ，；j = 1 ，2 ，- - ，“ 2异方差截尾线性回归模型参数m l e 的e m 算法 c o v ( e i j ，e i 7 ) = 0 ，i i ，j j 在试验设计中，这类模型也是一种很重要的模型，比如在区组设计中，不同区组中的 数据显然是异方差的因此，考虑此模型非常有意义 同方差数据情形近年来的工作 假设y ，y z ，y 。是经过变换后相互独立的响应变量的观测值，且服从正态分布， 则同方差右截尾线性回归模型可表示为： f 玑2 地十口矗，1 i n 2 k i n ( o ，1 ) ，矗一l i d ( 1 ) l 胁= z ：卢= e 岛。奶x i 0 ；1 、j = 0 在( 1 ) 式中口 0 ，z ：和卢是p 维向量，矗为随机误差，r l 为试验次数，( p 口) 为未 知参数观测值y 是精确值，1 i m ；一o o 0 ，( p ，o - 1 ，o - 2 ，( 9 - m ) 为未知参数，5 。，为随机误差， 总的试验次数n = 篓1 批且一o o a js b 玎 0 ，目c 0 当f _ l o g l 时，定理21 仍然是有效的并且他们指出其定理的条件只是一个充 分条牛，其中条件( i ) 保证了参数m l e 的存在性但对( 5 ) 式中参数m l e 的存在性 问题一直没有研究 本文在h a r t l e y 和j a y a t i l l a k e ( 1 9 7 3 ) 的基础上，考虑异方差截尾线性回归模型( 5 ) 的参数m l e 的存在性问题，给出了正态情形下参数m l e 存在的充分必要条件，且 引言 这些条件可以转化为线性规划可行解的存在性问题之后给出了用似然方程求参数 卢，呒2 的估计p ，霹的两步程序，接着用e m 算法来给出参数的m l e ，证明了e m 算 法与此似然方程法的一致眭由d e m p s t e r ，l a i r d ，r u b i n ( 1 9 7 7 ) 和w u ( 1 9 8 3 ) 文献中的 结论保证了迭代的收敛性并且e m 算法给出的回归系数的m l e 与广义最小二乘估 计是相同的，但对于方差的估计不同，还考虑了它们估计的偏度问题最后通过模拟 试验说明估计的精度和截尾数据的截尾时间有关 5 6 = 一量童耋垄垒些竺皇望堡型圭塞丝堡塑垒丝苎童 第一章m l e 的存在性 1 1 m l e 存在的充分必要条件 考虑模型( 5 ) 中，对任意ls i m ，每组观测值蜘包含如下四种情形： ( a ) 一o o = a 0 驹sb 日 。，1 jsr 1 ( 左截尾) ； ( b j o 。 o 巧s sb o = o 。，尼l + 1s js 丑2 ( 右截尾) ； ( c ) 一o 。 0 记舻为p 维欧氏空间 记2 ( 日) = 一l o g l ( o ) ，参数0 的m l e 为0 + ，则其等价于0 + 使f ( 日) 最小化，即 3 = 目+ ：f ( 矿) = 啷n l ( 口) ) 是非空有界集 定义条件：对一切lsj 茎m ，存在非零的卢p ，使得 ( ) z ：，卢0 ，15j r 。】； ( ，) 。：，卢0 ，局1 + 1 j 最2 ； ( i i i ) z , f f 3 = 0 ，风2 十1sj 批 定理1 ：假设m 置，i = 1 ，2 ，m ：令q2 ( 。，y ， ，y 。_ 一，。) ： 存在卢使得z 台卢= ，i = 1 ，m ，j = r + 1 ， ，批 ，那么除了集合q 以外，参数 m l e 存在的充分必要条件是p 不满足 证明：先证其必要性： ( 反证法) 假定存在一个卢0 ，卢f p 满足，即： ( i ) z ：，卢0 ，1 j 皿】； ( 【i ) z ：，卢0 ，r 。1 + l j 咒2 ； ( i i i ) z ：，= 0 ，咒2 + 1sj 批； 则对任意的k 0 ，由正态分布的分布函数西的性质及口+ = + ：m ，怫， ， 新) 为参 数0 的m l e ，有： m l e 的存在性 ( i ) _ 一l o g ( i ) ( h ；( b 。- 一。；【卢+ + 卢j ) j 一l o g , t , ( h ：( b i j 一。0 + ) ) ，1 j 曼忍1 ； ( i i ) 一f 。g ( 1 一西( k ( 一z o ( 卢+ + 舢) s 山9 ( 1 一中( k ( n 甜一z 扩) ) ) ，r - + 1 j r 2 ； ( 1 1 1 ) 一l 。9 ( 西( 埘( 6 。一z o ( 卢+ + k 刚一西( k ( 一z ；( p + + k 卢) ) ) ) 一f 。9 ( 辩) 一f o g 曲i h ；( y t j 一 。o ( 卢+ 十即) ) ) = - f 0 9 ( 巾( m ( 岵一z 扩) ) 一面( 埘( 一z 扩) ) 一? ( 蟑) - f 0 9 ( h 池，一 z d t 卢f + ) ) ，r 。2 + 1 js 胍； 从而 2 ( p + + 女卢： ：，- - ， k ) = - f o g 壬( k ( 一z 妒十k 卢) ) ) + 登 乩，( 叫 * ( a i j - x ；( 外) ) ) + 卜2 叩( 1 一m ( ；( 卢+ + 卢) ) ) ) l + - l o g ( ( h * ( b 0 - - x o ( 卢4 + 俐一西( k ( 一z o ( 卢+ + 即) ) ) ) + _ f d 9 ( k ) _ f 0 鲫( m ( 如一z 妒+ 印) ) ) ) m ，且l ，、 - f o g r p ( h ；( b 订一z 矿) ) l i = 1 j = l + 登_ l o g ( z 叫丘( a ；j - - x 别) 见 ， + l f 0 9 ( 中( ：( 6 ，一屹_ 臼+ ) ) 中( k ( n i j = r 2 + l + 乩9 ( 埘) _ f 。g 庐( k ( 一。) 】) j = 忍+ i = f ( 矿) = 2 ( 卢+ ，k 一，) 如果f ( 矿q - k 卢，埘，h i ，7 嘞) = f ( 矿，啦， 新) ，则说明b 是无界的；如果 f ( 卢+ + 卢， j ， ；， h ) 0 ，存在一个区域nce ，在区域q 的外部有上 0 ，存在一个常数c 0 ，使得 知z 矿学) ( 卢， ) ：l ( z ， ) ) ( 卢h ) ：iz o g h d c ，i 一1 ，2 ，一，m ) 那么在集合h = h ：i 1 0 9 h , ：l c ，i = 1 ，2 ，一，m ) 上，当s + 不满足，则对于至少 一组1 tsm ，对任意的卢，一个或多个以下条件满足： ( n ) z ：k 卢 0 ，对于至少一个，1 k 茎r n ； ( b ) z 磊卢 0 时，由i l 卢| | 一。，| lz ；卢i i o ( 3 ，有 l i r as u pl ( f l ，h ) 怕_ | 一。h h 茎l i r as u p 圣( 也( 6 【一o 。h l l i r a 壬( d o ( 6 n i 俐l 一。 册m 商 ，jl 鬻 吕8 m 蚶 m 吕 l h 一 卵 ， m 蚶 m ：星卅喵 、f、l， 卢 吒 物 “ m 埘 m z 玑 e “ e m 蚶 m 消 z m l e 的存在性 情形( b ) 当至少有一个k 满足z 氛卢 0 时，由| i 卢| | 一。，l i 。0 卢| | 一 。，知 9 l l 口l l i l r a 。h s u h pl ( 卢，h ) 刮| i 慨r as u 川p ( 酬b t k - - x 叫姒叫圳i = 1 ，旦。蚋一驯) l i 微( 喇叫洲叫d 3 ( ( 2 t k - - x 冽) 驵m ，璺，螂( e “( 妨一z o 叫 儿 其中e 一。茎d o ，d 1 ，d 2 ，d 3 e 。是常数 当r 。2 + 1s ks r ，。曩卢 0 时，由i | 口| | 一。，l iz o pl l o 。，知 怕l i ms 磐u p ( 鼬) 叫l i ms u pf ( 去e 川舻) 2 2 ) 疆m ，童，( ( 曼| | 撬( 一去8 七1 1 m k 川2 归) 疆1 7 i d e ( m n i e 弋一 一z ) 。；卢) ) 当r + 1s 女n g z o 肘，同样有 怕l i ms u 。p 聊) 叫 i r as u ps f ( 、去8 “豇叱脚2 ) g i ih # ( m 姒q 叫) ，札 溉( 一去8 七1 q m q k 绑) 疆m ，曼螂( e _ c ( 物q 卢) ) 1 n中 h 一 二r l rl d 。h d 日 u u s h s h 昌8 暑8 0 ，存在常数c , d ，使得 ( 卢h i ，h m ) ：l ( z ，h i ，h m ) ) c “卢：h i ，h m ) ：| | 卢l l sd ，l l o g h i c ) 即当s + 条件不满足时，2 嗡l ( 目) = 0 ，参数的m l e 存在证明完毕 定义p + 条件t 如果对某一组观测1 tsm ，n t r t ，存在卢，使得 ( 1 ) a t j z ：，卢s0 ，当1 茎j 墨i t t l ； ( 2 ) b t j z 0 卢0 ，当r l + 15jsr t 2 ； ( 3 ) 。“z 0 卢曼b t j ，当r 挖+ 1sj 曼r t ( 4 ) y t j = z i 卢，当r c + 1 j j 则有如下定理： 定理2 ：如果s ”条件满足，则m l e 不存在 证明：假设存在卢满足条件p + ，那么对任意的h ? ( 0 ，o 。) ，i = 1 ，2 ，m ，任意的 k 0 ，有： ： 釜 - 1 0 9 西( ( n ? 撕。) ( b , j - z 圳 + r t 2 一f 。9 ( 1 一中( ( ? + 。) ( a t j - - x k ho p ) ) ) + 卜f d 9 ( 1 一中( ( 醒+ t) ) l + l 一2 。9 ( 西( ( ? + k h 。) ( 6 c j z ) 一圣( ( ? + k h t ) ( n “一z 0 声) ) ) f 十 一f 。9 ( ? + 女) 一2 。9 ( o ) ) j = r t + 1 m r l 十 乩。中( n 一z 拗 + 曼 一l o g ( 1 叫帆，一锎) + j 一一母( 碍( 一z ；声) ) ) i + 山。( 西( n 一。拗一中( n 一z ) + 山，( h ? ) 乩。( 一z o 酬) m l e 的存在性 考虑上式中第一部分 ) 中的内容，因为对1 曼t m ，l j r ，b t j z 0 卢 0 ，a t j - - z ：，声0 ，一t o a 曲( o ) 为常数，而第二部分对给定的卢，h o ，螋， 0 _ 1 ， l ， 为定值，所以，当k 0 0 ，2 ( 口础，明， 1 1 ，h + k h c ， o + ， ) 一一o 。故b 为空集或无界集，m l e 不存在证明完毕 1 2 定理中条件的验证办法 定理1 给出的s + 条件在实际中不容易判断，因此本文给出了基于线性规划的方 法，其很容易来验证条件是否满足，从而来判断参数m l e 的存在性 在s + 条件中令玩= ( z ：，) 为( 咒l m ) p 维矩阵， 1 zsm ，1 j r i l ，z 1 = ( z ：，) 为( ( 兄2 一忍1 ) m ) p 维矩阵，1 ism ，r l + 1 j r m 邑= ( 蜀一z 1 ) ，z 3 = ( z ：，) 为( ( 肌一砬2 ) m ) x p 维矩阵，i = 1 i m ，咒2 + 1 j 肌， 非零向量”= 卢，则s + 条件可以表达如下： 存在 0 ，使得易 s0 ，z a v = 0( 1 2 ) s i l v a p u l l e 和b u r r i d g e ( 1 9 8 6 ) 的文章中给出了用线性规划的方法来决定是否存在 一个方向满足( 1 2 ) 式本文同样可以用它来处理，简述如下：假设z = ( z l ，蜀) ，且 秩r ( z ) = p ，如果r ( z a ) = p r ，那么这个问题就可以简化为： 寻找= ( w 1 ，w 2 ，w ，) 0 ，使得c w 0( 1 3 ) c 为其相对应的矩阵，这样就可以用线性规划的方法来处理这个问题 首先引入一非零的松弛变量8 来代替无约束的向量w ，使这个问题变为 寻找s 芝0 ，使得d s 之0 如果存在0 ，满足( 1 3 ) 式，则存在松弛变量s + 0 ，使得( g 中i 代表单位矩阵又因为c 是r 满秩的，则上式可以写为： c 1 1 1 0 忖s l 。 ( 1 4 ) 州s + ) = 0 ，其 其中c 1 为非奇异的r r 矩阵，s ：为r 1 向量求解( 1 , 5 ) 式，有 w = 一e i s j ，s ；= 岛g i ls i ( 1 5 ) f 16 1 1 2 异方差截尾线性回归模型参数u l e 的e m 算法 在( 1 4 ) 式中，如果令d = c 2 c ；1 ，s = s i ，则( 14 ) 与( 1 3 ) 式相同除了平凡的情形， ( 14 ) 式转化为的线性规划问题为：m i n e s i = s 7 c ，满足约束条件c s 一b ，s 0 ，b = ( 1 ，1 ) 7 ，则它用标准的线性规划的方法很容易来判断线性规划的解是否存在故若 线性规划问题存在解，则本文所考虑模型的参数m l e 不存在 同样定理2 的s ”条件可以转化为线性规划的求解问题 参数的估计 考虑模型( 5 ) ，似然函数为 第二章参数的估计 m ” t 7 础) - 。f l 闰f i ( 西( 警( 攀) ) 瑚f i i 1 曲( 学) ) ( 2 ) 记对数似然函数为l o g l ( o ) ，则基于似然等式的求解卢，一；的极大似然估计的两步 程序如下： 第一步t 对数似然函数在正态情形下关于卢的偏导数， 酉a l o g l = 姜髂霎器端蛳+ 毒，塞。c m n t = e 。i ( 功一卢。) * i = 1 1 i = 1 其中 a i j z ：，卢 2 = 学 # 。巧) z 玎 j ( 2 2 ) 锄= z ；口+ a s ( u i j ，v i i ) j = 1 ，兄 = j = r ：+ 1 ，批 ，= 畿蒜器 令( 22 ) 式等于零，则卢的极大似然估计可以写成一个广义最小二乘的形式，其截尾 的数据被硒+ 吼s ( ，v 4 i j ) 所代替即 p = ( 善m 孑1 墨训弋善m 嘉倒钏 ( 。3 ) 其中 咒= ( x i l ，j t i n i ) ，毛= ( z i v ，z i n ，) 第二步：对o - i 求偏导， 1 3 一o l o g l n i - r i + 至兰堡! ! ! 挲二丝芏+ 三争生堕! ) ( 堕二丝l 重塑1 2 堕二! 堡! ，仉叽 i 丁?。霄笔j 由( j ) 一西( 札“) 一( m 一毋) 一是。业业穹最渊划 一丐r 一 1 4异方差截尾线性回归模型参数m l e 的丑算法 令其偏导等于零，则有 辞= 墨r + 。( 一o ，) 2 ( 川一r ) 一墨。业趔号老骞等麓掣 ( 2 4 ) ( 2 5 ) 由等式( 2 3 ) 和( 2 5 ) 得到的迭代极大似然估计的过程如下： 迭代o ：把截尾数据的截尾时间当作其实际失效时间，与其他的精确数据构成一伪完 全数据集，基于此数据集用广义最小二乘来拟合模型记得到的参数卢，盯 的初始估 计为声( 0 1 ) 或( 0 ) 接着用声( 0 、，民( o 来计算z i ，构成新的伪完全数据集 迭代1 ：把迭代0 中得到的铴伪完全数据集代入( 2 3 ) 来计算新的卢的估计，记为 声( ”接着计算f 岛，代入( 2 5 ) 式得到以的新估计魂( ” 接下来的迭代：继续上面的迭代过程直到其收敛为止 本文用e m 算法来求m l e ，它保证了其迭代的收敛性这个算法的本质是先计算 出在完全数据( 精确值) 情形下参数的充分统计量，然后在统计量中截尾数据在每一 步迭代中用它的条件期望代替( e 步) ；这些被代替的统计量可看作完全数据，用来获 得新的参数估计( m 步) 对于正态数据情形，e m 算法的两步如下： e 步：在完全数据情形参数的充分统计量为兰，兰。z i k y 。和差，磅，则它的条 件期望为 跚e ( n i 曼b o ，卢，爵) + 蛳 。= lj = l# 1j = r + j 和 鱼m e ( 增场s ，卢，。；) + 孵 其中 e ( 场i 。玎茎k ，兰b i j ，卢，砖) = z 0 卢+ 盯。s ( 札缸，弘j ) e ( 瑁i a i j 曼b i j ，妒，一；) = ( z 苦卢) 2 + 。； 妒( u ) ( 。嵇+ z ：，卢) + 以可荔 ( 见附录证明) 兰! 堕! 堕堕星1 4 p ( u _ f j ) m 步：因为卢的m l e 和广义最小二乘估计是相同的，所以第r 步迭代巾，用z i 卢十 一 7 s ( “乳谚) 代替截尾的玑得到第r + 1 步迭代中卢的新估计“在兰，( 蚴一 参数的估计 ，! ，) 2 中用z o p ? 、+ a p s ( “，”护) 代替截尾的，用( z o p ( r ) 2 + a 护2 + 一： 盟燮搿鬟产代雠2 ，有 m一+”2=壹j=1cz。p(r，2+。!呻2+a兰生璺兰垒堕j=-誓：荔；三善篆擎 2 茁0 声7 ( z 0 卢7 + 盯：n s ( u g 令”1 n = d 巾= 砰，有 茸= 兰丑+ 。( 蚴 ( m 一最) 一墨，业趔鼍焉窘岩黯掣 则井的m l e 与前面的( 25 ) 式有相同的形式 ( 26 ) p z 一 批一 十 户卢 ，+ q z+ 谚 1 6 异方差截尾线性回归模型参数m l e 的e m 算法 第三章方差估计的偏度 本文给出的回归参数卢的m l e 和广义最j , - 乘估计是相同的，但是其对于方差 一；的估计不同广义最小二乘估计在每一步迭代中，对于砰的估计都是用其在完全 数据情形的爵的无偏估计，其中截尾的观测用条件期望来代替 ( m p ) j ；= ( z ；口+ 函s ( ：壹( 删) 2 ) 吃声) 2 + 登( y i j ，- r j = r 。+ 1 眦 + ( j = r + 1 其中声是用广义最小二乘方法得到的p 的拟合值， 由广义最小二乘估计给出的砰的估计和( 2 5 ) 式给出的m l e 是不同的，且它们 的不同在“。，一一0 0 ，一o o 的极值情形很容易比较 当札4 _ 一。，_ 。时， s ( u u ，) _ o ，塑蔓堕 亳竽篙鼍学_ 0 ，则 m j ；= 吼- 2 + ( 扯，塞。鬯j 2 j 十1 ( m p ) 霹= ( 物一z 。tp ) 2 因为声和j 是相同的，则这两种方法给出的砰的估计的期望为 嘲= ；等挈e ( ) 耋。黔箬， ：尝砖 n 。一瓦 一r ，一p m p m 一兄一p n 。一p 簪 札甜 篇 肌斟耳， 方差估计的偏度 由此可见，广义最小二乘的方法和m l e 的方法得到的砰的估计都是有偏估计， 又因为在完全数据情形，的m l e 的偏度是通过调整其自由度来校正的，所以很 自然的在截尾情形我们建议也做相应的处理，即在口? 的m l e 中用m 一丑一p 来代 替批一忌并且砰的估计的偏度与总的观测中截尾观测所占的比例有关，这一点在 下一节模拟试验中可以看到也就是说在截尾观测个数较多的时候，如果不考虑口? 的估计的偏度，那么得到的估计的偏度可能就比较大，因此需要对其偏度进行调整 1 7 1 8异方差截尾线性回归模型参数m l e 的e m 算法 第四章模拟试验 本章中我们首先通过模拟试验来说明参数估计的精度和截尾数据观测的截尾时间 有关，且它随着截尾程度的加重而下降我们考虑试验设计中的一个例子，用一个包 含5 个因子( a ，b ，c ，d ，e ) ，每个因子有二个水平的1 6 次( 2 5 。) 试验来模拟试验的设 计矩阵见表4 1 表4 1 设计矩阵 在模型y i j = z ；卢+ 吼中，i = 1 ，2 ，j = 1 ，2 ，8 ，包含两个真正的效应( a ，b ) ， 相应的，? = ( 4 ，一3 ) ，截距项对应的阮= 2 ，( 7 1 = m 5 ，0 - 2 = 2 我们分别使响应变量在 1 ，3 , 6 ，9 时截尾，傲5 0 0 次模拟模拟结果见下表42 模拟试验 表42 估计精度的比较 未知在1 截尾情形 在3 截尾情形 参数均值标准差均方误差均值标准差均方误差 岛 1 5 7 3 6 1 5 2o 4 1 5 3 6 8 20 3 5 3 9 8 9 72 0 2 7 7 9 2 80 2 9 9 2 7 5 30 0 9 0 1 5 9 0 岛 3 4 7 2 3 5 9 20 2 9 1 1 6 0 60 3 6 3 0 0 9 74 0 1 3 5 6 8 20 3 2 9 6 5 7 50 1 0 8 6 4 0 8 岛 3 0 5 9 2 7 1 203 5 4 9 9 0 00 1 2 9 2 7 9 0- 29 7 3 6 5 002 7 6 4 8 4 400 7 6 9 8 5 0 盯l 0 4 3 8 8 8 0 101 7 0 8 8 8 000 3 2 8 7 9 904 3 0 2 5 3 8011 4 8 5 3 700 1 8 0 2 9 5 盯2 17 5 5 5 2 0 40 6 8 3 5 5 1 90 5 2 6 0 7 9 01 7 2 1 0 1 5 20 4 5 9 4 1 4 602 8 8 4 7 2 2 未知在6 截尾的情形在9 截尾的情形 参数 均值标准差均方误差均值 标准差均方误差 岛 19 6 5 6 6 8 00 2 5 4 9 0 7 800 6 6 0 2 6 719 8 0 1 0 5 8 01 9 5 9 8 4 700 3 8 7 2 9 0 卢l 3 9 4 0 4 3 0 402 8 0 3 3 7 00 0 8 1 9 8 0 23 9 5 3 3 8 0 20 1 9 5 7 7 8 400 4 0 4 2 5 9 昆 2 9 8 2 9 3 9 40 2 6 4 2 3 5 10 0 6 9 9 7 1 6、3 0 8 1 6 6 2 8 02 0 4 5 3 0 200 4 8 4 1 7 7 盯1 0 4 4 7 6 0 7 20 1 0 3 2 9 3 8 0 0 1 3 3 9 3 3 0 4 8 5 7 3 2 600 9 4 2 0 5 000 0 9 0 6 0 4 盯2 1 7 9 0 4 2 8 80 4 1 3 1 7 5 00 2 1 4 2 9 2 31 9 4 2 9 3 0 40 3 7 6 8 2 0 00 1 4 4 9 6 6 3 从表4 2 的模拟结果可很明显看出，从均方误差的意义上说，参数极大似然估计 的精度随着截尾程度的加重而下降并且我们发现在模拟中，其e m 算法的平均迭代 次数在9 截尾时为5 次，在6 截尾时为6 次，在3 截尾时为8 次，在1 截尾时为9 次，从而说明此迭代算法的收敛速度与截尾时间有关 表4 3 方差估计的偏度 1 9 参 在l 截尾的情形在3 截尾的情形 在6 截尾的情形在9 截尾的情形 数 0 一引mn l n 0 0 i 1 1 , ln l n 1 0 0 in 1 n l n 1 0 一引n in l n f y l 00 6 1 1 190 5 600 6 9 7 48050 ，0 5 2 3 940 2 50 0 1 4 2 620 1 2 5 0 - 2 02 4 4 4 70 2 了8 9 802 0 9 5 700 5 7 0 6 2 0异方差截尾线性回归模型参数m l e 的b m 算法 在表4 3 中，0 为方差的估计值，0 为真值，n 1 表示截尾观测的个数，n 表示 总的观测个数 从表4 3 中可以看出方差o - ；的估计的偏度与总的观测中截尾观测所占的比例，即 截尾观测的个数有关当响应变量在9 截尾时，截尾观测个数是总观测个数的1 25 ， 在1 截尾时是5 6 ，但在1 截尾的估计的偏度是9 截尾时的4 倍多 讨论 第五章讨论 5 1 b a y e s 方法 对截尾数据用基于似然的方法的一个问题就是m l e 不存在或者得到的解不是真 正的极大值点避开m l e 可估性问题的一个很自然的方法就是用b a y e s 方法来处 理，因为参数的估计可能很大但不是无穷大通过利用合适的先验分布。能够计算出 好的后验分布，并利用它们的众数或中位数来估计参数 w u 和h a m a d a ( 2 0 0 0 ) 给出了同方差情形下基于b a y e s 方法的对截尾和失效数据 的g i b b s 抽样算法，具体陈述如下： 第一步：对( 卢，一2 ) 的当前值从，( 吼慨a 2 ) 中抽样，对于每一个区间截尾观测值( o ，6 ) ， 从条件正态密度 0 - - 1 ( 矾一口 一一1 ) 西( 2 n 6 一z ：卢) 口一1 ) 一西( 2 n o z ：卢) a 一1 ) 】，l n a 亟 l n b 中抽取样本拈 第二步t 从，( 卢，a 妒) 中抽样，其中y + = ( ，口) ，y 表示实际失效时间数据将( 卢，0 - 2 ) 的抽样值作为下次迭代中所用的当前值 经过第一步和第二步之间的充分迭代后，g i b b s 样本将收敛，使得第二步中抽取的样 本是来自( 卢，0 - 2 ) 的后验分布 我们知道b a y e s 方法的结果依赖于其先验分布的选择，w u 和h a m a d a ( 2 0 0 0 ) 文 中0 - 2 的边际分布为一个参数为( u i 2 ，a s ；2 ) 的逆r 分布，即 在给定0 - 2 的条件下卢的条件分布是一个均值为卢协方差矩阵为一2 m “的多元正态 分布，即 卢j 。2 m ( 扫，0 - 2 m 。) 其中的先验常数的选择不同都可能导致其结果的不同 对于异方差情形，我们能否给出参数( 卢，砰) 的先验分布，以及怎样来给出基于 b a y e s 方法的g i b b s 抽样这是下一步我们要做的主要的工作 5 2 模型中非正态数据的情形 对随机误差是非正态的情形，只要知道其分布，理论上是可以得到结果的，但其是 否会有本文这样的结果，可以转化为广义最小二乘方法，以及算法的收敛速度问题， 2 1 2 2 异方差截尾线性回归模型参数m l e 的e m