寿命数据的统计分析在医学、工程学、生物学等许多领域中有着广泛应用. 上世纪五十年代, 为了提高工业产品的寿命, 很多统计学家开始分析各种类型的寿命数据. 其中, 通过截尾获得试验数据是一种经常采用的方法, 在生存分析、可靠性分析等领域中经常遇到这类数据. 传统的定时截尾、定数截尾以及逐步增加的定数截尾这三类截尾数据, 很多统计学家己经广泛研究过, 有比较成熟的理论和方法, 目前的研究成果较多, 如文献[1-8]. 2009年, Wu和Kus^[9]提出了逐步增加首失效截尾寿命试验方案, 后来受到很多统计学家的重视, 目前已经产生了很多基于这类截尾数据的研究成果, 如文献[9-13]. 随着实际的需要, 结合定时截尾和定数截尾各自的优点产生了混合Ⅰ型截尾试验和混合Ⅱ型截尾试验. 如果在0时刻把$n$个样品投入到寿命试验中, $X_{r:n}$表示第$r$个元件的失效时刻, $T$表示指定的截尾时刻, 混合Ⅰ型截尾试验选择的终止时刻为$T_1^\ast=\min(X_{r:n}, T)$, 混合Ⅱ型截尾试验选择的终止时刻为$T_2^\ast=\max(X_{r:n}, T)$, 基于这两类截尾试验数据进行统计推断的相关研究成果如文献[14-17]. 无论是混合Ⅰ型截尾还是混合Ⅱ型截尾都是一个定时截尾试验方案和一个定数截尾试验方案的混合. 实际上, 在定数截尾试验中, 如果在达到预先确定的失效样品数时, 进行的寿命试验没有达到规定的时间, 为了更多地了解产品的性能和提高估计的精度, 可以考虑继续进行试验重新确定样品的失效数. 因此在本文中首次提出了一种新的寿命试验方案, 即两个定数截尾试验方案的混合, 在这里被称为双定数混合截尾试验, 具体的方案将在第2节中介绍.

两参数Pareto分布被提出来后主要用来分析社会经济和自然现象, 后来该分布被应用到生存分析和可靠性理论中. 该分布的分布函数和概率密度函数分别为

若$X$表示产品的寿命, 则可靠度函数和失效率函数分别为

其中(1.1)–(1.2)式中的$\alpha(>0)$, $\theta(>0)$分别被称为尺度参数和形状参数.

本文假设受试元件的寿命服从两参数Pareto分布(1.1), 基于双定数混合截尾试验数据对未知参数、可靠度函数以及失效率函数进行统计分析.

在可靠性试验中双定数混合截尾就是两个定数截尾方案的混合, 其模型描述如下：假设把$n$个独立同分布的元件投入试验, 事先确定时刻$t_0$以及正整数$m_1, m_2$, 且满足$m_1<m_2\leq n$. 记第$m_1$个元件的失效时刻为$X_{m_1:n}$, 若$X_{m_1:n}\geq t_0$, 则试验在$X_{m_1:n}$时刻停止, 没有失效的$n-m_1$个元件退出试验, 其中$0<X_{1:n}\leq X_{2:n}\leq\cdots\leq X_{m_1:n}$为次序失效时刻; 如果$X_{m_1:n}<t_0$, 则试验继续进行到有$m_2$个元件失效后停止, 没有失效的$n-m_2$个元件退出试验, 其中$0<X_{1:n}\leq X_{2:n}\leq\cdots\leq X_{m_2:n}$为次序失效时刻. 从上面的描述可以看出, 该模型实际上就是在两个定数截尾方案之间的一种随机选择.

在本文中把上述两种情形分别表示为Case Ⅰ和Case Ⅱ, 并且称其为双定数混合截尾试验, 从而得到如下的观测数据

根据第2节中的寿命试验模型, 下面用极大似然法求两参数Pareto分布中未知参数$\alpha$和$\theta$的估计.

设

并设双定数混合截尾试验所得观测数据为$x^*=\bigl(x_{1:n}, x_{2:n}, \cdots, x_{k:n}\bigr)$.

基于上述观测数据, 省略与未知参数无关的常量, 则似然函数可以表示为

可以求得参数$\theta$和$\alpha$的极大似然估计分别为

当$\alpha$已知时, 根据似然函数(3.1)式可以得到

