复杂删失模型下指数威布尔分布的参数估计

Parameter Estimation of Exponential Weibull Distribution under Complex Censored Model

摘要: 针对服从指数威布尔分布的逐步II型删失模型，分别用极大似然方法和贝叶斯方法对指数威布尔分布的两个未知参数进行了估计。其中，使用Newton-Raphson方法获得了参数的极大似然估计值，使用Metropolis-Hastings算法获得了参数的贝叶斯估计值。通过蒙特卡洛方法对这两种估计方法的性能进行了评估，结果表明两种方法的估计效果都很好，且贝叶斯估计优于极大似然估计。最后，通过分析一组实际数据进一步证明了本文所提估计方法的精确性。

Abstract: The maximum likelihood method and Bayesian method are used to estimate the two unknown parameters of the exponential Weibull distribution for the progressively type II censoring model that follows the exponential Weibull distribution. The maximum likelihood estimation of param-eters is obtained by Newton-Raphson method, and the Bayesian estimation of parameters is ob-tained by Metropolis-Hastings algorithm. The performance of these two estimation methods is evaluated by Monte Carlo method. The results show that the estimation effects of the two methods are very good, and Bayesian estimation is better than maximum likelihood estimation. Finally, the accuracy of the proposed estimation method is further proved by analyzing a group of actual data.

1. 引言

威布尔分布的单调危险函数无法适应实际应用中经常出现的非单调危险率，为了解决此问题，Mudholkar等 [1] 提出了指数威布尔分布。作为威布尔分布的推广，指数威布尔分布被广泛应用于生物医学、精算、航空航天等诸多领域 [2]。近年已有大量文献对该分布在完全数据的情况下进行了参数估计的研究，桂国祥等 [3] 在线性损失下研究了指数威布尔分布的经验Bayes问题，并证明了所提的经验Bayes检验函数满足渐近最优性；杨冬霞等 [4] 在复合Mlinex损失函数下研究了指数威布尔分布的Bayes估计、E-Bayes估计和多层Bayes估计。

然而在大多数寿命试验中，考虑到效率和花费等问题，实验数据并不能被全部观测，研究者希望能撤出一部分未失效个体，所以就产生了删失数据方案。常见的删失数据方案有定时删失、定数删失、区间删失、I型删失和II型删失等。基于删失数据的统计推断问题一直被广泛研究，姜宁宁等 [5] 在定时截尾数据下，通过数据填充的方法，研究了威布尔分布的可靠度置信限；何朝兵等 [6] 在左截断右删失数据下，使用EM算法对几何分布进行了参数估计；魏艳华等 [7] 使用混合Gibbs算法，分别在分组数据和定数截尾场合对指数威布尔分布的参数进行了估计；王纯杰等 [8] 在右删失数据下对部分线性模型进行了贝叶斯P-样条估计，并将其应用于卵巢癌生存时间数据中；龙兵等 [9] 提出了一种新的双定数混合截尾方案，并在该截尾方案下研究了两参数Pareto分布的参数估计和可靠性指标估计。

在逐步删失方案下对寿命分布的研究有很多，但对指数威布尔分布的研究仍处于空白；提前给定删失方案有较大的局限性，随即移除更贴合现实情况；绝大多数关于参数估计的文献都是假设已知其中某一个参数来估计另一个未知参数，但在实际应用中，并不能提前已知某一个参数。因此，本文在带有随机移除的逐步II型删失方案下，研究了指数威布尔分布的两个未知参数的估计问题，本文其余部分内容如下：第二节介绍了本文所涉及的指数威布尔模型、逐步II型删失模型和带有随机移除的逐步II型删失模型；第三节计算了未知参数的极大似然估计；第四节使用Metropolis-Hastings算法计算了未知参数的贝叶斯估计；第五节使用蒙特卡洛模拟对本文所提的两种估计进行了比较；第六节对一组真实数据进行了分析，进一步证明了本文所提两种估计的精确性。

2. 模型简介

2.1. 指数威布尔模型

定义1 设x为随机变量，若其密度函数为：

$f\left(x\right)=\theta \alpha x^{\alpha -1}{e}^{-x^{\alpha }}{\left(1-e^{-x^{\alpha }}\right)}^{\theta -1},x>0,\theta >0,\alpha >0$ (1)

则称 $X$ 为服从参数为 $\theta$ 和 $\alpha$ 的指数威布尔分布，其分布函数为：

