成败数据下ZZ分布参数的估计与性质

Estimation and Properties of ZZ Distribution Parameters under Success or Failure Data

摘要: 本文在假定元件的寿命服从一种新型分布-ZZ分布的前提下，基于成败数据，运用频率方法和Bayes方法，分别给出了该分布参数的估计，并讨论了估计的性质，同时进行数值模拟。

Abstract: In this paper, assuming that the life of the components follows a new distribution ZZ distribution, based on the success or failure data, the estimation of the distribution parameters is given by using frequency method and Bayes method, and the properties of the estimation are discussed. At the same time, numerical simulation is carried out.

1. 引言

通常元件的寿命服从以下常见的几种分布：指数分布、威布尔分布、极值分布、对数正态分布等 [1]。1976年，Martz和Waller两位作者针对单参数指数寿命型分布在无失效数据条件下相关参数以及可靠性指标的评估和验证问题，提出使用Bayes方法 [2]。在接下来的一段时间里，人们的研究工作还只限于单参数指数分布在定数截尾条件下和有失效数据的定时截尾条件下进行。比如Umesh S.等人基于特殊置信度的设定对指数分布的参数进行贝叶斯估计 [3]；Nand K. S.等人对混合指数失效分布进行了相关的贝叶斯分析 [4]。1978年，Padget W. J.等人基于可靠性相关经验及前者的研究成果，将理论方法应用到正态分布中，对单参数正态分布失效模型进行贝叶斯估计 [5]；Soman R. P. and Misra K. B.、Calabria R. and Pulcini G.等人也着力于研究当存在失效数据时，基于可靠性相关经验，对威布尔分布中的相关参数进行了贝叶斯估计，具体研究内容及结果见文献 [6] [7]。近十几年来，很多的学者专家们也针对无失效数据下威布尔分布的可靠性评估问题进行了研究，Toskos等人给出基于贝叶斯可靠性方法的理论与仿真 [8]；Erto等人对双参数威布尔分布进行相关统计推断，对危险率等进行贝叶斯决策分析 [9] [10] [11]。

事实证明，有些元件的寿命用上述四种分布来刻画是合理的；但还有一部分元件的寿命若用上述分布来刻画会与事实不符。例如，用给定的样本数据对某品牌无人机的元件寿命进行检验时，发现四种寿命分布假设都不拒绝，但对指数分布、威布尔分布、对数正态分布的平均寿命的估计却与事实相差甚远；对于极值分布，平均寿命的估计略接近实际，但也远远高于实际寿命，这说明该类元件寿命分布并不属于熟知的这四种分布。

鉴于这种情况，2019年文献 [12] 首次提出了一种新型的寿命分布函数 $ZZ\left(m,\eta \right)$ 分布。在假设元件寿命服从 $ZZ\left(m,\eta \right)$ 分布的前提下，基于截尾样本，给出了参数的最优线性无偏估计与简单线性无偏估计，以及近似的置信区间，进一步又给出了平均寿命和可靠寿命的估计。

2020年，雷露、张国志等 [13] 针对无失效数据的情况，对 $ZZ\left(m,\eta \right)$ 分布进行可靠性分析，通过对 $ZZ\left(m,\eta \right)$ 分布可靠度函数进行变换，利用其凹凸性质，得到各检测时刻可靠度之间的关系。基于可靠度的贝叶斯估计，对各个时刻和该时刻下可靠度的相关关系进行研究，利用回归分析得到 $ZZ\left(m,\eta \right)$ 分布参数的估计。同时在无失效数据场合下，给出了 $ZZ\left(m,\eta \right)$ 分布可靠度函数的置信水平为 $1-a$ 的最优置信下限表达式，在几种特殊情况下，给出了便于使用的简化形式，最后又对结果进行了数值模拟。

本文将讨论成败数据下，ZZ分布参数的估计问题。

2. ZZ分布的参数估计

2.1. 首先介绍ZZ分布及其相关结论

定义1 [12] 若随机变量X的分布函数为

$F\left(t\right)=1-\exp \left\{1-{\text{e}}^{{\left(\frac{t}{\eta }\right)}^m}\right\}$， $t>0\left(m>0,\eta >0\right)$

称X服从参数为 $m,\eta$ 的ZZ分布，记为 $X~ZZ\left(m,\eta \right)$。

引理2 [14] 若 $X~Be\left(a,b\right)$，则 $\frac{b}{a}\frac{X}{1-X}~F\left(2a,2b\right)$。

引理3 设 $B_i$，则有

