Reliability Assessment of Heavily Censored Data Based on E-Bayesian Estimation

Abstract

The classic E-Bayesian estimation methods can only derive point estimation of the reliability parameters. In this paper, an improved E-Bayesian estimation method is proposed to evaluate product reliability under heavily censored data, which can achieve both point and confidence interval estimation for the reliability parameters. Firstly, by analyzing the concavity & convexity and function characteristics of the Weibull distribution, the value of product failure probability is limited to a certain range. Secondly, an improved weighted least squares method is utilized to construct the confidence interval estimation model of reliability parameters. Simulation results show that the proposed approach can significantly improve the calculation speed and estimation accuracy with just very few robustness reductions. Finally, a real-world case study of the sun gear transmission mechanism is used to validate the effectiveness of the presented method.

1. Introduction

The reliability evaluation method of the product is essentially a mathematical statistical method of specific censored data under specific distribution occasions. Normal distribution, exponential distribution and Weibull distribution are often used to describe the life distribution of products. Among them, Weibull distribution is the most flexible one. Weibull distribution was first proposed by Swedish scientist Weibull. The Weibull distribution model can be regarded as a generalization of the exponential distribution model. The introduction of shape parameters makes it more accurate than the exponential distribution model in describing the products whose failure rate increases or decreases with time. In addition, by adjusting the shape parameters, Weibull distribution is also capable of describing the normal distribution approximately. Many studies show that the life distribution of most mechanical products approximately follows Weibull distribution. Wang [1] proposed a new inference for constant-stress accelerated life tests with Weibull distribution. Wang [2], Joarder [3] and Kundu [4] proposed MLE methods for Weibull parameters. Denecke [5], Tan [6] and Krishnamoorthy [7] studied the method of solving lower confidence limit for Weibull distribution.

Reliability life test is an important way to collect life data. If a test stops before all samples fail, we call it censoring test, and the corresponding test data is truncated sample. The censoring test can be classified into four types, namely conventional Type-II censoring, progressive Type-II censoring [8], conventional Type-I censoring and progressive Type-I censoring [9].

With the rapid development of science and technology, products with long lifetime and high reliability are emerging. As a result, in many life tests only the heavily censored data are available for reliability estimation. Under some special conditions, including products with high reliability and very few test samples, no test data or only very few failure data can be obtained, which challenges the traditional reliability evaluation methods [10,11].

Jiang [12] studied the problem of reliability estimation with zero failure data, in which the Bayesian prior distribution model was constructed by using the “kernel” idea. Jiang [13] developed a hybrid censoring index to quantitatively describe censoring characteristics of a data set, proposed a novel parameter estimation method based on information extracted from censored observations, and evaluated the accuracy and robustness of the proposed method through a numerical experiment. Han [14,15,16] reviewed the reliability assessment methods under zero-failure data conditions, and summarized the corresponding advantages, disadvantages and application ranges. Li [17,18] proposed a reliability assessment method with revised confidence limits. It can effectively avoid the aggressive phenomenon of reliability estimation results. The hierarchical Bayesian method generalizes the current Bayesian method to provide the ability to perform parameter point estimation and confidence interval estimation at the same time, which effectively improves the credibility of the results [19]. Based on the modification of the hierarchical Bayesian estimation method, the E-Bayesian estimation method was proposed. Compared with the hierarchical Bayesian method, the E-Bayesian estimation method simplifies the mathematical model and improves the calculation efficiency [20,21,22]. Jia [23] proposed an improved method based on the Bayesian inference and least-squares method, and the four existing methods were com-pared with this method in terms of applicability, precision, efficiency, robustness, and simplicity.

The current E-Bayesian estimation methods mainly focus on point estimation of parameters, which do not solve the problem of confidence interval estimation. If the calculation method of confidence interval estimation is different from point estimation, the credibility of the results will reduce [24,25]. In addition, most E-Bayesian estimation methods assume that the range of failure probability is (0, 1), it does not consider the relationship between the failure probability of each truncation time, which causes the range of failure probability values to be too large, and affects the accuracy. Furthermore, in the process of solving the distribution parameters by using the distribution curve, the weight of the weighted least square method used does not consider the role of the failure sample, which makes the obtained reliability higher. In order to overcome the above problems, we propose a new reliability assessment method based on the E-Bayesian. Firstly, based on the E-Bayesian estimation and the distribution curve method, we develop a confidence interval method. Secondly, we use the unevenness of the distribution function and the function characteristics to determine the value range of the failure probability. Thirdly, improve the original weighted least square method with a new weight function, which consider the number of failed samples.

2. Heavily Censored Data Model

Due to time and cost constraints, for many products their life tests are often terminated before all units fail, and thus produce the so-called censored data. In the test, due to the high reliability of the product and the large number of units, the test observation time is relatively short. Therefore, in the collected data, the ratio of faults to the entire unit is usually less than 0.5 or very small. This collected data is called heavily censored data [13]. A model for heavily censored data is established, and it is the foundation for subsequent product reliability index evaluation.