$F\left(x\right)={\left(1-e^{-x^{\alpha }}\right)}^{\theta },x>0,\theta >0,\alpha >0$ (2)

显然，当 $\theta =1$ 时，指数威布尔分布与威布尔分布形式相同；当 $\theta =1$， $\alpha =1$ 时，指数威布尔分布与标准指数分布形式相同。

2.2. 逐步II型删失模型

假设有 $n$ 个相同的部件同时进行寿命测试实验，在实验进行前预先指定一个小于 $n$ 的数 $m$，第一个部件失效时，记录失效时间 $x_{1:m:n}$，并在剩余的 $n-1$ 个未失效部件中随机移除 $R_1$ 个，第二个部件失效时，记录失效时间 $x_{2:m:n}$，并在剩余的 $n-2-R_1$ 个未失效部件中随机移除 $R_2$ 个，以此类推，直至第 $m$ 个部件失效时，记录失效时间 $x_{m:m:n}$，并把剩余的 $n-m-R_1-R_2-\cdots -R_{m-1}$ 个未失效部件全部移除，其中删失方案 $R_i\left(\text{i}=1,2,\cdots ,m\right)$ 为提前给定的常数，记逐步II型删失样本为 $x=\left(x_{1:m:n},x_{2:m:n},\cdots ,x_{{}_{m:m:n}}\right)$。

2.3. 带有随机移除的逐步II型删失模型

逐步II型删失模型中删失方案 $R_i$ 是提前给定的常数，但是在实际应用中，删失方案通常会受到很多因素的影响，具有很大的随机性，假设删失方案 $R_i$ 为服从概率为 $p\left(0<p<1\right)$ 的二项分布的随机变量，其余实验流程按逐步II型删失模型的实验流程进行，即可得到带有随机移除的逐步II型删失样本。

3. 极大似然估计

假设 $x=\left(x_{1:m:n},x_{2:m:n},\cdots ,x_{{}_{m:m:n}}\right)$ 为服从指数威布尔分布的带有随机移除的逐步II型删失样本，删失方案为 $R=\left(R_1,R_2,\cdots ,R_{\text{m}}\right)$，方便起见，记删失样本为 $x=\left(x_1,x_2,\cdots ,x_m\right)$，删失样本 $x$ 基于删失方案 $R$ 的似然函数为：

$\begin{array}{c}L\left(\theta ,\alpha |x\right)=C{\displaystyle \underset{i=1}{\overset{m}{\prod }}f\left(x_i\right){\left[1-F\left(x_i\right)\right]}^{R_i}}\\ =C{\theta }^m{\alpha }^m{{\displaystyle \underset{i=1}{\overset{m}{\prod }}{x}_i^{\alpha -1}{e}^{-x_i^{\alpha }}{\left(1-e^{-x_i^{\alpha }}\right)}^{\theta -1}\left[1-{\left(1-e^{-x_i^{\alpha }}\right)}^{\theta }\right]}}^{R_i}\end{array}$ (3)

其中， $C=n\left(n-r_1-1\right)\left(n-r_1-r_2-2\right)\cdots \left(n-r_1-r_2-\cdots -r_m-m+1\right)$。

因此，对数似然函数为：

$\begin{array}{c}\ln L\left(\theta ,\alpha \right)=\ln C+m\ln \theta +m\ln \alpha +\left(\alpha -1\right){\displaystyle \underset{i=1}{\overset{m}{\sum }}\ln x_i}-{\displaystyle \underset{i=1}{\overset{m}{\sum }}{x}_i^{\alpha }}\\ +\left(\theta -1\right){\displaystyle \underset{i=1}{\overset{m}{\sum }}\ln \left(1-e^{-x_i^{\alpha }}\right)}+{\displaystyle \underset{i=1}{\overset{m}{\sum }}{R}_i\ln \left[1-{\left(1-e^{-x_i^{\alpha }}\right)}^{\theta }\right]}\end{array}$ (4)

对式 分别关于 $\theta$ 和 $\alpha$ 求偏导并令其等于0，可得：