则参数$\theta$的极大似然估计为

根据极大似然估计的不变性, 可靠度函数$R(x)$和失效率函数$H(x)$的极大似然估计分别为

令$Y=\ln X-\ln \alpha$, 不难得到$Y$的概率密度函数为

从(3.2)式可以看出$Y$服从均值为$1/\theta$的指数分布, 即$Y\sim\exp(\theta)$.

引理3.1^[18]  假设总体$Y\sim\exp(\theta)$, $Y_{1:n}, Y_{2:n}, \cdots, Y_{r:n}$为来自总体$Y$的前$r$个次序统计量, 令$T=\sum\limits_{i=1}^rY_{i:n}+(n-r)Y_{r:n}$, 则$2\theta T\sim\chi^2(2r)$.

对于双定数混合截尾样本$X_{1:n}, X_{2:n}, \cdots, X_{k:n}$, 令$Y_{i:n}=\ln X_{i:n}-\ln\alpha(i=1, 2, \cdots, k)$, 则$Y_{1:n}, Y_{2:n}, \cdots, Y_{k:n}$就是来自指数分布总体$\exp(\theta)$的前$k$个次序统计量, 根据引理3.1可得

即

因此, 当$\alpha$已知时参数$\theta$的置信水平为100(1-$\gamma$)%的置信区间为$(\theta_L, \theta_U)$, 其中

这里$\chi_{\gamma/2}^2(2k)$, $\chi_{1-\gamma/2}^2(2k)$分别表示自由度为$2k$的卡方分布下$\gamma/2$分位数和下$1-\gamma/2$分位数.

由于可靠度函数和失效率函数关于$\theta$分别是单调递减和递增的, 因此当$\alpha$已知时可靠度函数$R(x)$和失效率函数$H(x)$的置信水平为100(1-$\gamma$)%的置信区间分别为

取$\theta$的先验分布为Gamma分布, 其概率密度函数为

这里超参数$a>0, b>0$, $\Gamma(\cdot)$表示Gamma函数.

当$\alpha$已知时, 由(3.1)和(4.1)式, 根据Bayes公式可得到$\theta$的后验密度函数为

其中$A=\sum\limits_{i=1}^k\ln x_{i:n}+(n-k)\ln t-n\ln\alpha$.

引理4.1^[19]  设$x^*=\bigl(x_{1:n}, x_{2:n}, \cdots, x_{k:n}\bigr)$为来自某总体的双定数混合截尾数据, $\theta$为其未知参数, 则有

$\rm (i)$在平方损失函数$L_1(\theta, \delta)=(\theta-\delta)^2$下, $\theta$的Bayes估计为$\hat{\theta}_{B1}=E(\theta\mid x^*)$.

$\rm (ii)$在Linex损失函数$L_2(\theta, \delta)=e^{c(\delta-\theta)}-c(\delta-\theta)-1, (c\neq 0)$下, $\theta$的Bayes估计为

$\rm (iii)$在广义熵损失函数$L_3(\theta, \delta)\propto \bigl(\frac{\delta}{\theta}\bigr)^q-q\ln\bigl(\, \frac{\delta}{\theta}\, \bigr)-1$下, $\theta$的Bayes估计为

这里$\delta$是未知参数$\theta$的一个估计.

Linex损失和广义熵损失函数都是非对称损失函数, 过高估计和过低估计带来的损失是不一样的. 从引理4.1可以看出：对于Linex损失函数, 当$c<0$时, 过低估计的损失大于过高估计的损失, 当$c>0$时, 正好相反; 对于广义熵损失函数, 当$q>0$时, 过高估计的损失大于过低估计的损失, 当$q<0$时, 正好相反.

定理4.1  设$x^*=(x_{1:n}, x_{2:n}, \cdots, x_{k:n})$为来自两参数Pareto分布$(1.1)$的双定数混合截尾数据, 当$\alpha$已知时, 若$\theta$的先验分布由$(4.1)$式给出, 则有下列结论

$\rm (i)$在平方损失函数下, $\theta$的Bayes估计为

$\rm (ii)$当$c>-(A+b)$时, 在Linex损失函数下$\theta$的Bayes估计为

$\rm (iii)$当$q<k+a$时, 在广义熵损失函数下$\theta$的Bayes估计为