${\displaystyle \underset{0}{\overset{c}{\int }}{x}^{a-1}{\left(1-x\right)}^{b-1}\text{d}x}=\frac{\Gamma \left(a\right)\Gamma \left(b\right)}{\Gamma \left(a+b\right)}F\left(2a,2b,\frac{b}{a}\frac{c}{1-c}\right)$

其中 $F\left(2a,2b,x\right)$ 为 $F\left(2a,2b\right)$ 分布的分布函数。

证明 设 $X~Be\left(a,b\right)$，

2.2. 基于频率方法的参数估计及性质

模型

假设元件的寿命 $X~zz\left(m,\eta \right)$，即分布函数为

$F\left(t\right)=1-\exp \left\{{\text{e}}^{1-{\text{e}}^{{\left(\frac{t}{\eta }\right)}^m}}\right\}\triangleq 1-\exp \left\{1-{\text{e}}^{1-{\text{e}}^{\frac{t^m}{\beta }}}\right\}$

对元件做如下成败实验，时间节点为 $t_1,t_2,\cdots ,t_l$。对每个 $t_j$ 对应的试验投入 $n_j$ 个样品，到 $t_j$ 时刻有 $f_j$ 个失效。

对于上面模型，由于

$\ln \ln \left(1-\ln \left(1-F\left(t\right)\right)\right)=m\ln t-\ln \beta$

将时间节点代入得

$\ln \ln \left(1-\ln \left(1-F\left(t_j\right)\right)\right)=m\ln t_j-\ln \beta ,j=1,2,\cdots ,l$

用 $\frac{f_j}{n_j}$ 来估计 $P_j=F\left(t_j\right)$，如果 $\frac{f_j}{n_j}=0$，或 $\frac{f_j}{n_j}=1$，则用其修正估计替代，仍记为 $\frac{f_j}{n_j}$，那么有

$\ln \ln \left[1-\ln \left(1-\frac{f_j}{n_j}\right)\right]=m\ln t_j-\ln \beta +{\epsilon }_j,j=1,2,\cdots ,l$

其中 ${\epsilon }_j$ 为估计误差。

记

$y_j=\ln \ln \left[1-\ln \left(1-\frac{f_j}{n_j}\right)\right],x_j=\ln t_j,{\beta }_0=-\ln \beta ,{\beta }_1=m$

则

$y_j={\beta }_0+{\beta }_1{x}_j+{\epsilon }_j,j=1,2,\cdots ,l$ (1)

进一步令

$Y={\left(y_1,y_2,\cdots ,y_l\right)}^{\text{T}}$， $B={\left({\beta }_0,{\beta }_1\right)}^{\text{T}}$

$X=\left(\begin{array}{cc}1& x_1\\ 1& x_2\\ ⋮& ⋮\\ 1& x_l\end{array}\right)$， $\epsilon =\left(\begin{array}{c}{\epsilon }_1\\ {\epsilon }_2\\ ⋮\\ {\epsilon }_l\end{array}\right)$

则(1)式变为 $Y=XB+\epsilon$，则B的简单最小二乘估计为

$\hat{B}={\left(X^{\text{T}}X\right)}^{-1}{X}^{\text{T}}Y\triangleq {\left({\hat{\beta }}_0,{\hat{\beta }}_1\right)}^{\text{T}}$

那么有

$\hat{\beta }={\text{e}}^{-{\hat{\beta }}_0},\hat{m}={\hat{\beta }}_1$ (2)

可以证明该估计具有渐近正态性。记 ${\hat{p}}_j=\frac{f_j}{n_j}$

由中心极限定理 [15]

$\sqrt{n_j}\left({\hat{p}}_j-p_j\right)\stackrel{L}{\to }N\left(0,p_j\left(1-p_j\right)\right);j=1,2,\cdots ,l$

假设

$\frac{n}{n_j}\triangleq \frac{\underset{1\le j\le l}{\min }\left\{n_j\right\}}{n_j}\to {\lambda }_j^2\left(n\to \infty \right)$

则

$\sqrt{n}\left({\hat{p}}_j-p_j\right)\stackrel{L}{\to }N\left(0,{\lambda }_j^2{p}_j\left(1-p_j\right)\right)$

进而有

$\sqrt{n}\left(\left(\begin{array}{c}{\hat{p}}_1\\ {\hat{p}}_2\\ ⋮\\ {\hat{p}}_l\end{array}\right)-\left(\begin{array}{c}{p}_1\\ p_2\\ ⋮\\ p_l\end{array}\right)\right)\stackrel{L}{\to }N\left(0,\Sigma \right)$

这里 $\Sigma =diag\left({\lambda }_1^2{p}_1\left(1-p_1\right),{\lambda }_2^2{p}_2\left(1-p_2\right),\cdots ,{\lambda }_l^2{p}_l\left(1-p_l\right)\right)$。那么有