$\frac{\partial \ln L\left(\theta ,\alpha \right)}{\partial \theta }=\frac{m}{\theta }+{\displaystyle \underset{i=1}{\overset{m}{\sum }}\ln \left(1-e^{-x_i^{\alpha }}\right)}-{\displaystyle \underset{i=1}{\overset{m}{\sum }}{R}_i\frac{{\left(1-e^{-x_i^{\alpha }}\right)}^{\theta }\ln \left(1-e^{-x_i^{\alpha }}\right)}{1-{\left(1-e^{-x_i^{\alpha }}\right)}^{\theta }}=0}$ (5)

$\begin{array}{c}\frac{\partial \ln L\left(\theta ,\alpha \right)}{\partial \alpha }=\frac{m}{\alpha }+{\displaystyle \underset{i=1}{\overset{m}{\sum }}\ln x_i}-{\displaystyle \underset{i=1}{\overset{m}{\sum }}{x}_i^{\alpha }\ln x_i}\\ +\left(\theta -1\right){\displaystyle \underset{i=1}{\overset{m}{\sum }}\frac{e^{-x_i^{\alpha }}{x}_i^{\alpha }\ln x_i}{1-e^{-x_i^{\alpha }}}}-{\displaystyle \underset{i=1}{\overset{m}{\sum }}{R}_i\frac{\theta {\left(1-e^{-x_i^{\alpha }}\right)}^{\theta -1}{e}^{-x_i^{\alpha }}{x}_i^{\alpha }\ln x_i}{1-{\left(1-e^{-x_i^{\alpha }}\right)}^{\theta }}}=0\end{array}$ (6)

参数 $\theta$ 和 $\alpha$ 的极大似然估计值 ${\hat{\theta }}_{ML}$ 和 ${\hat{\alpha }}_{ML}$ 可通过同时求解非线性方程(5)和(6)获得，显然方程(5)和(6)很难得到解析解。为了解决这个问题，可在R软件中使用Newton-Raphson方法获得参数的极大似然估计。

4. 贝叶斯估计

与极大似然估计不同，贝叶斯估计不完全依赖于观测到的样本数据，在此基础上，贝叶斯估计结合了样本的先验信息进行推理和分析，因此，贝叶斯方法能更客观、合理地描述参数。选择适当的先验分布，可以提高贝叶斯方法下估计的准确性，一般来说，由于形式简单，共轭先验分布是首选。然而，当 $\theta$ 和 $\alpha$ 都未知时，参数的联合共轭先验不存在，鉴于以上情况，Nassar [10] 建议采用如下形式的二元先验密度函数：

${\pi }_1\left(\theta |\alpha \right)=\frac{{\alpha }^{-\text{v}}}{\Gamma \left(v\right)}{\theta }^{v-1}{e}^{-\frac{\theta }{\alpha }},\theta >0$

${\pi }_2\left(\alpha \right)=\frac{1}{b}{e}^{-\frac{\alpha }{b}},\alpha >0$

其中， ${\pi }_1\left(\theta |\alpha \right)$ 是假设 $\alpha$ 和 $v$ 都已知的伽马先验密度函数， ${\pi }_2\left(\alpha \right)$ 是假设 $b$ 已知的指数密度函数。因此， $\alpha$ 和 $\theta$ 的二元先验密度函数可以写成：

$\pi \left(\alpha ,\theta \right)={\left(b\Gamma \left(v\right)\right)}^{-1}{\alpha }^{-v}{\theta }^{v-1}{e}^{-\left({\alpha }^2+b\theta \right)/b\alpha },\theta >0,\alpha >0$ (7)

结合式(3)和式(7)可知参数 $\alpha$ 和 $\theta$ 的联合后验密度函数为：

$\begin{array}{c}\pi \left(\alpha ,\theta |x\right)=\frac{{\pi }_1\left(\theta |\alpha \right){\pi }_2\left(\alpha \right)L\left(\theta ,\alpha |x\right)}{{\displaystyle {\int }_0^{\infty }{\displaystyle {\int }_0^{\infty }{\pi }_1\left(\theta |\alpha \right){\pi }_2\left(\alpha \right)L\left(\theta ,\alpha |x\right)d\theta d\alpha }}}\\ \propto {\alpha }^{m-v}{\theta }^{m+v-1}{e}^{-\left({\alpha }^2+b\theta \right)/b\alpha }{{\displaystyle \underset{i=1}{\overset{m}{\prod }}{x}_i^{\alpha -1}{e}^{-x_i^{\alpha }}{\left(1-e^{-x_i^{\alpha }}\right)}^{\theta -1}\left[1-{\left(1-e^{-x_i^{\alpha }}\right)}^{\theta }\right]}}^{R_i}\end{array}$ (8)