The Weibull distribution is a common random distribution and widely used in modeling lifetime distribution for machinery and electronic products [26,27]. If we set the shape parameter of Weibull distribution to be 1, it reduces to an exponential distribution, which is only suitable to model electronic lifetime with the property of memoryless. Under the condition of heavily censored data, it is assumed that the product life $T$ follows a distribution function $T~WE\left(m,\eta \right)$, where there are $N$ ($N>k$) samples in the timing truncation test, a total of $k$ time timing truncation tests are performed, and the end time of each timing truncation test is used ${\tau }_i\left({\tau }_1<{\tau }_2<\cdots <{\tau }_k\right)$ to indicate, there are $n_i$ samples that terminate the test at ${\tau }_i\left(i=1,2,\cdots ,k\right)$, and ${\displaystyle \sum _{i=1}^k{n}_i}=N$.

Let $r_i$ be the number of failed samples during $({\tau }_{i-1},{\tau }_i]$, and $t_{ij}\left(j=1,2,\cdots ,r_i\right)$ be the corresponding failure time. Here, $i=1,2,\cdots ,k$ and ${\tau }_0=0$. Thus $e_i={\displaystyle \sum _{q=1}^j{r}_q}$ indicates the total number of failed samples in the test. The samples are numbered from 1 to $N$ at each end of the timing truncation test. After all timing truncation tests are completed, the life number data are collected corresponding to the sample number $t_i\left(i=1,2,\cdots ,N\right)$. Apparently, $d={\displaystyle \sum _{i=1}^k{e}_i}=0$ reduces to the special case of zero-failure life test.

In addition, $s_i={\displaystyle \sum _{j=i}^k{n}_j}$, where $s_i$ consist of two parts, the first part is the number of samples still participating in the test at ${\tau }_i\left(i=1,2,\cdots ,k\right)$, the second part is failure samples $r_i$ at the internal $({\tau }_{i-1},{\tau }_i]$. Setting $p_i=F\left({\tau }_i\right)=P\left(t\le {\tau }_i\right)\left(i=1,2,\cdots ,k\right)$ as the failure probability at each censoring moment. It is obvious that $p_0=0$ at time t = 0, $p_1\le p_2\le \cdots \le p_k$.

Here, the distribution curve method is used to estimate the distribution parameters. The estimator of the failure probability $p_i$ at each truncation moment ${\tau }_i\left(i=1,2,\cdots ,k\right)$ is the precondition, and it is the core link used by the distribution curve method.

3. E-Bayesian Estimation of $p_i$ under Heavily Censored Data

The E-Bayesian parameter estimation method can be regarded as the modified Bayesian estimation method [28,29].

Definition 1:

Let the hyperparameter$b\in D$. If the Bayesian point estimate${\widehat{\theta }}_B\left(b\right)$is continuous, the E-Bayesian point estimation of the parameters$\theta$is

$${\widehat{\theta }}_E={\displaystyle {\int }_D{\widehat{\theta }}_B\left(b\right)\pi \left(b\right)}db$$

(1)

where${\displaystyle {\int }_D{\widehat{\theta }}_B\left(b\right)\pi \left(b\right)}db$exists,$D=\left\{b:1<b<c\right\}$,$c$is a constant,$\pi \left(b\right)$is the density function of the hyperparameter$b$on the interval$D$, and${\widehat{\theta }}_B\left(b\right)$is Bayesian estimation of$\theta$.

Equation (1) indicates that the E-Bayesian point estimation of $\theta$ is the expectation of the Bayesian estimate ${\widehat{\theta }}_B\left(b\right)$ about the hyperparameter $b$, where ${\widehat{\theta }}_E=E\left[{\widehat{\theta }}_B\left(b\right)\right]$.

3.1. Prior Distribution of $p_i$

In the case of heavily censored data, it is evident that the prior distribution density function of the failure probability $p_i$ is a decreasing function. Some common distribution functions in constructing the failure probability prior distribution are: Beta distribution, Gamma distribution, Exponential distribution, and Logarithmic distribution, where Beta (a, b) is a conjugate prior distribution of parameter $p_i$. The conjugate prior distribution can conveniently integrate historical information and provide a reasonable premise for future test data analysis. Some parameters can also be better explained simultaneously [30]. This study uses the Beta (a, b) distribution as the prior distribution of $p_i$, and adjusts the Beta (a, b) value range according to the change of the failure probability value range.

$$\pi \left(p_i|a,b\right)=\frac{p_i^{a-1}{\left(1-p_i\right)}^{b-1}}{B_i\left(a,b\right)}$$

(2)

