（应用数学专业论文）生存数据的密度函数估计问题.pdf

英义摘要 生存数据的密度函数估计问题 摘要：生存分析是近二三十年来数理统计新分支，它是根据医学，生命科学，可 靠性科学，保险科学中的大量实际问题提出的，它可以广义地认为对生存时问( 非 负随机变量) 的一类统计分析技术，主要研究随机删失数据的统计分析。随机删 失数据是生命科学，医药追踪，可靠性寿命试验及其他一些实际问题中常常碰到 的一种重要类型的统计数据，其研究理论与研究方法不仅能应用于生命科学，医 药卫生，可靠性工程，而且在保险数学，犯罪学，社会学，市场学，环境科学， 航空航天科学，等高科技领域都有广泛的应用前景。 本文以非参数核估计为出发点，在数据未删失的情况下，用模拟方法研究使 用不同的核函数和不同的光滑因子对密度函数的估计效果，得出估计效果表现好 的核函数来进行计算，并初步探讨了光滑因子的不同选择方法的估计效果。在第 一部分的最后，本文还探讨了多维核函数的估计问题，并针对二维的情况及行了 模拟。 在第二部分，本文对比了在数据删失的情况下，生存数据的估计方法与非参 核密度估计的差异，说明了未删失情况下的非参核密度估计的局限性。进而，研 究针对生存数据的不同估计方法- k a p l a n m e i e r 乘积限和“薛留根”估计 法。最后，对于生存数据最最主要的困难合适的光滑参数，即窗宽选取进行 讨论，分别对比了窗宽的经验选择法和仞步变窗宽方法，紧邻估计法与交叉核实 法，并使用a m s e ( 渐进平均平方误差) 准则来对上述四个方法进行评价，最终 综合以上研究，结合初等变窗宽法与交错核实法制定了新的核密度估计方法。同 时这也是本文的最大的创新点之一。最后在总结中提出了变窗宽核估计的改进方 向。 关键词：核估计随机删失交叉核实法非参数估计 英文摘要 t h ep r o b l e mo fd e n s i t ye s t i m a t i o no fs u r v i v a ld a t a a b s t r a c t ：s u r v i v a la n a l y z eh a sb e e nd e v e l o p e da san e wb r a n c ho fs t a t i s t i c sf o r r e c e n t l yt w e n t yy e a r s i ti sb r o u g h to u tb yt h en e e do fm a n yp r a c t i c ep r o b l e m s ，s u c h a sm e d i c i n e ，l i f es c i e n c e ，r e l i a b i l i t ys c i e n c ea n di n s u r a n c e ，e t c i tc a nb er e g a r d e da s t h es t a t i s t i ct e c h n i q u eo fa n a l y s i sl i f e t i m ed a t a s u r v i v a la n a l y z e sf o c u s e so nt h e s t a t i s t i ca n a l y s i so fc e n s o r e dd a t a r i g h tr a n d o mc e n s o r s h i pi st h em o s ti m p o r t a n t s t a t i s t i cd a t ai nt e r m so fl i f es c i e n c e ，m e d i c i n a lp u r s u i n g ，s u r v i v a lt e s to fr e l i a b i l i t y a n do t h e rp r a c t i c a lp r o b l e m s ，w h o s er e s e a r c hm e t h o d sa n dt h e o r i e sn o to n l yc o u l d b e e na p p l i e di n t ol i f es c i e n c e ，p u b l i ch e a l t h ，r e l i a b i l i t ye n g i n e e r ，b u ta l s od e a lw i t ht h e p r o b l e m so fi n s u r a n c em a t h e m a t i c s ，c r i m i n o l o g y , s o c i o l o g y , c i r c u m s t a n c es c i e n c e ， a c a d e m i ca v i a t i o ns u b j e c t m e a n w h i l e ，t h i s s u b j e c th a s ae x t e n s i v e l y a p p l i e d f o r e g r o u n d a tf i r s t ，t h i s p a p e rd i s c u s s e dk e r n e ld e n s i t ye s t i m a t i o n ，w h i c hd e a l sw i t h u n c e n s o r e dd a t a t h e nw eu s es i m u l a t i v em e t h o dt oc o m p a r ed i f f e r e n tk e r n e lf u n c t i o n a n db a n d w i d t h s r e s u l t s ，i no r d e rt os e l e c tb e t t e rk e r n e l st oc a l c u l a t e e p a n e c h n i k o v k e r n e l 。w h a ti sm o r e ，w ea l s os i m p l yd i s c u s s e dt h es i m u l a t i o ne f f e c t so fd i f f e r e n t k i n d so fb a n d w i d t h s a tt h el a s to ft h i sp a r t ，t h ep r o b l e m so fn d i m e n s i o nk e r n e l e s t i m a t i o nh a db e e nr e s e a r c h e d ，a n dw es i m u l a t e dt w o - d i m e n s i o ns i t u a t i o n s e c o n d l y , u n d e rt h ec i r c u m s t a n c eo fc e n s o r e dd a t a ，w ec o m p a r e dt h ed i f f e r e n c e b e t w e e ne s t i m a t i v em e t h o do fs u r v i v a ld a t aa n dn o n p a r a m e t e rk e r n e le s t i m a t i o n ， w h i c ha c c o u n t e df o rw h yt h ef i r s tm e t h o di sn e c e s s a r yf o re s t i m a t i n gs u r v i v a ld a t a f a r t h e rm o r e ，k a p l a n - m e i e rp r o d u c t l i m i te s t i m a t o ra n do t h e rm e t h o d sg o ti n v o l v e d i nt h i sp a p e r f i n a l l y , t h em a i nd i f f i c u l t yo fk e r n e ld e n s i t ye s t i m a t i o n si np r a c t i c el i e s i ns e l e c t i n ga l la p p r o p r i a t ev a l u eo ft h es m o o t h i n gp a r a m e t e r , t h eb a n d w i d t h ，f o ru s e w i t ha g i v e nd a t as e t i nt h i sp a p e r , w ec o n s i d e rt h eb a n d w i d t hs e l e c t i o np r o b l e mf o r k e r n e l d e n s i t ye s t i m a t i o nf r o mr a n d o m l yf i g h t c e n s o r e dd a t ab a s e do n am e a n a b s o l u t ee r r o rc r i t e r i o n 1a s tb u tn o tt h el e a s t ，w h a ti st h ei n n o v a t i o no ft h i sp a p e r , w e s y n t h e s i z e dt w om e t h o d st op e r f o r mb e t t e r , e s p e c i a l l yo nb o u n d a r y k e y w o r d s ：k e r n e l e s t i m a t i o n ；r a n d o mc e n s o r e d ；c r o s sv a l i d a t i o n ； n o n p a r a m e t r i ce s t i m a t i o n v 参考义献 符号说明 ，【】表示某事件的示性函数 互，互，互表示非负独立同分布的表示寿命的随机变量 c l ，c 2 ，e 表示非负独立同分布的表示删失的随机变量 吃表示窗宽 m s e 表示均方误差 o ( n 五7 5 ) 表示，l 一纠5 的同阶变量 c v 表示平方拟合误差 芋一e x p ( 1 ) ；均值为1 的指数分布 鸭( o ，5 ) ；均值为0 ，方差为5 的正态分布 0 为零向量 i 为单位阵 f ( t ) 表示f p ) 的估计 北京化工大学位论文原创性声明 本人郑重声明：所呈交的学位论文，是本人在导师的指导下， 独立进行研究工作所取得的成果。除文中已经注明引用的内容外，本 论文不含任何其他个人或集体己经发表或撰写过的作品成果。对本文 的研究做出重要贡献的个人和集体，均已在文中以明确方式标明。本 人完全意识到本声明的法律结果由本人承担。 作者签名：堡垒垃生查：izh 1 以 日期：星qf 垒：! 呈鱼 关于论文使用授权的说明 学位论文作者完全了解北京化工大学有关保留和使用学位论文 的规定，即：研究生在校攻读学位期间论文工作的知识产权单位属北 京化工大学。学校有权保留并向国家有关部门或机构送交论文的复印 件和磁盘，允许学位论文被查阅和借阅；学校可以公布学位论文的全 部或部分内容，可以允许采用影印、缩印或其它复制手段保存、汇编 学位论文。 保密论文注释：本学位论文属于保密范围，在上年解密后适用 本授权书。非保密论文注释：本学位论文不属于保密范围，适用本授 权书。 作者签名：隍垒鱼墨童l 踟i 曲期：丞坠fq k 丞 导师签名：盔盘主登狮日期：导师签名：么盘墨亟型刈日期： 讪s 砂6 绪高 绪言 本文将着重非参数的介绍，特别是近一、二十年来这一领域的新进展，以下 我将介绍一些有关的概念，这可能有利于读者在阅读这本章之前对生存分析有一 个大致的了解。 1 生存时间 生存时间可以广泛地定义为一给定的事件的发生的时间。这个事件可以是疾 病的发生、一种处理( 治疗) 的反应、病情复发或死亡。因此，生存时间可以是 无肿瘤时间，从一种食疗开始到有反应的时间。生存数据可以反映生存时间、对 治疗的反应以及与反应、生存及疾病发生有关的病人特征。生存数据不仅出现在 生物医学中，而且出现在工业可靠性、社会科学和商业研究中。在这些领域的生 存数据的例子是：可靠性工程中电子设备( 元件或系统) 的寿命，犯罪学中重犯 人的假释时间，社会学中首次婚姻的持续时间，它也可以不是时间，它可以是汽 车车轮转动的圈数，也可以是市场学中报纸或杂志的篇幅或订费，甚至可能是保 险公司在某一索赔案中所付的保险费。 2 随机删失数据 设五，互，乙是非负独立同分布的表示寿命的随机变量，其分布函数为f ( f ) ； g ，c 2 ，e 是非负独立同分布的表示删失的随机变量，具有分布函数6 ( t ) 。在 随机右删失中，我们不能完全观测到z ，而仅能观测到 五= m i i l ( z ，c ：) ，谚= i t is c ：】，f = 1 ，2 ，n ， 其中，【】表示某事件的示性函数。显然，6 包含了删失信息。 下图表示随机右删失数据的例子： 这里，第一个病人在t = 0 时刻进入研究，但在研究停止前某时间死亡。第二 1 绪自 个病人在研究丌始某时刻进入研究，但在研究结束时，该病人仍然活着，因而产 生了一个删失观察；第三个病人在研究丌始后进入研究，但在研究截止前退出试 验，因而产生另一个删失观察。这种删失通常发生在医学实验及临床研究中。在 大部分临床研究试验中，研究时间通常是固定的，且病人通常在这段时间内的不 同时间进入研究。一些人在研究结束自“可能死亡这部分人生存时间已知；其他 人可能在研究结束之前退出试验或失去联系，或在研究结束时仍然活着。对于中 途退出或失去联系的病人生存时例至少是从进入研究到研究结束前这段时间；对 仍然话着的病人，生存时侧至少是从进入研究直到结束这段时闻，这两种观察就 是删失观察。既然进入实验的时间可能不同，因而被删失的时问也可能不同。例 如，假设6 个有急性白血病的病人，先后进入研究时间为1 年的临床试验，假设 6 个病人获得治疗并得到缓解。缓解时间见下图 控 一失女# i i 一研究停止 一复数 一研究停止 “”。2 = = 剖 斌 旨 病人a ，c ，e 获得缓解的时间分别是二月，四月，和九月，而复发的时间 分别是4 和5 及3 个月以后。病人b 在三个月开始获得缓解，但4 个月之后失 去联系，对b 的缓解时间至少是4 个月。病人d 和f 分别在第5 、第1 1 个月才 开始获得缓解。且在研究结束时还没有复发，因此他们的缓解时间至少是8 个月 和3 个月。这六个病人的缓解时间分别是4 ，4 + ，5 ，8 + ，3 ，a n d 2 + ，其中带“+ ”号的数 据表示删失数据。 3 生存时间的函数 生存时问的特征通常用下面3 个函数来刻画：( 1 ) 生存分布函数；( 2 ) 概率 密度函数：( 3 ) 失效函数。应当指出：这3 个函数在数学上是等价的。设非负随 机变量r 有密度，o ) 和分布函数f o ) ，则生存函数定义为f o ) 一1 - ，o ) ，失效 率为砸) - ，( i f ( r ) ) 。本文将探讨这些生存时间密度的具体估计程序及其 中变窗宽的选取方法。 2 4 随机删失回归 随机删失回归主要包含随机删失线性回归、半参数回归、非参数回归及c o x 回归等。先以线性回归为例说明。设有线性模型y ；x i r 卢+ 疋- - l 2 , - , n ) ，其中 卢是参数变量，是均值为零且方差有限的随机误差变量，墨是协变量向量，誓 是响应变量( i ；1 ，2 ，n ) 。在完全观测下，最小二乘法可以估计参数芦。然而在 随机删失下，z ( 通常是某寿命随机变量的对数) 可能被另一随机变量c 所删 失，使y 不能被完全观察，而观察到的数据是( z ，谚，五) ，其中 z ：= l t l i n ，g ) ，4 = ，嘭s g 】“- 1 , 2 , ，忍) 。这种响应变量被随机数据删失的回归 就成为随机删失回归。显然像标准方法，比如线性回归就不能直接应用。于是如 何利用随机删失数据对回归模型进行统计分析，正是随机删失数据回归所要探讨 的内容。 应当指出这里所介绍的随机变量的随机删失只是随机右删失，其他删失还有 左删失、双向删失及区间删失等，但在此不作一一介绍，因为本文主要集中在随 机右删失方面的内容。 3 第1 章1 r 参数密度核估计 第1 章非参数密度核估计 非参数模型的研究是当前计量经济学研究中的又一个重点方向，是继协整理 论之后，国际计量经济学的又一个热点的研究方向。在本章中，我将对核估计方 法的导出作一个简要的回顾和推导，阐述了核估计的估计原理，这将对核函数的 选取以及窗宽的选择的理解有所帮助。光滑参数的选取是后续章节的重点，并介 绍了一些选取窗宽的方法。 第1 1 节密度函数的非参数估计方法 1 1 1 核密度估计 本章讨论数据为删失情况下，密度函数非参数估计的具体方法。 假设墨，x 2 ，鼍是从服从密度函数为厂的总体中抽出的相互独立的样本。 k ( ) 为有有界支撑的对称的概率密度核函数。h 为窗宽，则厂o ) 在内点x 上的 核估计为 夕 ) = 忑1 善k 咩) 1 ) 刀疗左呷 ，l 厂 ) 表现的好坏与核函数k ( ) 和窗宽h 有关，其中窗宽的影响至关重要， 当窗宽较大时( 2 1 ) 式中起作用的与x 有关的五较多，使厂的估计较稳定， 方差波动较小，但正因为窗宽较大，考虑的f ) 较多，使得估计的偏差 f ( x ) zf ( x ) 较大。反之，当窗宽较小时，估计的偏差较小，但也因数据较少， 使得厂的波动方差较大。 综上，在密度函数厂 ) 的估计中，选取适当的窗宽是具体计算中的一个最 重要的问题。另外，由于在边界点的密度较少，样本出现的相对较稀疏，也同样 会影响到厂g ) 在边界点上厂 ) 的估计，这就引出在函数估计时的变窗宽的思想 的提出，从而使得相应的偏差较小。其中详细的步骤我将在后续的章节中给出。 1 1 2 局部线性密度估计 假设墨，五，瓦是从服从密度函数为厂的总体中抽出的相互独立的样本。 设鼍= 丁+ ( f - 1 2 ) b ，f = 1 ，2 ，g ，其中b 为根据样本容量变化的正数，并且 r ，g 和6 的选取能够保证没有数据落在区间【z ，t + 曲】之外。定义如下的示性 函数 笫1 带1 f 参数密度核估计 q 2 再h 2 删) ( x ，) ，m 玑g 2 使用下述的q 的表达式来近似厂 ) ，z 一场k _ 矿1 j 州乒+ b 2 2 厂 ) 一厂 ) 觚 ( 1 3 ) 然后使用局部线形回归的思想去拟合 哮薹( 圪。q 一蠢反( 鼍吖烀k ( 节肌f _ 堵一，咒 ( 1 4 ) 其中罗成( 而一x ) 为在x 点的，l 阶泰勒展开，h 为窗宽，k ( ) 为核函数，满 足k ( “) o ，r k ( u ) d u = 1 ，r k o , ) 以u = 0 。 设i = ( ：i 三三二! ! ；兰二) ，虫= k 呼) k ( 等- ) h 6 = c n h ，于是有计算公式 刚2 陋p 一度= 孛$ 这样，使用这种局部线性密度估计方法估计出的反代表了对应厂 ) 的七阶 导数。因此，使用此方法同时可以估计出厂 ) 其他导数的值。 当我们估计厂( ”时，就可以使用下列公式 呀n 善g 矿c i 卅2 k 呼m ，吼玑以 ( 1 6 ) 令，o ) 一卢，对上式中的卢求导取最小值，可得下式 抽= 去扣分l c f ) k ( 等) i 都，刀 ( 1 7 ) s 第1 章t 1 - 参数密度核估计 在密度函数估计的实际使用中，上式中的项n 一易k 通常是确定的数，因此上 述方法通常应用在对估计出的密度函数进行光滑，尤其是在变窗宽的生存数据的 估计问题中，虽然估计的精度较高，但是随之产生的是较高的方差，使密度函数 出现震荡，这时我们可以把此时估计出的密度函数数据当作非删失数据，使用上 述方法进行二次估计，光滑后的结果会令人满意。 第1 2 节核函数的选择 核估计的核心问题就是核权函数的选择和窗宽的选取。核权函数在核估计中起光 滑的作用，即消除扰动的随机因素，使所得曲线反映变量之间的实际经济联系。 对于就是变量是随机变量的情形，利用表1 1 ，o ) 的渐进均方误差可表示为： m s e = c r ( k ) ，l 一1 h 一1 + c 暑肛；( k ) j i l 4 ( 1 8 ) 其中g 与g 是与核函数无关的量，同时有 ：岱) = 降 弘， r ( k ) 5 严2 灿 表1 1 回归函数核估计的逐点渐近偏和方差( 内点) 渐近偏差渐进方差 h , 2 ( m f 驯+ 2 m 厂) ( 目 嘉) h f 、 7 最小化( 1 8 ) 式，即对f 式 ( 吃2 紫郴，卜蒜删 中的吃求导，使之最小可得到最佳窗宽： ”斟5 r ( k ) i 5 刀- 1 ，5 带入( 1 8 ) 式可得： m s e o e r = 厅刮5 ( o ) 4 巧( g ) 2 届( 4 1 巧+ 4 _ 4 届皿4 艏晖) 正巧( k ) 所以，最优的核函数是使 r 2 似) u ：衅) = ( 肛2 肋) 2p 2 k ( u ) d u 6 ( 1 9 ) ( 1 l o ) 第1 章1 参数密度核估计 达到最小的核函数k ) = 0 7 5 ( 1 一“2 ) + 第1 3 节窗宽的选择 窗宽是控制核估计精度的重要参数。太小的窗宽得到除了数据点外其它点的 函数值都为零的函数。所以，太小的窗宽会使得随机误差项产生的噪音没有被排 除，是没有意义的估计。而太大的窗宽得到过分光滑的曲线，接近于直线，此时 的估计也没有任何意义。核估计的实际应用中，如果回归函数的核估计接近于一 条直线，则窗宽肯定过大，参加局部加权的观察点过多，此时可减小窗宽。如果 回归函数很不光滑，则窗宽肯定过小，此时，随机误差项产生的噪音没有排除， 应该加大窗宽，使得在局部参加加权的观察点增多，从而更多地消除随机误差项 产生的噪音。 1 3 1 理论窗宽的最佳选择 对于解释变量是随机变量情形，核估计在解释变量在取值区间内点处的逐点 渐进偏差和渐进方差如表1 1 ，估计回归函数的渐进方差随着窗宽的减小而增大， 渐进偏差随着窗宽的减小而减小。所以，非参数估计就是在估计的偏和方差中寻 求平衡，使得均方误差达到最小，由均方误差的表达形式： m s e = c v r ( k ) n 一1 h 。1 + q a 2 ( k ) h 4 可直达到最小的最佳理论窗宽应具有形式 吃= 册一1 7 5 其中c 和n 无关，至于回归函数、解释变量的密度函数和核函数有关。此时，核 估计在内点处可达到最佳收敛速度为o ( n 硎5 ) 。 应用最佳理论窗宽k = c r t - 1 5 ，必须先估计c ，而对c 进行估计会产生偏差， 所以，在实际应用中最佳窗宽的选择是不断地调整c ，使得采用窗宽j l 。= 册- 1 5 的核估计达到满意的估计效果。 1 3 2 窗宽的经验选择方法 在上文中我已经提到窗宽选择过大，会使估计出的密度函数直观上看来过于 平缓，反之，过小的窗宽则让密度函数呈现锯齿状的大幅波动。因此，依据这两 个直观表现，我们在实际的密度函数的估计中，可以不断调整窗宽的大小，直观 的观察估计的效果，直到选取到满意的窗宽为止。 7 笫l 章非参数密度核估计 但是这种方法的主观性较大，究竟什么是令人满意的结果很难说清楚，密度 函数的光滑程度在实际的评价中也很难确定。即使对同一组数据，使用这种“经 验方法”作密度函数的估计，也很有可能出现不同的结果。 给予这种方法的缺陷，我们采用“初等窗宽修正”方法来弥补。实际上就是窗 宽的二次选取，这种方法的原理是：对于一组数据，它的数据点的分御是不同的， 而经验选取方法自始至终使用同一窗宽来估计密度函数，这样就会使数据点集中 的地方的窗宽取得偏大，由于参与计算的数据点较多同样使密度函数估计较大。 反之亦然。于是在窗宽的二次选取中，我们令 厶 。= ( 1 1 1 ) 1 1 ) ，z n = 百一 ( 1 。 厂 ) 这样，就可以使窗宽大的地方相应减小，同时也可以增大窗宽较大的地方。 达到对窗宽的初步调节，修正窗宽的经验选取方法中的误差，优化密度函数估计 效果。 1 3 3 样本窗宽的交错鉴定选择方法 交错鉴定方法是选择窗宽的一个常用方法，其基本思路是：在每一个局部 观察点z = 五，首先，在样本中删除观察点( 五，￥) ，其次，将剩下的n - 1 个观 察点在x = 置处进行核估计： m n , - i ( x i ) f f i 川w 打j ( x f ) y j ( 1 。1 2 ) 最后，通过比较平方拟合误差 c v ( h n ) - - - - 1 。1 一m t ( 五) ) 2 w 伐) ( 1 1 3 ) 胃 选择使得平方拟合误差达到最小的窗宽，其中w ( x ) 苫0 为某函数，可以是 某权函数。在核函数估计中，这个权函数则是在密度函数估计中的核函数。该方 法的关键是在样本中剔除观察点( 五，￥) 。如果不这样的话，由于该权函数 ) 0 在观察点x = 置处达到最大值，就会使得x = 五的重要性过分夸大而 其他观察点数据的重要性程度降低。所以采用交叉鉴定方法就避免了因没剔除观 察点( 互，￥) 而将有用的数据排除在外的情况。 虽然这种方法可以在理论上得到比较高的精度，但是无疑，它的计算量将是 巨大的，这种方法将增加一个数量级的计算量来计算每一个点的窗宽i i i 。，在后文 中，我将对此种方法的缺陷加以改进，同时这也是本文的创新点之一。 8 第1 章托参数密度核估汁 第1 4 节核估计的例子 上文介绍的非参数回归模型的核估计即适合于解释变量是随机变量的情形， 也适用于解释变量是确定性的情形，还适合非参数的时间序列模型。本小节先分 别使用均值为1 的指数分布，标准正态分布，均值为o ，方差分别为5 和o 5 的 正念分布生成随机数据，应用窗宽的经验选择法暂定h n = 0 8 以- 1 ，5 和理论上估 计效果最好的e p a n e c h n i k o v 核对随机数据进行核密度估计，做出图像观察估计 结果。然后根据估计结果，使用交叉核实法窗宽的交错鉴定法改进估计 不甚理想的密度函数。本节例子的核估计都是由软件m a t l a b 完成的，在估计 过程中，并没有实现核估计的致信区间和变窗宽估计的计算，这方面的讨论将要 在后续章节中讨论。本章节所给出的程序只要稍加修改就可应用于任何数据集合 的密度函数的估计。 1 4 1 窗宽选取对估计效果的影响 解释变量为确定变量，随机生成的数据分别满足，均值为1 的指数分布 亭一e ) ( p ( 1 ) ；均值为0 ，方差为1 的j 下念分布够一( 0 ，1 ) ；均值为0 ，方差为5 的 正态分布 魄一( o ，5 ) ；，均值为0 ，方差为0 5 的正态分布仍一( o ，o 5 ) ；应 用窗宽 。= 0 8 ，1 1 5 和e p a n e c h n i k o v 核对随机数据进行核密度估计，得到如下 效果图： 图表l 1 核密度估计效果图 9 第1 章1 r 参数惭度核估计 从下图中可以看出，对于锡- n ( o , o 5 ) ；的估计效果最好，对于 仍一n ( 0 ，5 ) ；估计的密度函数震荡得很严重，这说明窗框取得过小。另外，对 于亭e x p ( 1 ) ；，在内点估计效果较为满意但是在不为零的边界点，出现了差距 较大的情况，这也是影响核估计的边界效应问题，这种情况主要是依靠变窗宽的 方法加以改善，这方面的讨论将在后面的章节中，即生存数据密度估计的时候进 行研究。 1 4 2 交叉核实法 使用上文提到的交叉核实法改进对仍一( o ，5 ) ；的估计进行改进改进的 结果如下图所示。 图表1 2 交叉核实法估计效果图 可以看出，由于交叉核实法使在内点的平方拟合误差达到最小，使这点的估 计效果明显改善了，从图形上来看，方差和偏差都可以达到令人满意的结果。 第1 5 节多元密度函数的核估计 设p 维随机向量x 的密度函数厂( x ) = 厂瓴，) 未知， 他的一个独立分布的样本，则l ( x ) 的核估计为： 五o ) ；赤羹k f x f i - x ) 1 0 x l ，x 2 , - - - , x n 是 ( 1 1 4 ) 第l 章1 r 参数崭度梭僻汁 其中k ( ) 是满足 k ( u ) z 0 ，严( u j d u 2 1 俨( u ) u d u o , f k ( u ) u u l d u 。p z 衅) i 的多元函数其中0 为零向量，i 为单位阵。满足这些条件的最常用的核函 数为 刖= 警铲( 1 一砰) + ( 1 - - 5 ) 其中s d - h ”r ( d 2 ) 可以证明，多元密度函数的核估计在内点处一致性和渐近j 下态性。当是严格 平稳序列时，该结论仍然成立。 图表l 3 模拟数据点分布图 我选取二元标准正态分布来产生模拟数据，上图蓝色的点显示模拟数据服从 二元标准正态分布。使用上述方法来对这些数据进行密度函数的估计，在这里的 核函数我并没有选择上文中常用的多元核，而是使用多元高斯核，也就是标准正 态核 ，卜丽万1 叶南 堕一堑皋掣+ 盟茅 其中tq - 呸- 1 ；a am a 2 o ；r - 啦为二元标准正态核。 第1 非参数* 度核# 计 依题意，应该选用欧式空问距离作为不同点之间的窗宽，把两点的坐标分别 代到高斯核的表达式中，在宙宽范围之的数据给与的权重几乎为零，是比较方便 的核函数，不需要再编程序的时候给与过多的判断。 得到的估计效果如下图所示，直观来看估计效果还算理想。 图表1 4 估计的密度函数效果图 第2 章生存数据密度估计 第2 章生存数据密度估计 在生存数据统计分析中，生存数掘的概率密度函数是主要的研究对象，而密 度函数的估计是非参数统计学中的必不可少的组成部分。自从上世纪6 0 年代以 来，密度函数的估计已经被广泛地研究，同时提出了许多估计方法，例如，核估 计光滑方法和紧邻估计方法，见r o s e n b l a t t ( 1 9 5 6 ) ，l o f l g e r d e n 和q u e s e n b e r r y ( 1 9 6 5 ) 1 , 2 1 。 在本章中，我将首先对比非参数统计核估计与生存数据密度估计方法在估计生存 数据密度函数时的效果，直观说明删失估计方法在处理生存数据时的必要性。接 下来再下一小节中，我将详细介绍两种针对删失数据的统计方法，即 k a p l a n m e i e r 乘积限估计和“薛留根”估计法。然后，我将提出四种针对删失数 据的窗宽选取方法，并分别把定窗宽与初等变窗宽，紧邻估计法与交叉核实法进 行分组比较，并画出偏差图讨论它们的估计效果。结合上述讨论的结果，我提出 了一种综合上面各种方法的综合估计法，在改进上述方法缺点的同时，试计算量 得到控制。为了验证每种方法的稳定性，同时考虑到我的计算机的计算能力，本 章中的每种方法我们分别模拟了1 0 0 次，计算出的它们的渐近平方平均误差，比 较它们估计效果的优劣。最后介绍一些关于生存数据密度估计的其他的方法。 第2 1 节生存数据统计方法 2 1 1k a p l a n m e i e r 乘积限估计 设丁表示生存函数，f o ) = p ( 丁s f ) 表示丁的分布函数，则确) = 1 一f o ) 定义z 的生存分布函数，它实际上是个体生存时间长于t 的概率。易知，f ( o 是 非增函数，且，( o ) = 1 ，f ( o o ) = o 。 函数f i j f ) 也叫累积生存率，它的图形叫做生存曲线。陡峭的生存曲线表示 低的生存率；较平坦的曲线表示高的生存概率。 在实践中，如果数据被完全观察，生存函数可用生存时间长于t 者所占的比 例来估计： 郧) = 塑苇铲 ( 2 1 ) 这里f o ) 表示f o ) 的估计。当属具有删失时，是( 3 1 ) 的分子一般不能确 定。例如考虑下面的生存数据：4 ，6 ，6 + ，1 0 + ，1 5 ，2 0 ，其中带“+ ”号的数据表 示是删失数据。利用式( 2 1 ) 可得f ( 5 ) = 差= o 8 3 3 ，但不能得至u f 0 1 ) ，因为 生存时问长于1 1 的病人是不知道的，第三个病人或第四个病人的生存时问可能 长于也可能小于1 1 。因此，一旦有删失数据，用式( 2 1 ) 估计f ( f ) 是不合适的。 于是在随机删失下，构造，p ) 的估计是本节的内容。 构造f q ) 主要有两种方法。第一种是生命表分析法，这种方法适合于样本 量很大( 例如数以千计) 或数据是按区间分组等情况；第二种方法是k a p l a n 与 m e i e r 所提出的估计生存函数的乘积限方法。由于计算机应用越来越广泛，这个 方法可用于小样本，中样本和大样本各种情形。生命表估计与乘积限估计实质上 是一样的。很多作者也把乘积限估计称作寿命表估计，两者的差别是：乘积限估 计是基于一个一个的数，而寿命表估计是基于按区问的分组数据，因而乘积限估 计是寿命表估计在各个区间只含一个观察值的一种特殊情形。由于人们主要致力 于乘积限的研究，并获得一系列丰富的成果，因此在此只介绍乘积限估计方法。 设互，互，互是非负独立同分布表示寿命的随机变量，其分布函数为f ； g ，c 2 ，e 是非负独立同分布表示删失的随机变量，具有分布函数g 。在随机 右删失模型中，我们不能完全观测到i ，而仅能观察到 工，= n l i n ( 巧，g ) 睡= ，晖s c ：】，f = 】，2 ，l ， 显然五有分布函数日( f ) 一p ( xs f ) = 1 一( 1 一f ( f ) ) ( 1 一g ( f ) ) 。k a p l a n 和m e i e r 针对这一随机删失数据提出了生存分布f q ) 的乘积限估计。 我们观察到的数据对是，4 ) ，哦) ，色) ，设x ( 。) x ( ：) x ( 。) 是 五，五，瓦的次序统计量，k a p l a n 与m e i e r 所定义的乘积限估计是 赫h 薪职b n i r = 职( 11 ) 4 ， 自从r o s e n b l a t t ( 1 9 5 6 ) 及p a r z e n ( 1 9 6 2 ) 1 2 1 的工作以来，在完全样本下核 密度估计可能是最流行的估计之一，并得到深入的研究。然而对于随机右删失， 概率密度核估计的第一篇文章是b l u m 与s u s a r l a ( 1 9 8 0 ) 。基于 ( x ，4 ) ( i 一1 ，2 ，1 ) ，b l u m 与s u s a r l a ( 1 9 8 0 ) 定义f ( t ) 的如下估计： 删嵋赁k ( 管弦 ( 2 3 ) 其中吃是趋于零的常数序列，k ( ) 是核函数，f 。一1 一f 。，e 是k a p l a n 与 m e i e r 所定义的乘积限估计。 1 4 笙! 里竺堡鍪塑堕丝! 查生 一 2 1 2 概率密度函数的“薛留根”核估计 设置，噩，鼍是独立同分布的随机变量，具有分布函数，和密度函数厂。 x ，k ，k 时表示删失的独立同分布随机变量，具有连续的分布函数g 。假定诸 五独立于y 。在生存分析领域的研究中，x ，e 常表示寿命，它往往因随机 右删失而不能完全观测到，即仅能观察到数据( 五，4 ) ，f = 1 ，2 ，z 其中 z ： = m i n ( x f ，￥) ，4 = ，s ￥) ，这里，( ) 时表示某事件的示性函数。易知 z 1 ，z 2 ，z ：仍是独立同分布的随机变量，具有分布函数h = i - 0 - f ) ( 1 - g ) 。 问题是要用样本( z l ，磊) ，皈，瓯) 对厂o ) 做出估计。当g 己知时，定义下面的 核估计 胁击纠z i 吃- - xi j l - g ( 多i 葡 其中 ) 是趋于零的常数列，k ( ) 是某个b 。r e l 可测函数【3 l 。 实际中g 一般是未知的，此时六i x ) 还不能够被使用，因此，当g 未知时， 采用g 的k a p l a n - m e i e r 估计 3 1 ，即 g ) = 以 1 - l = l 1 屈o ) ，z 功， 屈 ) = ，( 谚= o ，互sz ) 。用g ：l 代替g 即可得 到l ( x ) 的估计【3 】 m 去薹k 降z i _ - - x i j l - g 6 i 两 g 。6 ) 接下来我分别使用针对生存数据的核估计方法和可应用于非删失数据的非 参数统计方法对删失模拟数据进行密度函数的估计，生存数据的核估计方法使用 b l u m 与s u s a r l a 的k a p l a n m e i e r 乘积限估计，而非删失数据的非参数方法使用 定窗宽核估计法。此次模拟中使用的数据先是由均值为1 指数函数生成，再由均 值为3 4 的指数分布截断而产生的模拟删失数据，为更好地观察和对比估计结果， 第2 章生存数据惭度估计 这两种窗宽都是使用e p a n e c h n i k o v 核作为核估计函数，定窗宽选为 玩= 0 8 n 一1 7 5 得到的估计结果如下图所示 图表2 l 删失数据核估计与非参数核估计估计效果对比图 如图示指出，红色的曲线表示删失数据的真实密度函数，绿色曲线为为删失 估计方法，蓝色曲线为删失估计方法，下面的红色圆点指示出了删失数据点。从 上图中我们可以直观地看出，绿色线在生存数据的估计方法较为精确，在函数估 计的前半部分，非删失方法得到的密度估计值较大，那是这种方法把所有的数据 当作真实数据来处理，在估计一点的密度时，把窗宽内的所有数据加权后求和， 当我们观察删失数据估计方法的密度函数估计式后发现这种方法“剔除”删失数 据后再加权求和的，所以跟接近真实密度函数。而这种方法后半部，未删失估计 法估计的密度函数的值明显小于k - p 估计法的值，这是因为删失估计法的估计 式中有修正项队估计值进行适当的修正，所以生存数据的密度函数估计效果在整 体上较为令人满意。综上，在对含有删失数据的密度函数进行估计时，乘积限这 种生存数据统计方法还是比较必要的。 另外，我们也可以看到两种方法对边界点的估计都不甚理想，尤其是当自变 量趋近于零时的密度函数估计效果两种方法都不理想。这种“边界效应”我将在下 面的章节中加以讨论和改善。 1 6 第2 章生存数据密度估计 第2 2 节窗宽选取方法 非参数密度估计是一种强大的数学分析工具，它能够用于鉴别基本分布函数 种类和形式的调查。在多种可用的非参数密度估计中，由于核估计的简单形式和 光滑的性质，使它成为密度估计中最普遍的方法。在实际中，核估计的最主要的 困难是在给定的数据集合中，选取合适的光滑参数，即窗宽的选取问题。在下文 中，我将要给出四种窗宽的选取方法，并对他们的估计效果进行比较。 2 2 1 定窗宽与初等变窗宽 定窗宽的选取方法与未删失情况下的非参数核估计经验选取方法基本类似。 通过调整j i l 。= c n - 1 5 中的c 来使得到的密度函数有较好的光滑度，但又不至于使 估计出的结果过于平坦以至于失去实际意义。c 通常在( 0 ，1 1 之问取值。但是这种 窗宽选择方法显得过于粗糙，只是依靠直观观察，虽然在大样本性质上可以达到 理想的收敛速度，但c 的选取依靠估计效果来判断，且含有很大的个人主观性在 其中，可以说并不是一个统一和可靠的方法。初等变窗宽改进了经验窗宽选择中 最佳窗宽选取的困难，这种方法采用变窗宽的思想，使用的窗宽定义如下： 办= 厂 其中，为根据经验定窗宽选出的估计效果相对较好的窗宽。在点x 的窗 宽为这点的定窗宽与密度函数的比值。这样，由于数据较为集中的区域密度 函数估计值相对较大，在这些店的估计窗宽就相应变小，反之，在数据点稀疏的 区域，甚至是边界点，就可增多参与估计这点密度值的样本，提高了估计的精度， 这种方法的估计原理在前面的非参数估计方法中已有详细介绍，在此不再叙述。 其中的模拟数据使用均值为1 的指数分布生成，用均值为3 4 的指数分布截 断尾而产生删失数据。选区的窗宽与上文中非参数估计法相同，都是 吃= 0 8 ，l - 1 5 。下面的子图指示出了各估计方法与真实密度函数的偏差值，从 中我们可以直观看出不同估计方法的估计效果。 下面是使用m a t l a b 程序编制的两种方法的对同一组删失数据的对比图： 1 7 图表2 2 定窗宽与初等变窗宽估计效果对比图 从偏差对比图中可以看出，在内点初等变窗宽的效果明显好于经验定窗宽估 计法。 2 2 2k 近邻估计法 令1 k n ，记 l 乒= f ：而是离支最近的忌个观测值之一) m i e l n i c z u k ( 1 9 8 6 ) 【2 1 定义下面包一紧邻估计： o ) 2 页南荟哌r ，忌) ，j j x 乒 ( 2 。7 ) 其中比， ，七) 是某核函数，满足 形； ，k ) =一箭卜乒 g 0 ， i 圣j i , k 其中c 是上文中所定义的k a p l a n m e i e r 乘积限估计，尺0 ) 是f 到第包格最 近的非删失的距离。1 2 1 实际上这种方法采用定窗宽的思想，却属于变窗宽的方法，它只是使用估计 第2 章生存数据密度估计 点左右两边一定数量的数据去估计这一点的密度函数值，窗宽却随着估计点