其中$A=\sum\limits_{i=1}^k\ln x_{i:n}+(n-k)\ln t-n\ln\alpha$.

证  (i)在平方损失函数下, $\theta$的Bayes估计为其后验分布的均值, 因此

(ii)

因此在Linex损失函数下, $\theta$的Bayes估计为

(iii)

因此在广义熵损失函数下, $\theta$的Bayes估计为

定理4.1得证.

利用定理4.1的结论, 要得到参数$\theta$的Bayes估计值就需要确定超参数$a, b$的值. 如果没有先验信息可以利用, 为了得到较稳健的Bayes估计可以考虑求参数$\theta$的E-Bayes估计. 根据文献[20]中关于E-Bayes估计的定义, 需要利用超参数$a, b$的先验密度函数并对$\theta$的Bayes估计再求一次数学期望.

根据文献[20]的研究成果, 可取超参数$a, b$的联合先验密度函数为

因此在平方损失函数下, $\theta$的E-Bayes估计为

当$c>-A$时, 在Linex损失函数下$\theta$的E-Bayes估计为

当$q<k$时, 在广义熵损失函数下$\theta$的E-Bayes估计为

定理4.2  设$x^*=(x_{1:n}, x_{2:n}, \cdots, x_{k:n})$为来自两参数Pareto分布$(1.1)$的双定数混合截尾数据, 当$\alpha$已知时, 若$\theta$的先验分布由$(4.1)$式给出, 则有下列结论

$\rm (i)$在平方损失函数下, 可靠度函数$R(x)$的Bayes估计为

$\rm (ii)$在Linex损失函数下, 可靠度函数$R(x)$的Bayes估计为

$\rm (iii)$当$q<\frac{A+b}{\ln x-\ln\alpha}$时, 在广义熵损失函数下可靠度函数$R(x)$的Bayes估计为

其中$A=\sum\limits_{i=1}^k\ln x_{i:n}+(n-k)\ln t-n\ln\alpha$, $x>\alpha.$

证  (i)在平方损失函数下, 可靠度函数$R(x)$的Bayes估计为

(ii)

因此在Linex损失函数下, 可靠度函数$R(x)$的Bayes估计为

(iii)

因此在广义熵损失函数下, 可靠度函数$R(x)$的Bayes估计为

定理4.2证毕.

定理4.3  设$x^*=(x_{1:n}, x_{2:n}, \cdots, x_{k:n})$为来自两参数Pareto分布$(1.1)$的双定数混合截尾数据, 当$\alpha$已知时, 若$\theta$的先验分布由$(4.1)$式给出, 则有下列结论

$\rm (i)$在平方损失函数下, 失效率函数$H(x)$的Bayes估计为

$\rm (ii)$当$c>-(A+b)x$时, 在Linex损失函数下失效率函数$H(x)$的Bayes估计为

$\rm (iii)$当$q<k+a$时, 在广义熵损失函数下失效率函数$H(x)$的Bayes估计为

其中$A=\sum\limits_{i=1}^k\ln x_{i:n}+(n-k)\ln t-n\ln\alpha$, $x>\alpha.$

证  (i)

(ii)

(iii)

则

定理4.3证毕.

取参数$\alpha$的无信息先验分布为$\pi(\alpha)=\frac{1}{\alpha}, \quad \alpha>0.$

取参数$\theta$的先验分布为指数分布, 概率密度函数为

其中超参数$\lambda>0$.

根据Bayes公式可得参数$\alpha$, $\theta$的联合后验分布为

其中$B=\sum\limits_{i=1}^{k}\ln x_{i:n}+(n-k)\ln t.$

根据(5.2)式可以得到$\theta$的后验分布为Gamma分布, 概率密度函数为

$\alpha$的后验密度函数为

因此在平方损失函数下, 参数$\theta$的Bayes估计为

在平方损失函数下, 参数$\alpha$的Bayes估计为

由于(5.5)式中的定积分不易求出, 因此$\alpha$的Bayes估计没有显式表达式, 但是可以进行数值计算求出近似值.

下面介绍一种采用随机模拟求参数$\alpha$的Bayes估计的方法.

当$\theta$给定时, 参数$\alpha$的条件后验分布为

在平方损失下, $\alpha$的Bayes估计为

由于(5.6)式中的$\theta$是未知的, 因此并不能求出具体值, 可以借助$\theta$的后验分布(5.3)式利用统计模拟计算出近似值, 具体步骤如下:

