龙沁怡, 徐丽平*, 龙 兵
(1.长江大学信息与数学学院, 湖北 荆州 434023; 2.荆楚理工学院数理学院, 湖北 荆门 448000)
Pareto分布作为一种统计模型,在工资收入、保险精算、城市规划以及可靠性工程等领域有着极为广泛的应用. 因此,对 Pareto产品进行可靠性试验具有实际意义. 在本文中假定Pareto分布有如下的累积分布函数和概率密度函数:
(1)
(2)
其中,α(>0),θ(>0)分别被称为尺度参数和形状参数.为书写简便,将F(x;α,θ)和f(x;α,θ)分别简记为F(x)和f(x).
对于截尾样本来说,根据观测数据预测未失效产品的信息是一个至关重要的问题.在试验的早期阶段可以告诉我们测试的成本有多高,以及是否需要重新设计试验方案.目前关于预测方面的研究成果较多[13-17].本文假设受试产品的失效时刻服从Pareto分布,在双定时混合截尾样本下对产品的失效时刻进行了预测.
在这里将对文献[12]中提出的双定时混合截尾方案进行再次描述,该试验方案如下.
假设在0时刻把n个独立同分布的样品投入试验,它们的失效时刻分别为X1,X2,…,Xn,正整数m 记 即在时刻t1之前的失效样品数为m1个,而在时刻t2之前的失效样品数为m2个,m1,m2为随机变量.如果m1≥m,则在时刻t1停止试验,没有失效的n-m1个样品退出试验,其中0 记 基于上述试验数据,省略正则化常数后似然函数为 (3) 其中,I(·)为示性函数. 根据(3)式,利用极大似然法可得参数θ和α的极大似然估计分别为 (4) (5) [1-F(y)]n-k-s[1-F(t)]-(n-k), (6) 其中,y≥t. 对于Pareto分布(1)式,可以把(6)式变形为 (7) 在似然预测方法中,将运用极大似然法得到未知参数的估计和未来失效时刻的预测.基于观测数据可以得到Y和(α,θ)的预测似然函数为 根据(3)和(7)式,则Y和(α,θ)的预测似然函数为 (8) 显然α的极大似然估计为 (9) 取α=x1:n,去掉常数项后,对数预测似然函数为 nθlnx1:n-θ(s-1)lnt-[θ(n-k)+ (10) 根据(10)式,分别求对数预测似然函数关于θ和y的偏导数,则θ,Y的预测似然方程为 (11) (12) 由方程(11)和(12)可以得到Y=X(s+k):n的极大似然预测为 下面求出Y的条件分布的中位数,以此作为对Y的预测,并被称为条件中位数预测. 根据(7)式,利用二项式展开 因此Y=X(s+k):n的条件概率密度函数可以变形为 (13) dj=n-k-s+j+1. 即 由条件概率密度函数(13),可以得到 (14) 由条件概率密度函数(13)不难得到Y=X(s+k):n的100(1-γ)%等尾预测区间为(L1,U1),其中L1和U1分别满足下面的2个方程, (15) (16) 在Bayes统计分析中,未知参数先验分布的选取是非常重要的,很多文献中选择Gamma分布作为先验分布.在这里也选取θ的先验分布为Gamma分布,其概率密度函数为 (17) 这里,超参数a>0,b>0,Γ(·)表示Gamma函数. 由(3)和(17)式,利用Bayes公式可得θ的后验概率密度函数为 (18) 由(13)和(18)式可得Y=X(s+k):n的Bayes预测密度函数为 b]-(k+a+1). (19) 根据(19)式可以得到Bayes预测生存函数为 (20) 因此Y=X(s+k):n的100(1-γ)%的Bayes等尾预测区间的下限L2和上限U2分别满足 (21) (22) 可以用数值方法求出方程(21)和(22)的解,从而得到Bayes预测区间(L2,U2).事实上,利用Bayes预测密度函数(19)式也可以求出Y=X(s+k):n的Bayes预测值. 采用Bayes方法可以借助先验信息来提高统计推断的精度,如果没有先验信息可以利用,也可以取无信息先验分布,在这里取θ的无信息先验分布为 (23) 由(3)和(23)式,可得θ的后验密度函数为 (24) 实际上在(18)式中取a=0,b=0就得到后验概率密度函数(24)式,因此Y=X(s+k):n的Bayes预测概率密度函数为 (25) 由(25)式可得到Bayes预测生存函数为 因此Y=X(s+k):n的100(1-γ)%的Bayes等尾预测区间为(L3,U3),其中下限L3和上限U3分别满足 根据Z的概率密度函数(2)式,可得生存函数为 当α已知时,θ的极大似然估计为 根据后验概率密度函数(18)式,可以得到第n+1个失效数据Z=Xn+1的后验预测密度函数为 (k+a)(A+b)k+az-1× (lnz-lnα+A+b)-(k+a+1),z≥α. (26) 根据(26)式得到的Bayes预测生存函数为 (27) 由(27)式可得Z=Xn+1的Bayes中位数预测值为 zU=exp[(A+b)(γ-1/(k+a)-1)+lnα]. 下面将对文献[18]中的数据集进行分析,把这些数据按照从小到大的顺序排列如下:0.500 9、0.504 0、0.514 2、0.522 1、0.526 1、0.541 8、0.547 3、0.583 4、0.609 1、0.625 2、0.640 4、0.649 8、0.675 0、0.703 1、0.709 9、0.716 8、0.791 8、0.846 5、0.903 5、1.114 3. 借助上面的完全数据可以得到不同的双定时混合截尾数据.已知α=0.5,利用文中的结论能够计算出被截尾样品未来次序失效时刻的预测值和预测区间(γ=0.05),相关结果列于表1中. 表1 Y=X(s+k):n的预测值和预测区间(a=0.5,b=0.2) 根据不同的双定时混合截尾数据也可以得到Z=Xn+1的预测值和预测区间(γ=0.05),相关结果列于表2中. 表2 Z=Xn+1的预测值和预测区间 从表1可以看到,Y=X(s+k):n的3种预测区间的下限比较接近,当t1,t2,m和s固定时Y=X(s+k):n极大似然预测小于条件中位数预测,经典预测区间的长度最小,两种Bayes预测区间的长度比较接近,真正的观测值都落在3个预测区间的内部.在表2中对于超参数a,b的2种不同的取值及不同的双定时混合截尾试验数据,Xn+1的预测值和预测区间都比较接近,当t1,t2和m固定时,中位数预测值小于Bayes中位数预测值.从这个例子中也可以看到利用文中的结论所得到的失效时刻的预测值和真实的失效时刻还是比较接近的. 本文根据观测到的双定时混合截尾试验数据,讨论了未来失效数据的预测问题.当试验样品的失效时刻服从Pareto分布时,利用经典方法得到了被截尾样品未来失效时刻的预测值.在取两种不同的先验分布时计算了未来失效时刻的预测区间.对于独立同分布的任一产品,对它的失效时刻进行了预测.通过数值例子的计算,所得到的结果也是符合实际情况的.利用本文中的方法可以随时根据观测到的失效数据对未来的失效时刻进行预测,以便于估算试验的费用,进而考虑是否需要重新确定试验的终止时刻,在控制费用的情况下尽可能提高统计推断的精度.另外在双定时混合截尾试验数据下,也可以讨论其它寿命分布模型产品失效时刻的预测问题.2 极大似然估计及条件分布
3 未来观测值的经典预测
3.1 似然预测方法
3.2 条件预测方法
3.3 经典预测区间
4 未来观测值的Bayes预测
4.1 Gamma先验分布
4.2 无信息先验分布
5 任一观测值的预测
6 数据分析
7 结论