从式(8)可以看出，贝叶斯估计量不能以闭合形式获得，所以使用数值方法来计算贝叶斯估计值。Gibbs抽样是一种常用的估计特定分布属性的方法，可将复杂且难以计算的期望简单化。Gibbs抽样下的Metropolis-Hastings (MH)算法更容易从后验分布生成样本，该算法建立了一个稳定分布的马尔科夫链，并使其达到稳定状态，所以MH算法易于实现，并且有助于降低高维分布的操作复杂性。本文将使用MH算法获得参数的贝叶斯估计值，算法具体过程如下：

步骤1 给定参数的初始值： $\left({\theta }_0,{\alpha }_0\right)$ ；

步骤2 取建议分布为正态分布，利用二元正态分布 $N_2\left(\left({\theta }_{q-1},{\alpha }_{q-1}\right),{\displaystyle \sum {}_{MH}}\right)$ 生成候选值 $\left({\theta }^{\ast },{\alpha }^{\ast }\right)$，其中 ${\displaystyle \sum {}_{MH}}$ 是参数的方差–协方差矩阵， $q=1,2,\cdots ,N$ ；

步骤3 构造判断标准： $\gamma =\min \left\{1,\frac{\pi \left({\theta }^{\ast },{\alpha }^{\ast }|x\right)}{\pi \left({\theta }_{q-1},{\alpha }_{q-1}|x\right)}\right\}$，其中 $\pi \left(\cdot |x\right)$ 是对应的联合后验密度函数；

步骤4 利用均匀分布 $U\left(0,1\right)$ 生成 $u$ ；

步骤5 使用判断标准进行比较和筛选：如果 $u<\gamma$，则接受该提案，并令 $\left({\theta }_q,{\alpha }_q\right)=\left({\theta }^{\ast },{\alpha }^{\ast }\right)$，否则，令 $\left({\theta }_q,{\alpha }_q\right)=\left({\theta }_{q-1},{\alpha }_{q-1}\right)$ ；

步骤6 将步骤2~5重复 $N$ 次，获得样本 $\left({\theta }_1,{\theta }_2,\cdots ,{\theta }_N\right)$ 和 $\left({\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_N\right)$ ；

步骤7 为了保证收敛性并消除初值选择的影响，去除前 $N_0$ 个样本；

步骤8 使用剩余的 $N-N_0$ 个样本计算贝叶斯估计值： ${\hat{\theta }}_{MH}=\frac{1}{N-N_0}{\displaystyle \underset{q=N_0+1}{\overset{N}{\sum }}{\theta }_q}$， ${\hat{\alpha }}_{MH}=\frac{1}{N-N_0}{\displaystyle \underset{q=N_0+1}{\overset{N}{\sum }}{\alpha }_q}$。

5. 数值模拟

蒙特卡洛方法是一种重要的统计分析方法，在本节中使用该方法评估前面章节所提估算方法的性能和效率。首先，按照Wang等 [11] 提出的方法，推导了服从指数威布尔分布的带有随机移除的逐步II型删失样本的算法：

步骤1 利用二项分布 $R_i\sim B\left(n-m-{\displaystyle \underset{j=1}{\overset{i-1}{\sum }}{R}_j},p\right)$ 生成每次有部件失效时移除的部件数 $R_1,R_2,\cdots ,R_m$ ；

步骤2 利用均匀分布 $U\left(0,1\right)$ 生成 $m$ 个独立同分布的随机变量 $U_1,U_2,\cdots ,U_m$ ；

步骤3 计算 $V_i=U_i^{1/\left(i+R_m+R_{m-1}+\cdots +R_{m-i+1}\right)}$， $i=1,2,\cdots ,m$ ；

步骤4 计算 $W_i=1-V_m{V}_{m-1}\cdots V_{m-i+1}$， $i=1,2,\cdots ,m$ ；

步骤5 对于给定的 $\theta$ 和 $\alpha$，计算 $x_i=F^{-1}\left(V_i\right)={\left[-\ln \left(1-y^{1/\theta }\right)\right]}^{1/\alpha }$， $i=1,2,\cdots ,m$，则 $x_1,x_2,\cdots ,x_m$ 为服从指数威布尔分布的带有随机移除的逐步II型删失样本。

基于以上随机样本，通过R软件可得到参数的极大似然估计值(MLE)和贝叶斯估计值(BE)。本文选择了不同的样本量和有效样本量来更全面地研究所提方法的性能。为了比较极大似然估计和贝叶斯估计的性能，分别计算了估计值的平均绝对偏差(AB)和均方误差(MSE)。估计过程重复M次，则平均绝对偏差为 $\frac{1}{M}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\left|\hat{\theta }-\theta \right|}$，均方误差为 $\frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\left(\hat{\theta }-\theta \right)}^2}$。取参数 $p=0.3$ 、 $b=2$ 和 $v=2$，将极大似然估计过程重复1000次，考虑MH算法的复杂性和算法程序运行的效率，将贝叶斯估计过程重复100次。模拟结果见表1。在模拟中，用来生成删失方案的参数 $p$ 对估计值的影响不大，所以表1只展示 $p=0.3$ 时的估计值。