Step 1: 产生$N$个服从Gamma分布(5.3)式的随机数, 记为$\theta^{(1)}, \theta^{(2)}, \cdots, \theta^{(N)}$;

Step 2: 计算出参数$\alpha$的Bayes估计为$\hat{\alpha}_B\approx\frac{1}{N}\sum\limits_{i=1}^{N}\frac{n\theta^{(i)}}{n\theta^{(i)}+1}x_{1:n}.$

另外, 在平方损失下, 可靠度函数$R(x)$的Bayes估计为

在平方损失下, 失效率函数$H(x)$的Bayes估计为

对于多参数分布总体, 要求出参数及可靠性指标的Bayes估计, 通常要计算较复杂的多重积分, 在很多情况下都没有显式表达式, 可以选择采用Gibbs抽样法、Lindley方法和TK方法求出Bayes估计的近似值. 实际上, 超参数$\lambda$是未知的, 从而导致上述Bayes估计不能得到具体的值, 因此需要确定超参数的值. 本文利用参数的极大似然估计给出超参数$\lambda$的估计. 根据$\theta$的先验分布(5.1)式可得$E(\theta)=1/\lambda$, 而$\theta$的极大似然估计为

令$\hat{\theta}=1/\lambda$, 则可以得到超参数$\lambda$的估计为

由于分布函数$F(x;\alpha, \theta)\sim U(0, 1)$, 则$1-F(x;\alpha, \theta)\sim U(0, 1)$. 设$r_1, r_2, \cdots, r_n$是来自均匀分布$U(0, 1)$的相互独立随机数, 给定参数$\alpha, \theta$的值, 根据直接抽样法, 则$x_i=\alpha r_i^{-1/\theta},$$i=1, 2, \cdots, n,$即为来自两参数Pareto分布(1.1)的随机数. 在小样本量($n=20$)、中样本量($n=30$)和大样本量($n=50$)下给定$m_1, m_2, t_0$的值; 若$x_{m_1:n}\geq t_0$, 可以得到定数截尾样本$(x_{1:n}, x_{2:n}, \cdots, x_{m_1:n})$; 若$x_{m_1:n}<t_0$则得到定数截尾样本$(x_{1:n}, x_{2:n}, \cdots, x_{m_2:n})$. 超参数的值取$(a, b)=(2, 1)$且在损失函数中取$c=1, q=1,$利用所得到的定数截尾样本计算出各种估计值. 以上过程重复10000次, 可以得到各种点估计的均值和平均相对误差(括号中), 其中相对误差的计算公式为$RE(\hat{\beta})=|\hat{\beta}-\beta|/\beta$. 模拟计算结果如表 1–4所示.

比较表 1–4中的平均相对误差可以看出, 各种估计量的Bayes估计都要优于相应的极大似然估计, 特别是在小样本场合Bayes估计的优势更加明显. 随着样本量$n$的增大, 参数及可靠性指标的估计的平均相对误差逐渐减小.

下面分析文献[21]中第3.1节中的例子, 其中的失效数据服从Ⅱ型Pareto分布, 经过变换后可以得到来自两参数Pareto分布(1.1)的失效数据, 按从小到大的顺序排列如下:

根据本文中的结论, 我们能够计算得到未知参数、可靠度及失效率的估计, 相关结果列于表 5中.

从表 5可以看到, 尺度参数$\alpha$的极大似然估计与其Bayes估计比较接近, 随着$m_1, m_2$的增大, 形状参数$\theta$的极大似然估计与其Bayes估计的差距逐渐变小.

由于两参数Pareto分布在可靠性工程中的重要地位, 近年来产生了许多关于此分布的研究成果. 在本文中首次提出了双定数混合截尾试验方案, 并在此类截尾样本下讨论了两参数Pareto分布的参数及可靠性指标的估计问题. 利用经典统计理论计算出极大似然估计值, 根据Bayes理论在三种不同损失函数下得到了参数及可靠度指标的估计. 通过随机模拟得到各种估计量的均值和平均相对误差, 结果表明各种估计的Bayes估计都要优于相应的极大似然估计, 估计量的精度随着样本量的增大而逐渐提高. 另外也可以利用本文中提出的截尾试验模型对其它的寿命分布进行统计分析.