$\sqrt{n}\left(y_j-\left({\beta }_0+{\beta }_1{x}_j\right)\right)\stackrel{L}{\to }N\left(0,{\sigma }_j^2\right)$

其中

$\begin{array}{c}{\sigma }_j^2={\left({\frac{\partial \ln \ln \left(1-\ln \left(1-x\right)\right)}{\partial x}|}_{x=p_j}\right)}^2{\lambda }_j^2{p}_j\left(1-p_j\right)\\ =\frac{{\lambda }_j^2{p}_j\left(1-p_j\right)}{\left[\left(1-p_j\right)\left(1-\ln \left(1-p_j\right)\right)\right]\ln \ln {\left[1-\ln \left(1-p_j\right)\right]}^2}\end{array}$

所以

$\sqrt{n}\left(\left(\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_l\end{array}\right)-\left(\begin{array}{c}{\beta }_0+{\beta }_1{x}_1\\ {\beta }_0+{\beta }_1{x}_2\\ ⋮\\ {\beta }_0+{\beta }_1{x}_l\end{array}\right)\right)\stackrel{L}{\to }N\left(0,{\Sigma }^{\ast }\right)$

其中

${\Sigma }^{\ast }=diag\left({\sigma }_1^2,{\sigma }_2^2,\cdots ,{\sigma }_l^2\right)$

则

$\sqrt{n}\left(\hat{B}-B\right)\stackrel{L}{\to }N\left(0,{\left(X^{\text{T}}X\right)}^{-1}{X}^{\text{T}}{\Sigma }^{\ast }X{\left(X^{\text{T}}X\right)}^{-1}\right)$

令

$\left(\begin{array}{c}m\\ \beta \end{array}\right)=\left(\begin{array}{c}{\beta }_1\\ {\text{e}}^{-{\beta }_0}\end{array}\right)\triangleq \left(\begin{array}{c}{g}_1\left({\beta }_0,{\beta }_1\right)\\ g_0\left({\beta }_0,{\beta }_1\right)\end{array}\right)$

则

$\left(\begin{array}{cc}\frac{\partial g_1}{\partial {\beta }_0}& \frac{\partial g_1}{\partial {\beta }_1}\\ \frac{\partial g_2}{\partial {\beta }_0}& \frac{\partial g_2}{\partial {\beta }_1}\end{array}\right)=\left(\begin{array}{cc}0& 1\\ -{\text{e}}^{-{\beta }_0}& 0\end{array}\right)\triangleq Q$

故

$\sqrt{n}\left(\left(\begin{array}{c}\hat{m}\\ \hat{\beta }\end{array}\right)-\left(\begin{array}{c}m\\ \beta \end{array}\right)\right)\stackrel{L}{\to }N\left(0,Q{\left(X^{\text{T}}X\right)}^{-1}{X}^{\text{T}}{\Sigma }^{\ast }X{\left(X^{\text{T}}X\right)}^{-1}{Q}^{\text{T}}\right)$

即该估计具有渐近正态性。

2.3. 基于Bayes方法的参数估计及性质

在模型中，不妨设时间节点为 $t_1<t_2<\cdots <t_l$，对每个 $t_j$ 对应的试验投入 $n_j$ 个样品，到 $t_j$ 时刻有 $f_j$ 个失效，无失效个数为

$r_j=n_j-f_j,r_j\le n_j,j=1,2,\cdots ,l$

在 $t_j$ 时刻，可靠度 $R\left(t_j\right)$ 的估计自然可取

$\hat{R}\left(t_j\right)=\frac{r_j}{n_j}\triangleq {\hat{R}}_j$

当然希望 ${\hat{R}}_1\ge {\hat{R}}_2\ge \cdots \ge {\hat{R}}_l$，但由于样本的随机性，这种序关系未必成立，这种情况称为数据倒挂。采用数据倒挂的一般方法，对每个 ${\hat{R}}_j$ 进行修正，使之序关系成立。为方便起见修正后的数据依然记为 ${\hat{R}}_j$，此时序关系 ${\hat{R}}_1\ge {\hat{R}}_2\ge \cdots \ge {\hat{R}}_l$ 成立。记

$Q_j=\left[{\hat{R}}_j{n}_j\right]\left(j=1,2,\cdots ,l\right)$

符号 $\left[x\right]$ 表示取x的整数部分。 $Q_j$ 表示 $r_j$ 的修正，可认为

$Q_j~B\left(n_j,R_j\right)\left(j=1,2,\cdots ,l\right)$