Table 1. Maximum likelihood estimation and Bayesian estimation of parameters

表1. 参数的极大似然估计和贝叶斯估计

由表1的模拟结果可得出以下结论：

1) 绝大多数估计值的平均绝对偏差和均方误差都很小，说明极大似然估计和贝叶斯估计的效果都很好。

2) 由均方误差可知 $\alpha$ 的估计性能优于 $\theta$ 的估计性能；由平均绝对偏差可知 $\alpha$ 的估计值比 $\theta$ 的估计值更接近真实值。

3) 随着n或m的增加，所有估计值的平均绝对偏差和均方误差均减小，即有效数据越多，参数的模拟效果越好。

4) 贝叶斯估计优于极大似然估计。贝叶斯估计值有更小的平均绝对偏差和均方误差，这是因为与极大似然方法相比，贝叶斯方法不仅考虑了样本信息，还考虑了未知参数的先验信息。

5) 对于不适当的先验分布和超参数的选取，贝叶斯估计可能会有较大的偏差，模拟结果表明，本文对先验分布和超参数的选取是合理的。

6. 实例分析

本节用Nichols等 [12] 提供的一组真实数据来说明本文提出的估算方法的实际应用价值。数据给出了碳纤维断裂应力的100个观察值(单位：Gba)：3.7，2.74，2.73，2.5，3.6，3.11，3.27，2.87，1.47，3.11，4.42，2.41，3.19，3.22，1.69，3.28，3.09，1.87，3.15，4.9，3.75，2.43，2.95，2.97，3.39，2.96，2.53，2.67，2.93，3.22，3.39，2.81，4.2，3.33，2.55，3.31，3.31，2.85，2.56，3.56，3.15，2.35，2.55，2.59，2.38，2.81，2.77，2.17，2.83，1.92，1.41，3.68，2.97，1.36，0.98，2.76，4.91，3.68，1.84，1.59，3.19，1.57，0.81，5.56，1.73，1.59，2，1.22，1.12，1.71，2.17，1.17，5.08，2.48，1.18，3.51，2.17，1.69，1.25，4.38，1.84，0.39，3.68，2.48，0.85，1.61，2.79，4.7，2.03，1.8，1.57，1.08，2.03，1.61，2.12，1.89，2.88，2.82，2.05，3.65。Pal等 [13] 证明了该组数据与指数威布尔分布的拟合度很高，并使用求解似然方程的方法对指数威布尔分布的参数进行了估计，估计值为： ${\theta }_p=1.17262$， ${\alpha }_P=2.57902$。

基于这组真实数据，采用本文提出的两种估算方法对指数威布尔分布的参数进行估计，估计过程重复100次，参数的估计值取其平均值，分别在有效样本为50、65、80和95的情况下进行参数估计，取参数 $p=0.3$ 、 $b=2.5$ 和 $v=0.5$，模拟结果见表2。

Table 2. Parameter estimation of carbon fiber fracture stress data under different scenarios

表2.碳纤维断裂应力数据在不同方案下的参数估计

表2中的估计值与Pal计算出的估计值差距很小，说明了估计的精确性；同时表2中极大似然估计值与贝叶斯估计值相差不大，这也侧面印证了本文所提估计方法的正确性。最后，表2进一步验证了数值模拟所得的结论。

基金项目

国家自然科学基金项目(11801488)；新疆师范大学教学研究与改革项目(SDJG2020-30)；新疆师范大学科研发展专项项目(XJNUZX202001)。