 Statistical and Application 统计学与应用, 2012, 1, 20-25 http://dx.doi.org/10.12677/sa.2012.12005 Published Online December 2012 (http://www.hanspub.org/journal/sa.html) Fitting of Twofold Mixed Model for Two-Parameter Weibull Distribution* Xia ol i n g Xu1, Ronghua Wang2, Beiqing Gu1, Shengrong Wu3 1Business Information Management School, Shanghai Institute of Foreign Trade, Shanghai 2Mathematics and Science College, Shanghai Normal University, Shanghai 3Yixing Import and Export Inspection and Quarantine Bureau, Yixing Email: xlxu@shift.edu.cn Received: Nov. 6th, 2012; revised: Nov. 14th, 2012; accepted: Nov. 27th, 2012 Abstract: Considering the data of survival time of raying mice, we can know the reason of death by dissec- tion after the mice died. One reason is thymus lymphoma, and the other is reticulum cell sarcoma. In this pa- per, we give the statistical method to separate the survival data caused by two reasons without dissection. It is also considered that this kind of mice survival time can be fitted by twofold mixed model of two-parameter Weibull distribution. Besid es, the optimum selection criterion is proposed, and the inverse moment estimates and interval estimates of parameters are obtained. The calculated results indicate that the method is simple and feasible. Keywords: Weibull Distribution; Twofold Mixed Model; Thymus Lymphoma; Reticulum Cell Sarcoma; Inverse Moment Estimate; Interval Estimate 两重两参数威布尔分布混合模型的拟合* 徐晓岭 1，王蓉华 2，顾蓓青 1，吴生荣 3 1上海对外贸易学院商务信息学院，上海 2上海师范大学数理学院，上海 3宜兴出入境检验检疫局，宜兴 Email: xlxu@shift.edu.cn 收稿日期：2012 年11 月6日；修回日期：2012 年11 月14 日；录用日期：2012 年11 月27 日 摘 要：考虑被辐射的老鼠的存活时间的数据，当老鼠死之后，可以经解剖知道其死亡原因，一类原 因为胸腺淋巴瘤，另一类原因为网状细胞肉瘤。本文给出在不经解剖的情形下如何将两种原因引起的 存活数据进行分离的统计方法，认为这类老鼠存活时间可用两重两参数威布尔分布混合模型来拟合， 同时提出了一个优选准则，并得到了参数的逆矩估计和区间估计。计算结果表明本文方法的简便易行。 关键词：威布尔分布；两重混合模型；胸腺淋巴瘤；网状细胞肉瘤；逆矩估计；区间估计 1. 引言 文献[1]给出了一个有关于 60只被辐射的老鼠的 存活时间的数据。当老鼠死之后，经解剖可以知道 其死亡原因。其考虑两类死亡原因所引起的死亡数 据。一类原因为胸腺淋巴瘤 (包含 22 个数据)，另一类 原因为网状细胞肉瘤(含有38 个数据)，数据见表 1、 表2。 *资助信息：本项目得到国家自然科学基金(11141002)；上海对外贸 易学院 2011~2012 年度“085”重点学科建设项目(Z08511009)；上 海对外贸易学院 2011~2012 年度“085”重点专业建设——《全球 通用商科人才培养本科专业建设》项目；上海市星系与宇宙学半解 析研究重点实验室(SKLA1101)；上海高校选拔培养优秀青年教师 科研专项基金(swm10010)资助。 Copyright © 2012 Hanspub 20  两重两参数威布尔分布混合模型的拟合 Table 1. Thymus lymphoma (22 data) 表1. 胸腺淋巴瘤(22 个数据) 159 189 191 198 200 207 220 235 245 250 256 261 265 266 280 343 356 383 403 414 428 432 Table 2. Reticulum cell sarcoma (38 data) 表2. 网状细胞肉瘤(38 个数据) 317 318 399 495 525 536 549 552 554557 558 571 586 594 596 605 612 621 628631 636 643 647 648 649 661 663 666 670695 697 700 705 712 713 738 748 753 上述数据是通过解剖得到的。如果要做很多次这 样的试验，那么解剖老鼠的工作将耗费较大的人力、 物力和财力。那么有没有不经过解剖，就能确定这些 存活时间数据中哪些是由胸腺淋巴瘤引起的，哪些是 由网状细胞肉瘤引起的呢？这正是本文所要研究的 问题。考虑受辐射老鼠的存活时间可用两重两参数威 布尔分布混合模型来拟合。关于混合模型的表述如 下： 当一个总体中每一个个体都是k个不同类个体中 的某一个，第 i类个体在总体中所占比例为 ，诸 i pi p 01 1 1 k i i p i Rx 1 k ii i pRx i Rt 1,2,, 满足和 ，这时就产生了混合模型。 i p 假定第 i类的个体服从生存函数为 的分布，从 这总体中任选一个个体，其生存函数为： Rx , 这种模型被称为“混合模型”。在不能区分不同类的 个体时，常常用此模型。诸 常取自同一参数族， 当然这非必需的。 令 1, ii F xR xi k ，即第 i类的个体 的分布函数，则混合模型的分布函数 F x 1 k ii i 为： 11 11 kk iii i ii F xpRxpR xpFx 2k 在此，称其为 k重混合模型。如果 ，则称为两重 混合模型。 混合模型是一个很灵活而且是强有力的统计建 模工具。在实际中混合模型已成为分析复杂现象的一 个重要工具，它几乎涵盖了各个学科，其广泛应用于 证券分析、声学研究、生物学、医学、经济、金融、环 境科学、植物学、农业和可靠性等领域。有关混合模 型的研究历史，第一篇关于混合模型的文章是 Pearson (1894)用两重正态混合分布拟合 Weldon 的数据， Pearson 使用了矩方法，上世纪七十年代前对混合模型 的研究大都是使用矩方法及图形方法。从七十年代开 始，对混合模型的研究更多地使用了极大似然方法， 虽然第一篇使用极大似然方法的是 Rao(1948)用 Fisher 得分方法研究两重等方差正态混合模型的参数 估计。随着高速计算机的发展以及迅速普及，混合模 型的研究大量运用极大似然方法，尤其是EM 算法的 出现真正促使极大似然方法更广泛地应用于混合模 型。Dempster(1977)为EM 算法建立了一个一般的理 论，从而解决了应用极大似然方法在混合模型计算上 的障碍。混合模型的前期理论研究成果可从文献[2-6] 中找到。 混合模型研究已经成为统计中的一个热点。主要 有两方面的原因：一是由于从数学上看，混合模型有 许多新的性质，从而带来了新的理论问题；二是在众 多应用问题中，混合模型提供了一种方便灵活的统计 建模方法，成为应用研究者建模时一种方便的选择， 具体地说混合模型比单一模型能更好地拟合现实数 据，再者，正如Pearson(1894)所遇到的问题那样，现 实大量问题本身就是一个有明确物理意义的混合模 型，这样混合模型中各参数都有明确的物理意义，这 时混合模型便是最自然的选择。在实际中，混合模型 有其物理背景。譬如，某些产品有缺陷，那就容易引起 早期失效，而没有缺陷的产品是逐渐地走向失效。有 关混合模型的理论和实际应用已有大量文献做了研 究。 指数分布的混合模型的研究见 Mendenhall 和 Hader(1958)即文献[7]，Tallis 和Light(1968)即文献[8]。 威布尔分布的混合模型的研究见 Kao(1959)即文献 [9]、Falls(1970)即文献[10]。伽玛分布的混合模型的研 究见 Dickinson(1974)即文献[11]，蒋仁言和左明健在 文献[12,13]较为系统地总结了威布尔混合模型的应用 实例，其应用包括：电子管、大气数据、电器设备、老 化数据、互补型金属氧化物半导体试验、作物产量估 计、仪器仪表、机械零件、品面型二氧化硅击穿数据、 筛选等等。Cohen 在文献[14]所研究过的混合正态分 Copyright © 2012 Hanspub 21  两重两参数威布尔分布混合模型的拟合 布可应用于各种风速和大量生产的器件的物理尺寸 的研究。Kao 在文献[9]所研究过的混合威布尔分布在 可靠性研究中是有用的，特别在那些涉及电子管的可 靠性问题的研究。金国騿、胡文忠在文献[15]中研究 了风速频率分布的混合模型——指数分布和威布尔分 布的混合，即在零风速和低风速时指数分布起主要作 用；在中风速段和高风速段威布尔分布应起主要作 用。文献[16]研究了混合指数分布的参数估计问题。文 献[17]给出了一个双胞胎数据的有限混合模型。文献 [18]给出了有限泊松混合的稳健统计推断。 2. 两重两参数威布尔分布混合模型的 参数的逆矩估计 设12 n ,,, X XX n 是来自总体 X的样本容量为 n的 一个简单随机样本， 12 X XX 为其次序统 计量。总体 X的分布函数为 F x，密度函数为 f x 12 1 ， 其中 F xpFxpFx , 12 1 f xpfxpfx01p, 而 i F x，fx 分别为两参数威布尔分布的分 布函数及密度函数，即 ,1,2 ii 1e i Fxxp i m i x 1exp i i m i im i m i mx x fx 1, 2i 0 i0 i i , 其中 m称为形状参数， 1 称为刻度参数。在 此假定 2 , n 。 上述即为两重两参数威布尔混合模型。所谓混合 威布尔分布表示一个总体由几个威布尔子总体组成。 早在 30 多年前，人们就认识到以混 合威布尔分布作为 具有一个以上失效起因的单元的寿命模型是合适的。 蒋仁言在文献[12]、蒋仁言和左明健在文献[13]中采用 WPP 图较为详细地研究了二重、三重两参数威布尔混 合模型的点估计。 在两重两参数威布尔混合模型中，可以认为 12 ,, X XX中有 np 个是来自总体分布函数 1 F x 的一个样本，有 nnp个是来自总体分布函数 2 F x 的一个样本。于是只要能估计出p，其它参数的估计 就很容易得到。 首先，假定 p的估计，取 ˆ,2,3,,2 j pj n n 12 作 为参数 p的估计，于是， X j XX为来自 总体分布函数为 1 F x 12 的容量为 j的前 j个次序统计 量，而 X j jn XX 为来自总体分布函数 为 2xj的容量为 n F 的前 n次序统计量。为保 证估计尽量地精确，采用王炳兴在文献[19]提出的逆 矩估计的思想可得参数 1 m j 个 1 ,；2 , 2 m 的 1 ˆ , 逆矩估计，其 依赖于 j，记逆矩估计为 1 ˆ； ˆ , 22 ˆ j j m j j m 。 1 ˆ 其中， j m为如下方程的根(方程有唯一正根)： 1 1 2, 3 ln 1, 2 jj ii Sjj j S 1 注：文献[19]指出当 j较小时，求解 m的估计可用上 述方程。事实上，文献[20]通过大量的 Monte-Carlo 模 拟说明当j较大时也是合适的。 1 1 1 ˆ ˆ 11 1 ˆ j j jm m ji i X j 2 ˆ m为如下方程的根： j 1 1 2, 3 ln 1, 2 nn ij i nj jn S jn S 2 2 1 ˆ ˆ 21 1 ˆ j j nm m ji ij X nj 11 1 ,1,2,, imm ili l SXjiXi j 22 1 ,1,2,, imm ji jl ji l SXnjiXi nj 其中 12 mm 特别地，若假定 ，并同设为 m，则参数 m的 逆矩估计 为如下方程的根： ˆ m 11 11 lnln 3 jn jn iij ii SSn SS 1 ,1,2,, imm ili l SXjiXi j 1 ,1,2,, imm ji jl ji l SXnjiXi nj 其中 Copyright © 2012 Hanspub 22  两重两参数威布尔分布混合模型的拟合 3. 两重两参数威布尔分布混合模型的 参数估计的优选准则 余下的问题是如何选择j，j，下 面 给出一个优选准则，称为极大似然准则。 2,3, ,2n ˆ 优选准则(极大似然准则)：设似然函数为 L，其 为参数的函数，而 L的依赖于 j的估计记为 j L 022 max ，选择 这样的 ，使得： 0 jˆˆ j j jn LL，或 0 ln 22 ˆˆ max ln j j jn LL n ，其中 1 !1 i Ln pfx 12 ii pfx 1 2 ˆ ˆ1 1 ˆ 2 xp ˆ xp ˆ j j m m ii j m ii j x x 1 1 2 2 1 ˆ 11 ˆ1 2 ˆ 2 ˆ ˆ!e ˆ ˆ 1e ˆ j j j j nj jm ij m j m j mx j Ln n mx j n 112 2 ,,,,pm m 注：通过大量模拟发现，在使用极大似然准则时，为 保证估计的精确性，在n很大时，通常 j的取值不应 太小或太大。另外，应根据参数的估计值作适当甄别， 例如形状参数过大等等。 于是，可得参数 11 ˆ ,,pm 的点估计 2 2 ˆˆ ˆ ,,m ˆ ： 0 ˆj p 1 ˆˆ n，0 1 j mm， ˆˆ 10 1 j 2 ˆˆ ，0 2 j mm ，ˆˆ 0 22 j 若给定置信水平1 ，则可得参数 的置信 区间分别为： ，，其中 分别为如下方程的根： 12 ,mm 00 2122 ˆ , jj mm 00 111 ˆ , j mm 2j ˆˆ 11 ˆ, j mm 00 12 ˆj 0 1 12 ln 2 1 2 ii j S 01 01 j S2 j 0120 2 12 1 jj Sj 00 2122 ˆˆ , jj mm 0 1 ln 2 ii S 分别为如下方程的根： 120 12 1 nn Snj 01 12 ln 2 ij i S 120 1 ln2 1 nn Snj 12 ˆ mmm 00 12 ˆˆ , jj mm 012 2 ij i S 若有，可得参数m的置信区间分别 为： 00 12 ˆˆ , jj ，其中 mm分别为如下方程 的根： 00 0 112 1 11 2 1 lnln2 2 2 jn jn iij ii SSn SS 00 0 112 11 2 1 lnln2 2 2 jn jn iij ii SSn SS 4. 老鼠受辐射的生存时间的实例分析 首先检验表 1和表 2的两组数据都服从两参数威 布尔分布。检验方法采用：范蒙特福特(Van-Montfort) 检验统计量。具体步骤表述如下： 设生存时间 X的分布函数为 F x12 ,,, n ， X XX 12 n 为来自总体X的容量为n的一个简单随机样本，而 X XX 0 为来自总体 X的容量为 n的前 n 个次序统计量。要检验假设 H ： 1exp m x Fx 0m0 其中 为形状参数， 为刻度参数。 令1 1 lnln ,1,2,,1 ii iii XX in EZ EZ H ,1,2,, i EZ in ，其中 为标准极小值总体容量为 n的次序统 计量的数学期望，其值可查阅文献[21]。 构造范蒙特福特(Van-Montfort)检验统计量： 1 1 1 1 n i ir r i i H r Wnr H ，其中 2 n r 0 。 21,2 H 成立时，W渐近地服从 F 当nr r 。 于是给定置信水平1 ，当 12 2 21,2 21,2 F nr rWFnr r 0 时， 接受 H ，即认为 12 ,,, n X XX为 22n 来自两参数威布尔分 布总体。 1) 针对第一组胸腺淋巴瘤的22 个数据： , 11 2 n r , 1 1 1 1.0747 1 n i ir r i i H r Wnr H 10.95, Copyright © 2012 Hanspub 23  两重两参数威布尔分布混合模型的拟合 1 20.975 21,2FnrrF 20,220.4109 2 0.025 21,220,FnrrF 22 2.3890 由于 21,2 12 2 21,2 F nrrWF nr r 0 所以接受 H ，即认为第一组胸腺淋巴瘤的 22个数据 来自两参数威布尔分布总体。 2) 针对第二组网状细胞肉瘤的 38个数据： 38n, 19 2 n r , 1 11.0372 n i i H H 10.95 1 1 ir r i r Wnr , 1 20.975 21,,380.5181Fnr 2 36rF 2 0.025 21,236,381.9190FnrrF 由于 12 2 21,2 21,2 F nr rWFnr r 所以接受 0 H ，即认为第二组网状细胞肉瘤的38 个数 据来自两参数威布尔分布总体。 因此，由胸腺淋巴瘤引起的存活数据和由网状细 胞肉瘤引起的存活时间都可认为服从两参数威布尔 分布。那么受辐射老鼠的存活时间可用两重两参数威 布尔分布混合模型来拟合。 由于胸腺淋巴瘤引起的生存时间共22 个数据， 网状细胞肉瘤引起的生存时间共有38 个数据，所以 参数 p的真值为 22 0.3667 22 p 11 , 38。 若采用胸腺淋巴瘤这 22 个数据可得参数 m 的 逆矩估计为： ， 181 ˆ2 ˆ3.474m311.018 m ；若采用网 状细胞肉瘤这38 个数据可得参数22 , 的逆矩估计 为： ， 2 ˆ7.4442m2 ˆ645.8610 。 若给定置信水平10.95 ，此时， 0.975 42 25 22 0.025 .99866,42 61.77676 99.67835 , 22 0.975 0.025 74 52.10283,74 , 则可得参数 的置信区间为： 1 m 2.8557,5.5116 2 ；参数 的置信区间为： m 6.4330,10.9231 。 由文献[1]中得到的参数估计如下(称为 E-J方法)： ˆ0.5p ， m1 ˆ5.0 ，ˆ1301.2 ； m，ˆ 2 ˆ11.62672.7 若给定置信水平10.95，此时， 22 0.975 0.025 58 38.84351,58 80.93559 1 m , 则可得参数的置信区间为： 2.2129,3.8913 2 ，从中 可以看到区间估计没有包含其点估计，在此可以认为 E-J 方法得到的参数估计精度较差；参数 m的置信区 间为： 9.2455,16.2217 24j 。 在形状参数不等的假定下，利用本文方法可解 得：当 时， 最大，此时，， 24 ˆ L1 ˆ3.7305m 1 ˆ311.2029 2 ˆ657.1054 2 ˆ9.5753m； ； ， 24 ˆ0.4 60 p 0.95 。 若给定置信水平1，此时， 22 0.975 0.025 46 29.16005,46 66.61653 , 22 0.975 0.025 70 48.75756,70 95.02318 1 m , 则可得参数 的置信区间为： 2.7085,5.1238 2 m ；参数 的置信区间为： 7.3338,12.4593 ˆ0.31p 。 文献[12]假定形状参数相等，并认为是恰当的， 用WPP 图方法得到的参数 p的估计更接近于真值， 参数估计如下： ，m， ˆ7.61 ˆ250 2 ˆ650 ， 。 若给定置信水平10.95，此时， 22 0.975 0.025 11688.0848,116 147.699 , 则可得参数 m的置信区间为： 5.0106,7.5743 17j ，从中 可以看到区间估计没有包含其点估计，在此可以认为 假定形状参数相等不一定是不合适的。 在形状参数相等的假定下，利用本文方法可解 得：当 时， L最大，此时， ， 17 ˆˆ6.1428m 1 ˆ257.3367 2 ˆ628.8087 ， ，17 ˆ0.2833 60 p 0.95 。 若给定置信水平1，此时， 22 0.975 0.025 11688.0848,116 147.699 ， 则可得参数 m的置信区间为： 4.9630,7.5913 法得到的p的估计值为 0.4，与真值的误差为0.0333 ； ，从中 可以看到区间估计包含其点估计。 评述：1) 若在形状参数不等的假定下，由本文方 Copyright © 2012 Hanspub 24  两重两参数威布尔分布混合模型的拟合 Copyright © 2012 Hanspub 25 值的误 3.8913 , 30,10.9231 , ,16.2217 ; 用本文方法得到的 085,5.1238 , .5116 , .1238 338,12.4593 , 10.9231 , .4593 由此可以看出本文的方法更合适。 WPP 图方法得到 的p 法得到的 参考文献 (References) 而采用 E-J 方法得到的 p的估计值为0.5，与真 差为 0.1333，从中可以看到本文方法是比较精确的。 另外，用E-J 方法得到的 15 2.212m [1] R. C. Elandt-Johnson, N. L. Johson. Survival model and data analysis. New York: John Wiley & Sons, 1980: 193. [2] G. J. MaLachlan, D. Peel. Finite mixture models. New York: Wiley, 2000. 2 ˆ9, ˆ11.6 6.43 7.44429.2455 m [3] B. G. Lindsay, Mixture model: Theory, geometry and application. Hayward: Institute for mathematical Statistics, 1995. [4] G. J. MaLachlan, K. E. Basford, Mixture models, inference and application to clustering. New York and Basel: Marcel Dekker, 1988. [5] D. M. Titterington. Statistical analysis of finite mixture distribu- tion. New York: John Wiley & Sons Ltd., 1985. [6] K. C. Lanf. Introduction to the special issue on finite mixture models. Social Research Methods, 2001, 29(3): 275-281. 1 ˆ3.73m05 2.7 3.7305 2.8557,5 3.4748 2.7085,5 [7] W. Mendenhall, R. J. Hader. Estimation of parameters of mixed exponential distributed failure times from censored life test data. Biometrical, 1958, 45(3-4): 504-520. [8] G. M. Tallis, R. Light. The use of fractional moments for esti- mating the parameters of a mixed exponential distribution. Tech- nometrics, 1968, 10(1): 161-175. 2 ˆ9.5753 7.3 9.5753 6.4330, 7.4442 7.3338,12 m [9] J. H. K. Kao. A graphical estimation of mixed Weibull parame- ters in life testing of electron tubes. Technometrics, 1959, 1(4): 389-407. [10] L. W. Falls. Estimation of parameters in compound Weibull dis- tributions. Technometrics, 1970, 12(2): 399-407. [11] J. P. Dickinson. On the resolution of a mixture of observations from two gamma distributions by the method of maximum like- lihood. Metrika, 1974, 21(1): 133-141. 2) 若假定形状参数相等，采用 [12] 蒋仁言. 威布尔模型族 特性 、参 数估计 和应 用[M]. 北京: 科 学出版社, 1998: 62-65, 122-125, 172-174, 199-202. 的估计值为0.31，与真值的误差为–0.0567，用 本 文方法得到的 p的估计值为 0.2833，与真值的误差为 –0.0834；而用 WPP 图方法得到的 ˆ7.6 5.0106,7.5743m ，用本文方 ˆ6.1428m4 [13] 蒋仁言, 左明健. 可靠性模型与应用[M]. 北京: 机械工业出 版社, 1999: 20, 124, 126-127, 169-170. [14] A. C. Cohen. Estimation of mixtures of poisson and mixtures of exponential distributions. NASA Technical Memorandum NASA TMX 53235, 1965. [15] 金国騿, 胡文忠. 风速频率分布混合模型的研究[J]. 太阳能 学报, 1994, 15(4): 353-356. .9630,7.59 形状 参数相等不一定是不合适的。 结论 本文利用被辐射老鼠的存活时间的已有数据，运 用两 13 ，在此可以认为假定 5. [16] 朱利平, 卢一强, 茆诗松. 混合指数分布的参数估计[J]. 应 用概率统计, 2006, 22(2): 137-150. [17] G. N. Michael. A finite mixture distribution model for data col- lected from twins. Twin Research, 2003, 6(3): 235-239. [18] D. Kalis, E. Xekalaki. Robust inference for finite Poisson mix- tures. Journal of Statistical Planning and Inference, 2001, 93(1- 2): 93- 115. [19] 王炳兴. Weibull分布的统计推断[J]. 应用概率统计, 1992, 8(4): 357-364. 重两参数威布尔分布混合模型来拟合该数据，并 提出极大似然优选准则对参数进行逆矩估计和区间 估计，得到的估计结果比原有方法更优越一些。另外， 本文的方法具有一般性，还可应用于电器装备到失效 时所使用的周期数[22]，某型号继电器进行寿命试验的 寿命检验数据[23]等其他实际应用问题。 [20] 王蓉华, 徐晓岭. 三参数Weibull 分布的统计推断[J]. 数理统 计与应用概率, 1997, 12(1): 86-100. [21] 中国电子技术标准化研究所. 可靠性试验用表(增订本)[M]. 北京: 国防工业出版社, 1987. [22] J. F. Lawless. Statistical models and methods for lifetime data. New York: John Wiley & Sons, Inc., 1982. [23] D. B. Kececioglu. Reliability Engineering Handbook (Volume 1). Endlewood Cliffs: Prentice Hall Inc., 1996. |