where, $B_i\left(a,b\right)={\displaystyle {\int }_{l_i}^{u_i}{x}^{a-1}}{\left(1-x\right)}^{b-1}dx$

As this paper focuses on the case of heavily censored data in the timing data in the timing truncation test, the failure probability $p_i$ is more likely to be small value in [0, 1]. And the prior distribution density function $\pi \left(p_i|a,b\right)$ is derived:

$$\begin{array}{l}\frac{d\pi \left(p_i|a,b\right)}{d{p}_i}\\ =\frac{p_i^{a-2}{\left(1-p_i\right)}^{b-2}}{B_i\left(a,b\right)}\left[\left(a-1\right)\left(1-p_i\right)-p_i\left(b-1\right)\right]\le 0\end{array}$$

(3)

Equation (3) is established when $a\in (0,1]$ and $b\ge 1$. According to the property of the probability density of the Beta distribution function, when $a\in (0,1]$, the larger the value of $b$, the thinner the right tail of the probability density function. From the perspective of the robustness of the Bayesian estimation, the thinner the right tail of the prior distribution, the worse is the robustness of the Bayesian estimation. Therefore, there must be an upper limit on the value of $b$, which is $c$. In practical applications, $c$ is generally selected in an interval [2, 8]. It is very difficult to further determine the exact value of $a$ and $b$. Rather than using expert experience to determine a uniquely determined value, it is better to assume that $a$ and $b$ are uniformly distributed within their respective values, so the prior distribution density function of $a$ and $b$ is

$$\pi \left(a,b\right)=\frac{1}{c-1},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0<a<1,\hspace{0.17em}\hspace{0.17em}1<b<c$$

(4)

3.2. E-Bayesian Estimation of Failure Probability

In the model of heavily censored data, the likelihood function of the failure probability $p_i$ is

$$L\left(s_i,e_i|p_i\right)=C_{s_i}^{e_i}{p}_i^{e_i}{\left(1-p_i\right)}^{s_i-e_i}$$

(5)

Taking the prior distribution function of the failure probability $p_i$ in Equation (2) and the likelihood function in Equation (5) into the Bayesian formula, the poster distribution density function of failure probability $p_i$ can be obtained

$$\begin{array}{ll}h\left(p_i|r_i\right)& =\frac{\pi \left(p_i|a,b\right)L\left(r_i|p_i\right)}{{\displaystyle {\int }_{l_i}^{u_i}\pi \left(p_i|a,b\right)L\left(r_i|p_i\right)d{p}_i}}\\ & =\frac{p_i^{a+e_i-1}{\left(1-p_i\right)}^{b+s_i-e_i-1}}{{\displaystyle {\int }_{l_i}^{u_i}{p}_i^{a+e_i-1}{\left(1-p_i\right)}^{b+s_i-e_i-1}d{p}_i}}\end{array}$$

(6)

With the squared loss function, the Bayesian point estimation ${\widehat{p}}_i\left(a,b\right)$ of the failure probability $p_i$ is calculated as

$$\begin{array}{ll}{\widehat{p}}_i\left(a,b\right)& ={\displaystyle {\int }_{l_i}^{u_i}{p}_{i}h\left(p_i|r_i\right)}d{p}_i\\ & ={\displaystyle {\int }_{l_i}^{u_i}\frac{p_i{p}_i^{a+e_i-1}{\left(1-p_i\right)}^{b+s_i-e_i-1}}{{\displaystyle {\int }_{l_i}^{u_i}{p}_i^{a+e_i-1}{\left(1-p_i\right)}^{b+s_i-e_i-1}d{p}_i}}d{p}_i}\\ & =\frac{{\displaystyle {\int }_{l_i}^{u_i}{p}_i^{a+e_i}{\left(1-p_i\right)}^{b+s_i-e_i-1}d{p}_i}}{{\displaystyle {\int }_{l_i}^{u_i}{p}_i^{a+e_i-1}{\left(1-p_i\right)}^{b+s_i-e_i-1}d{p}_i}}\end{array}$$

(7)

when $l_i=0$ and $u_i=1$, Equation (6) becomes

$$\begin{array}{ll}h\left(p_i|r_i\right)& =\frac{\pi \left(p_i|a,b\right)L\left(r_i|p_i\right)}{{\displaystyle {\int }_0^1\pi \left(p_i|a,b\right)L\left(r_i|p_i\right)d{p}_i}}\\ & =\frac{p_i^{a+e_i-1}{\left(1-p_i\right)}^{b+s_i-e_i-1}}{B\left(a+e_i,b+s_i-e_i\right)}\end{array}$$

(8)

Equation (7) becomes

$$\begin{array}{ll}{\widehat{p}}_i\left(a,b\right)& ={\displaystyle {\int }_0^1{p}_{i}h\left(p_i|r_i\right)}d{p}_i\\ & =\frac{B(a+e_i+1,b+s_i-e_i)}{B(a+e_i,b+s_i-e_i)}=\frac{a+e_i}{a+b+s_i}\end{array}$$