其中 $R_j=R\left(t_j\right)$。

由于 ${\hat{R}}_1\ge {\hat{R}}_2\ge \cdots \ge {\hat{R}}_l$，一旦有 ${\hat{R}}_1=1$ (无失效数出现)，若用 ${\hat{R}}_1$ 来估计 $R\left(t_1\right)$ 显然不合适，所以借鉴处理无失效数据的办法，取

${\tilde{R}}_1=\{\begin{array}{l}\frac{Q_1+0.5}{n_1+1}=\frac{r_1+0.5}{n_1+1},{\hat{R}}_1=1\\ {\hat{R}}_1,{\hat{R}}_1<1\end{array}$

由文献 [10] 可知

$R_2<\exp \left\{1-{\left(1-\ln R_1\right)}^{\frac{t_2}{t_1}}\right\}$

记

$R_2^{*}=\exp \left\{1-{\left(1-\ln {\tilde{R}}_1\right)}^{\frac{t_2}{t_1}}\right\}$

当对 $R_2$ 没有更多信息时，那么可取 $R_2$ 的先验分布是均匀分布 $U\left(0,R_2^{*}\right)$，根据 $Q_2~B\left(n_2,R_2\right)$，可得到 $R_2$ 的后验密度为

$\pi \left(R_2|Q_2\right)=\frac{R_2^{Q_2}{\left(1-R_2\right)}^{n_2-Q_2}}{{\displaystyle \underset{0}{\overset{R_2^{*}}{\int }}{R}_2^{Q_2}{\left(1-R_2\right)}^{n_2-Q_2}\text{d}{R}_2}},\left(0<R_2<R_2^{*}\right)$

所以 $R_2$ 的贝叶斯估计 ${\tilde{R}}_2$ 为

${\tilde{R}}_2=\frac{{\displaystyle \underset{0}{\overset{R_2^{*}}{\int }}{R}_2^{Q_2+1}{\left(1-R_2\right)}^{n_2-Q_2}\text{d}{R}_2}}{{\displaystyle \underset{0}{\overset{R_2^{*}}{\int }}{R}_2^{Q_2}{\left(1-R_2\right)}^{n_2-Q_2}\text{d}{R}_2}}$

由引理3可知

$\begin{array}{l}{\displaystyle \underset{0}{\overset{R_2^{*}}{\int }}{R}_2^{Q_2+1}{\left(1-R_2\right)}^{n_2-Q_2}\text{d}{R}_2}\\ =\frac{\Gamma \left(Q_2+2\right)\Gamma \left(n_2-Q_2+1\right)}{\Gamma \left(n_2+3\right)}F\left(2\left(Q_2+2\right),2\left(n_2-Q_2+1\right),\frac{n_2-Q_2+1}{Q_2+2}\cdot \frac{R_2^{*}}{1-R_2^{*}}\right)\end{array}$

$\begin{array}{l}{\displaystyle \underset{0}{\overset{R_2^{*}}{\int }}{R}_2^{Q_2}{\left(1-R_2\right)}^{n_2-Q_2}\text{d}{R}_2}\\ =\frac{\Gamma \left(Q_2+1\right)\Gamma \left(n_2-Q_2+1\right)}{\Gamma \left(n_2+2\right)}F\left(2\left(Q_2+1\right),2\left(n_2-Q_2+1\right),\frac{n_2-Q_2+1}{Q_2+1}\cdot \frac{R_2^{*}}{1-R_2^{*}}\right)\end{array}$

所以

${\tilde{R}}_2=\frac{Q_2+1}{n_2+2}\frac{F\left(2\left(Q_2+2\right),2\left(n_2-Q_2+1\right),\frac{n_2-Q_2+1}{Q_2+1}\cdot \frac{R_2^{*}}{1-R_2^{*}}\right)}{F\left(2\left(Q_2+1\right),2\left(n_2-Q_2+1\right),\frac{n_2-Q_2+1}{Q_2+1}\cdot \frac{R_2^{*}}{1-R_2^{*}}\right)}$

又由于

$R_3<\exp \left\{1-{\left(1-\ln R_2\right)}^{\frac{t_3}{t_2}}\right\}$

记

$R_3^{*}=\exp \left\{1-{\left(1-\ln {\tilde{R}}_2\right)}^{\frac{t_3}{t_2}}\right\}$

那么可取 $R_3$ 的先验分布是均匀分布 $U\left(0,R_3^{*}\right)$，根据 $Q_3~B\left(n_3,R_3\right)$，同理可得到 $R_3$ 的贝叶斯估计 ${\tilde{R}}_3$ 为

${\tilde{R}}_3=\frac{Q_3+1}{n_3+2}\frac{F\left(2\left(Q_3+2\right),2\left(n_3-Q_3+1\right),\frac{n_3-Q_3+1}{Q_3+1}\cdot \frac{R_3^{*}}{1-R_3^{*}}\right)}{F\left(2\left(Q_3+1\right),2\left(n_3-Q_{i3}+1\right),\frac{n_3-Q_3+1}{Q_3+1}\cdot \frac{R_3^{*}}{1-R_3^{*}}\right)}$

依此类推可得到 ${\tilde{R}}_{i4},{\tilde{R}}_{i5},\cdots ,{\tilde{R}}_{i{l}_i}$ 的贝叶斯估计。

易知，这些贝叶斯估计满足 $1>{\tilde{R}}_1\ge {\tilde{R}}_2\ge \cdots \ge {\tilde{R}}_l>0$。由于元件的可靠度函数为

$R\left(t\right)=\exp \left\{1-{\text{e}}^{\frac{t^m}{\beta }}\right\}$

那么有

$\ln \ln \left(1-\ln R\left(t\right)\right)=-\ln \beta +m\ln t$

用 ${\tilde{R}}_j$ 去估计 $R\left(t_j\right)$ 其误差记为 ${\epsilon }_j$，则得到

$\ln \ln \left(1-\ln {\tilde{R}}_j\right)=-\ln \beta +m\ln t_j+{\epsilon }_j\left(j=1,2,\cdots ,l\right)$ (4)

记

$Y_j=\ln \ln \left(1-\ln {\tilde{R}}_j\right)$， $Y={\left(Y_1,Y_2,\cdots ,Y_l\right)}^{\text{T}}$， $B={\left(-\ln \beta ,m\right)}^{\text{T}}$

$X=\left(\begin{array}{cc}1& \ln t_1\\ 1& \ln t_2\\ ⋮& ⋮\\ 1& \ln t_l\end{array}\right)$， $\epsilon =\left(\begin{array}{c}{\epsilon }_1\\ {\epsilon }_2\\ ⋮\\ {\epsilon }_l\end{array}\right)$

则(4)式可表示为 $Y=XB+\epsilon$，由此 $B_i$ 的最小二乘估计为

$\hat{B}={\left(X^{\text{T}}X\right)}^{-1}{X}^{\text{T}}Y$

记 ${\xi }_1={\left(1,0\right)}^{\text{T}},{\xi }_2={\left(0,1\right)}^{\text{T}}$，那么有

$\{\begin{array}{l}\hat{\beta }={\text{e}}^{-{\xi }_1^{\text{T}}\hat{B}}\\ \hat{m}={\xi }_2^{\text{T}}\hat{B}\end{array}$ (5)

2.4. 算例

参数的真值和样本容量设定数如表1。

Table 1. True value of parameters and number of test samples

表1. 参数真值和试验样品数

由Monte-Carlo模拟法，在每个设定的时间节点分别产生样本总数为n，寿命服从ZZ分布的成败数据，如表2。

Table 2. Failure frequency table of each time node

表2. 各时间节点失效频率表

表3基于Bayes方法，模拟列出各时间节点的 ${\hat{R}}_j,Q_j,{\tilde{R}}_j$ 数值。

Table 3. Value of each R ^ j , Q j , R ˜ j

表3. 各 ${\hat{R}}_j,Q_j,{\tilde{R}}_j$ 的值

由此得到参数m和 $\eta$ 的估计如表4。

Table 4. Estimation of parameters m and η

表4. 参数m和 $\eta$ 的估计

3. 结论

本文在成败数据下，研究了复杂系统可靠度的估计问题，基于两种方法得到具体结果如下：

1) 基于成败数据，采用频率估计概率的方法，针对模型 ，首先利用数据的特点及回归分析方法给出了各元件寿命分布参数的估计。然后利用复杂系统可靠度解析表达式得到了该系统的可靠度估计，并证明了该估计具有渐近正态性。

2) 基于成败数据，相比前一种方法，本方法深入考虑了可靠度函数单调非增这一特点，对可能出现“数据倒挂”现象进行了预处理。然后利用Bayes方法及回归分析方法，给出了各元件寿命分布参数的估计。两种估计方法各有长处，方法一简单且估计具有渐近正态性；方法二充分利用了可靠度函数单调性的特点，估计会好一些，但估计的分布不易获得。

3) 对上述结果进行了数值模拟，结果显示两种方法可行有效。

致谢

本论文的工作是在我的老师张国志教授的悉心指导下完成的，在此向老师表示衷心的谢意。

基金项目

黑龙江省自然科学基金资助项目，编号：A2018006。