(9)

Taking the prior distribution density function of $a$ and $b$ in Equation (4) into Equation (1), the E-Bayesian estimation of failure probability $p_i$ is

$$\begin{array}{ll}{\widehat{p}}_{iEB}& ={\displaystyle \int {\displaystyle {\int }_D{\widehat{p}}_i\left(a,b\right)\pi \left(a,b\right)}}dadb\\ & =\frac{1}{\left(c-1\right)}{\displaystyle {\int }_1^c{\displaystyle {\int }_0^1\frac{{\displaystyle {\int }_{l_i}^{u_i}{p}_i^{a+e_i}{\left(1-p_i\right)}^{b+s_i-e_i-1}d{p}_i}}{{\displaystyle {\int }_{l_i}^{u_i}{p}_i^{a+e_i-1}{\left(1-p_i\right)}^{b+s_i-e_i-1}d{p}_i}}dadb}}\end{array}$$

(10)

3.3. Value Range of Failure Probability

Motivated by reference [31], the distribution function of the life $T_i$ obeying the Weibull distribution is

$$F\left(t\right)=P\left(T\le t\right)=1-\exp \left[-{\left(\frac{t}{\eta }\right)}^m\right]$$

(11)

where $T_i~WE\left(m,\eta \right)$, $m$ is the shape parameter, and $\eta$ is the scale parameter. The corresponding reliability function is

$$R\left(t\right)=\exp \left[-{\left(\frac{t}{\eta }\right)}^m\right]$$

(12)

Regarding the determination of the failure probability $p_i$, the value range are $l_i$, and $u_i$, the current research is mainly divided into two methods. The first method takes $l_i=0$ and $u_i=1$. The second method uses the unevenness of the distribution function and the function characteristics to determine the value range of the failure probability. The method of setting the failure probability value range $l_i=0$ and $u_i=1$ is too conservative, and the characteristics of heavily censored data are not considered, so the calculation accuracy is poor. In order to improve the accuracy of the failure probability value range, this study uses method two to set the value range of the failure probability.

The second derivative of the Weibull distribution function is

$$\begin{array}{ll}\frac{d^{2}F\left(t\right)}{d{t}^2}& =\frac{m{t}^{m-2}\exp \left(-{\left(\frac{t}{\eta }\right)}^m\right)}{{\eta }^m}\left(m-1-\frac{m{t}^m}{{\eta }^m}\right)\\ & =\frac{m^2}{{\eta }^2}{\left(\frac{t}{\eta }\right)}^{m-2}\left(\frac{m-1}{m}-\frac{t^m}{{\eta }^m}\right)\exp \left[-{\left(\frac{t}{\eta }\right)}^m\right]\end{array}$$

(13)

when $m\le 1$, it is easy to prove that $\frac{d^{2}F\left(t\right)}{d{t}^2}<0$ by Equation (13), and thus the distribution function $F\left(t\right)$ is a convex function. According to the properties of the convex function, we can obtain

$$\frac{p_1}{t_1}>\frac{p_2}{t_2}>\cdots >\frac{p_{i-1}}{t_{i-1}}>\frac{p_i}{t_i}>\cdots >\frac{p_k}{t_k}$$

(14)

Because $p_1\le p_2\le \cdots \le p_k$, we can get the interval of $p_i$ as follows

$$p_{i-1}\le p_i<\frac{t_i}{t_{i-1}}{p}_{i-1}\begin{array}{c}\end{array}(i\ge 2)$$

(15)

The first step is to limit the value range of $p_1$ is $0<p_1<p_{\lim }$, where $p_{\lim }$ is the upper limit of $p_1$. The value of $p_{\lim }$ is based on the actual situation and the experience of an expert. For heavily censored data [13], $p_{\lim }$ is usually set to be less than 0.5.

When $m>1$, it is difficult to determine the concavity and convexity of distribution function $F\left(t\right)$. Transforming the Weibull distribution function according to the function characteristics, Equation (16) is obtained

$$\ln \ln \frac{1}{1-p}=m\ln t-m\ln \eta$$

(16)

Taking $p_{i-1}$ and $p_i$ into Equation (16), we have

$$\ln \frac{\ln (1-p_{i-1})}{\ln (1-p_i)}=m\ln \frac{t_{i-1}}{t_i}<\ln \frac{t_{i-1}}{t_i}$$

(17)

Then, Equation (17) is transformed into

$$\frac{\ln (1-p_{i-1})}{\ln (1-p_i)}<\frac{t_{i-1}}{t_i}$$

(18)

From Equation (18), $p_i>1-{\left(1-p_{i-1}\right)}^{\frac{t_i}{t_{i-1}}}$, it is known that $1-{\left(1-p_{i-1}\right)}^{\frac{t_i}{t_{i-1}}}$ is the lower limit of $p_i$ and $p_{upp}$ is the upper limit of $p_i$. The value of $p_{upp}$ is based on the actual situation. According to the unevenness and function characteristics of the distribution function, the value range of $p_i$ can be determined as follow

$$l_i=\{\begin{array}{c}{p}_{i-1}\\ 1-{\left(1-p_{i-1}\right)}^{\frac{t_i}{t_{i-1}}}\end{array}\begin{array}{c}m\le 1\\ m>1\end{array}\begin{array}{c}\end{array},u_i=\{\begin{array}{c}\frac{t_i}{t_{i-1}}{p}_{i-1}\\ p_{upp}\end{array}\begin{array}{c}m\le 1\\ m>1\end{array}$$

(19)

4. E-Bayesian Point Estimation of Reliability

In order to obtain distribution parameter estimates, points $\left({\tau }_i,p_i\right)$ first need to be obtained, and then a distribution curve based on the points needs to be assigned. As the exact value of the failure probability $p_i$ cannot be obtained during the calculation, the estimation value ${\widehat{p}}_{iEB}$ is used to describe the failure probability $p_i$. After all the ${\widehat{p}}_{iEB}$ are obtained, the distribution curve method is used to form a distribution curve based on these points $\left({\tau }_i,{\widehat{p}}_{iEB}\right)$. The current commonly used method for fitting a distribution curve is the weighted least square method. Reference [19] using the weighted least squares method to fit the distribution curve, the weight is $w_i=\frac{n_i{t}_i}{{\displaystyle \sum _{j=1}^k{n}_j{t}_j}}$. This method only considers the influence of the number of test samples and test time on the estimated value of the distribution parameter, and ignores the role of failed samples in the value of the test, so the method can only be used to estimate the distribution parameter of the sample without failure data. For test samples containing failure data, the method has obvious shortcomings. In order to consider the number of test samples, test time, and number of failed samples in the distribution curve, the original weighted least square method is improved. An improved least square method is proposed, the improved weight is $w_i=\frac{n_i{}^{r_i+1}{t}_i}{{\displaystyle \sum _{j=1}^k{n}_j{}^{r_i+1}{t}_j}}$, and the method is used to fit the distribution curve.

Let $y=\ln \left[-\ln \left(1-p\right)\right]$, $x=\ln t$, $b=m\ln \eta$, and then Equation (11) is transformed into

$$y=mx-b$$

(20)

With the improved weighted least squares method, the point estimates of distribution parameters $m$ and $\eta$ are obtained by fitting $\left(x_i,{\widehat{y}}_{iEB}\right)$, and then the point estimate of reliability is obtained. By minimizing the square of fitting errors, we have

$$Q={\displaystyle \sum _{i=1}^k{w}_i{\left({\widehat{y}}_{iEB}-m{x}_i+m\ln \eta \right)}^2}$$

(21)

where $w_i=\frac{n_i{}^{r_i+1}{t}_i}{{\displaystyle \sum _{j=1}^k{n}_j{}^{r_i+1}{t}_j}}$, $k\ge 2$.

For the derivative calculation of Equation (21), let $A={\displaystyle \sum _{i=1}^k{w}_i{x}_i}$, $B={\displaystyle \sum _{i=1}^k{w}_i{x}_i^2}$, $C={\displaystyle \sum _{i=1}^k{w}_i{\widehat{y}}_{iEB}}$, and $D={\displaystyle \sum _{i=1}^k{w}_i{x}_i{\widehat{y}}_{iEB}}$, giving

$$\hat{m}=\frac{D-AC}{B-A^2} \widehat{\eta }=\exp \left(\frac{AD-BC}{D-AC}\right)$$

(22)

Taking the estimated value of the distribution parameters obtained from Equation (22) into Equation (19), the point estimation ${\hat{R}}_{EB}\left(t\right)$ under the Weibull distribution can be obtained

$${\hat{R}}_{EB}\left(t\right)=R_{EB}\left(t;\hat{m},\widehat{\eta }\right)=\exp \left[-{\left(\frac{t}{\widehat{\eta }}\right)}^{\hat{m}}\right]$$

(23)

5. E-Bayesian Confidence Interval Estimation

Currently, all applications of the E-Bayesian estimation method focus on the point estimation of parameters, and no further research about parameter confidence interval calculations on this basis. For the estimation of parameter confidence intervals, the other methods should be used, such as the optimal confidence limit method [32] and bootstrap method [33]. Because the methods used for the parameter confidence interval estimation and the point estimation are different, the result will have “disjointed” problem. In fact, on the basis of the E-Bayesian estimation method and the distribution curve method, the confidence interval of the parameters can be obtained.

Under the confidence level value $\left(1-\alpha \right)$, the first step is to estimate upper confidence value ${\widehat{p}}_{iEB}^u$ of the failure probability $p_i$. Based on the poster distribution density function of $p_i$ in Equation (6) and the E-Bayesian estimate of the failure probability $p_i$ in Equation (10), $G\left(x,s_i,e_i\right)$ is defined as

$$G\left(x,s_i,e_i\right)=\frac{1}{\left(c-1\right)}{\displaystyle {\int }_1^c{\displaystyle {\int }_0^1\frac{{\displaystyle {\int }_{l_i}^x{p}_i^{a+e_i-1}{\left(1-p_i\right)}^{b+s_i-e_i-1}d{p}_i}}{{\displaystyle {\int }_{l_i}^{u_i}{p}_i^{a+e_i-1}{\left(1-p_i\right)}^{b+s_i-e_i-1}d{p}_i}}dadb}}$$

(24)

Similar to the estimation of upper confidence limit for the hierarchical Bayesian estimation method, Based on the definition of the upper confidence limit, we have

$$G\left({\widehat{p}}_{iEB}^u,s_i,e_i\right)=1-\alpha$$

(25)

Since $G\left(l_i,s_i,e_i\right)=0$ and $G\left(u_i,s_i,e_i\right)=1$, Equation (25) is solved by the binary search method, and then the upper confidence limit ${\widehat{p}}_{iEB}^u$ of failure probability $p_i$ is obtained. Based on the linearization process to Equation (18), ${\widehat{y}}_{iEB}^u=\ln \left[-\ln \left(1-{\widehat{p}}_{iEB}^u\right)\right]$ and $x_i=\ln t_i$ are obtained.

According to the improved weighted least square method, the failure probability upper confidence limit curve can be obtained by fitting points $\left(x_i,{\widehat{y}}_{iEB}^u\right)$. From Equation (21), we have

$$Q^u={\displaystyle \sum _{i=1}^k{w}_i{\left({\widehat{y}}_{iEB}^u-m{x}_i+m\ln \eta \right)}^2}$$

(26)

As $m$ has a limited value range and it has a small effect on the reliability index, so the confidence interval of $m$ is not assumed. $\hat{m}$ is the point estimation value of $m$, and then taking $\hat{m}$ into fitting the curve with $\left({\tau }_i,p_{iEB}^u\right)$, obtaining the point estimation value $\hat{m}$. The next step is to obtain the lower confidence limit ${\eta }_L$ of the estimation value $\eta$. Taking $\hat{m}$ into Equation (26), we have

$$Q^u={\displaystyle \sum _{i=1}^k{w}_i{\left({\widehat{y}}_{iEB}^u-\hat{m}{x}_i+\hat{m}\ln \eta \right)}^2}$$

(27)

The lower confidence limit ${\eta }_L$ of $\eta$ can be obtained when $Q^u$ have the smallest value:

$${\eta }_L=\exp \left[{\displaystyle \sum _{i=1}^k{w}_i\left(x_i-\frac{y_{iEB}^u}{\hat{m}}\right)}\right]$$

(28)

The point estimation value $\hat{m}$ and lower confidence limit ${\eta }_L$ are taken into the Weibull reliability calculation formula, and then the lower confidence limit of reliability $R$ at confidence level $\left(1-\alpha \right)$ can be obtained at any time $t$, it is

$$R_L\left(t\right)=\exp \left[-{\left(\frac{t}{{\eta }_L}\right)}^{\hat{m}}\right]$$

(29)

Similarly, the upper confidence limit of the failure probability is

$$F_u\left(t\right)=1-R_L\left(t\right)=1-\exp \left[-{\left(\frac{t}{{\eta }_L}\right)}^{\hat{m}}\right]$$

(30)

The upper confidence limit of the failure rate function is

$${\lambda }_u\left(t\right)=\frac{f_u\left(t\right)}{R_L\left(t\right)}=\frac{\hat{m}}{{\eta }_L}{\left(\frac{t}{{\eta }_L}\right)}^{\hat{m}-1}$$

(31)

6. Simulation Verification

In order to verify the accuracy of the proposed method, simulation experiments were used to generate simulation samples. Based on the generated simulation samples, the classic estimation method [34], the hierarchical Bayesian estimation method [19], and the E-Bayesian estimation method are used to evaluate the reliability of the sample. The results are compared from multiple aspects such as robustness, average time consumption and distance from the true value. Zhang [35] describes the process of generating simulation samples:

(1) Setting the true value of the shape parameter $m$ and scale parameter $\eta$ in the Weibull distribution, and then the failure time $t_i\left(i=1,2,\cdots ,N\right)$ is randomly generated through MATLAB, that is, there are $N$ experimental samples from the simulation experiment.

(2) $N$ uniformly distributed random numbers $u_i$, $u_i\in \left(0,1\right)$.

(3) The truncation level CL is determined, based on that, $r$ ($r=N\times CL$) failure times from all $t_i\left(i=1,2,\cdots ,N\right)$ are arbitrarily selected, denoted as $t_{fi}\left(i=1,2,\cdots ,r\right)$. The remaining failure time is represented by $t_{si}\left(i=1,2,\cdots ,N-r\right)$.

(4) $t_{sj}$ is replaced by $u_i{t}_{sj}$, which is used to indicate zero-failure truncation time $t_{sj}$.

(5) Based on the value of $t_{fi}\left(i=1,2,\cdots ,r\right)$ and $t_{sj}$, the truncation time is grouped and the minimum truncation time is taken as the truncation time of the censored experiment of the group.

For this simulation, 11 samples were selected for 5 censored experiments. In the Weibull distribution, the shape parameter $m=2.2$ and scale parameter $\eta =2400$. Taking the truncation level CL = 0.1, the simulation samples were generated, they were shown in

Table 1. Simulation samples.

$\mathbf{Number}\mathit{i}$	1	2	3	4	5
Censored time $t_i$	682	813.5	1460	2090.7	2390.7
Samples $n_i$	2	1	1	2	5
Failed samples $r_i$	0	0	0	1	0
Summary of failed samples $e_i$	0	0	0	1	1
Remaining samples $s_i$	11	9	8	7	5

.

From Equation (14), when $m>1$, the upper limit $p_{upp}$ of $p_i$ needs to be determined, since even in a censored test, the value of $p_{upp}$ will continuously change with time $t$. It will increase the calculation complexity and cumulative error, so generally the value of $p_{upp}$ is set to 1 and the value range of failure probability $p_i$ is $\left[1-{\left(1-p_{i-1}\right)}^{\frac{t_i}{t_{i-1}}},1\right]$.

Based on the data in Table 1, the classical estimation method [34], the hierarchical Bayesian estimation method [19], and the E-Bayesian estimation method proposed in this study were used to solve the failure probability estimates. The calculation results of the failure probability at time $t_5$ are shown in

Table 2. The failure probability $p_i$ at c.

$\mathit{c}$	2	3	4	5	6	7	8	Range	Average Running Time/s
Method proposed in this study	0.5615	0.5428	0.5258	0.5103	0.4961	0.4829	0.4706	0.0909	102
Method proposed in [19]	0.2221	0.2092	0.1981	0.1887	0.1805	0.1734	0.1672	0.0549	1038
Method proposed in [34]	0.0833							/	$\approx 0$
True value	0.6290							/	/

.

Table 2 gives the failure probability estimation results of three different methods. The estimation range under different $c$ values and average running time are also compared in the table. Compared with [19], although the robustness of the proposed method is slightly lower, and the algorithm efficiency is significantly improved. Furthermore, compared with [19,34], our method has higher calculation accuracy due to modify the value range of failure probability.

$c=5$ is used in subsequent calculation, and reliability ${\hat{R}}_{EB}\left(t\right)$ calculated by the method proposed in this study is compared with ${\hat{R}}_{HB}\left(t\right)$ and ${\hat{R}}_B\left(t\right)$ which are calculated by methods proposed in [19,34] respectively, the results are shown in Figure 1.

In the Figure 1, Method 1 is the hierarchical Bayesian estimation method [19], Method 2 is the classic estimation method [34], and Method 3 is the E-Bayesian reliability evaluation method proposed in this study. From Figure 1, it can be seen that E-Bayesian reliability assessment methods with heavily censored data have the closest estimated results to the true values. The classic estimation method [34] does not consider the effect of censored time on the failure probability in the timing truncation test, which makes the reliability overestimated. In the process of determining the failure probability value range, the hierarchical Bayesian estimation method [19] fixed the interval to [0, 1], without considering the sequence relationship between the failure probability values. The determination of the failure probability value range is too conservative. In addition, in the process of solving the distribution parameters using the weighted least squares method, the weight model does not consider the effect of the number of failed samples on the weight value. Due to the above reasons, the obtained reliability estimation value is too large.

In order to solve the lower confidence limit of reliability, the first step is to calculate the upper confidence limit $p_i^u$ of failure probability $p_i$ under the confidence level of 0.9, and then getting the lower confidence limit of the reliability parameter. As reference [34] does not further analyze the confidence interval of the reliability parameter, the method proposed in this study is only compared with the method in [19]. The corresponding lower confidence limit (LFM) comparison results are shown in Figure 2.

In Figure 2, Method 1 is the hierarchical Bayesian estimation method [19] and Method 3 is the E-Bayesian reliability evaluation method proposed in this study. From the comparison results, it can be seen that the method proposed in this study is better than the hierarchical Bayesian method proposed in [19].

7. Case Study

7.1. Transmission Mechanisms

Most of the transmission mechanisms in the integrated transmission are planetary transmission mechanisms. This type of planetary transmission mechanism refers to a transmission mechanism in which the structure includes a planetary gear that rotates around a central axis. Planetary gears are internally connected parts, which are connected to the various elements that rotate around the central axis, including the sun gear, ring gear, and planet carrier, forming various planetary differentials (planetary rows). It can be used alone as a differential, reducer, or power splitting and converging mechanism in the drive train, or it can be used in series to form a speed increasing and reducing mechanism in the drive train. The flexible application of planetary mechanisms has changed a lot and has become the basis for many new developments in transmissions. Among them, the sun gear is a key component of the planetary transmission mechanism, which affects whether the planetary transmission mechanism can be normally used.

In this study, the sun gear of a transmission mechanism is used as the research object, and the timing truncation test is carried out under its normal working conditions. The test data is shown in

Table 3. Sun Gear Test Data.

$\mathbf{Number}\mathit{i}$	1	2	3	4	5
Censored time $t_i$	1000	1500	2000	2500	3000
Samples $n_i$	3	3	3	3	3
Failed samples $r_i$	0	0	0	1	0
Summary of failed samples $e_i$	0	0	0	1	1
Remaining samples $s_i$	15	12	9	6	3

.

Where $c=5$, the failure probability E-Bayesian point estimate value ${\widehat{p}}_{iEB}$ at every censored moment $t_i$ is calculated, and the results are shown in

Table 4. Failure probability and upper confidence limit point estimation.

$\mathbf{Number}\mathit{i}$	1	2	3	4	5
Failure probability estimated value ${\widehat{p}}_{iEB}$	0.02689	0.08861	0.17639	0.31400	0.48287
The upper confidence limit of failure probability ${\widehat{p}}_{iEB}^u$	0.0474	0.1193	0.2143	0.3767	0.5604

.

Then, the point estimation (solid line in Figure 3) curve of the reliability under the Weibull distribution is obtained.

In this study, we set the confidence level $\left(1-\alpha \right)=0.8$. With Equation (25), the upper confidence point estimation value ${\widehat{p}}_{iEB}^u$ (dash line in Figure 3) of failure probability $p_i$ is obtained. The results are shown in Table 4.

The lower confidence limit can be obtained at any $t$ moment when the reliability $R$ at the confidence level 0.8, and the results are shown in Figure 3.

With Equation (31), the upper confidence limit of the failure rate function ${\lambda }_u\left(t\right)$ is obtained, and the result is shown in Figure 4.

7.2. GPS Receiver

The measurement and control subsystem, it is an important subsystem of the satellite platform. It cooperates with the ground measurement and control system to realize the speed, distance and angle measurement of a satellite. In addition, this subsystem is also used to collect relevant parameters inside the satellite. GPS receiver is one of the components of measurement and control subsystem.

We collected 10 times telemetry data, and there are no failure data of all GPS receivers. The receiving times are shown in Figure 5.

Firstly, sort the samples in descending order, namely $t_{10}>\cdots >t_1$,

With the GPS receiver data, we can obtain the reliability point estimation and remaining samples $s_i$. Set the confidence level to be $\left(1-\alpha \right)=0.9$. Finally, the reliability point estimation and lower confidence limit of GPS receiver for 1 to 8 years are calculated in Figure 6.

8. Conclusions

In this paper, an improved E-Bayesian estimation method based on Weibull distribution is proposed. Compared with the traditional ones, the main contribution is twofold: (1) We extend the application range of E-Bayesian estimation from zero failure data to extreme few failure data; (2) A novel confidence interval estimation approach for E-Bayesian is provided.

Simulation results show that although the proposed method sacrifices few robustness, it outperforms the traditional ones in the following aspects: (1) Compared with the hierarchical Bayesian estimation method, the E-Bayesian estimation method proposed in this study has simpler mathematical models and calculations. The efficiency is significantly improved. (2) The value range of the failure probability is improved by using the unevenness and function characteristics of the distribution function, and the improved weighted least square method is used to solve the distribution parameters. Compared with the hierarchical Bayesian estimation method, the E-Bayesian estimation method proposed in this study has higher calculation accuracy.

Finally, the method proposed in this study is applied to the data analysis of the sun gear timing truncation test of a transmission mechanism and GPS receiver. It further verified that the method proposed in the study is more conducive to practical application in engineering. The results show that the heavily censored data based on E-Bayesian estimation method can solve the reliability evaluation problem under the condition of non-failure data and few failure data. Therefore, the method proposed in this paper has a wider range of engineering applications.

Reliability assessment methods for heavily censored data has not been studied thoroughly and there are still many issues to be resolved. In practice, we can only have very few failure data of products, it is difficult to collect product performance data. In addition, due to the limited test time and samples, it is generally impossible to collect sufficient life data. In the future, we will pay more attention to the combination of life data and other reliability data, as well as the utilization of historical data.