Statistical Analysis of a Lifetime Distribution with a Bathtub-Shaped Failure Rate Function under Adaptive Progressive Type-II Censoring

In this paper, the estimation problem of two parameters of a lifetime distribution with a bathtub-shaped failure rate function based on adaptive progressive type-II censored data is discussed. This censoring scheme allows the experiment to be more efficient in the use of time and cost while ensuring the statistical inference efficiency based on the experimental results. Maximum likelihood estimators are proposed and the approximate confidence intervals for two parameters are computed using the asymptotic normality. Lindley approximation and Gibbs sampling are used to obtain Bayes point estimates and the latter is also used to generate Bayes two-sided credible intervals. Finally, the performance of various estimation methods is evaluated through simulation experiments, and the proposed estimation method is illustrated through the analysis of a real data set.


Chen(β, θ) Distribution
In practical analysis, the failure rate function is the key for researchers to study the life phenomenon of a product. In many industries, products are more likely to have a bathtub-shaped failure rate curve than a fixed failure rate. For the bathtub-shaped failure rate function, the product failure rate begins with a decreasing process, then a constant phase, and finally an ascending phase. This model has wide application prospects. For example, some treatments or drugs may work better in middle-aged men than in other age groups. That is, the failure rate of the treatment for children is high, and it will continue to fall to a constant level, then it will maintain that constant level until a certain age, later will rise again as the growth of the age gradually. The exponentiated Weibull extension proposed by Chen [1] (written as Chen(β, θ) distribution and abbreviated as the Chen distribution throughout the paper) is a lifetime distribution with bathtub-shaped failure rate function. The confidence regions of the two unknown parameters of the Chen(β, θ) distribution have closed forms, which are important for practical applications. Chen [1] proposed an exact confidence interval for two parameters based on type-II censored data. Many researchers have also been interested in this new lifetime distribution and have conducted a series of studies. Based on the same censored data, Wang [2] explored some analytical properties of this distribution, added moment estimation, and evaluated the performance of the proposed estimation methods through intensive Monte Carlo simulations. Sanku Dey [3] focused on frequentist analyses and Mohammad Z Raqab [4] further studied bayesian analyses with the record data. Under the progressive type-II censored data, Wu [5] proposed maximum likelihood estimation and a method to construct exact confidence intervals as well as exact joint confidence regions of parameters. Manoj Kumar Rastogi [6] further discussed the estimation of its reliability function as well as hazard function based on the same censoring scheme. Bahman Tarvirdizade [7] considered the estimation of the stress-strength reliability based on upper record values. Essam Kh Ahmed [8] further used the Markov Chain Monte Carlo method (MCMC) for Bayes estimation with the progressive type-II censored data. Zhang [9] discussed acceleration factors based on the tampered failure rate model.
The probability density function (pdf) of Chen(β, θ) is given by and cumulative distribution function (cdf) is where β > 0, θ > 0 are the parameters. Some of its characteristics are illustrated below. f (x) is decreasing for 0 < θ ≤ 1, β ≤ 1; f (x) is unimodal for β > 1. When θ > 1, f (x) is either decreasing, or decreasing first and then increasing and decreasing finally (see [2] in details). Figures 1-4 illustrate the curves of f (x) under several settings of two parameter values.
The failure rate function of the Chen(β, θ) is Taking the derivation, λ (x) = β θ · x β−2 · e x β · (β − 1) + βx β , from which we can conclude that λ(x) is in a bathtub-shaped curve when 0 < β < 1 and increases when β ≥ 1. Figures 5 and 6 show the curves of λ(x) under various settings of β and θ. Since the failure rate at the bottom of the bathtub is stable, which is closely related to the working stage of the actual product, it is of great significance to analyze this kind of lifetime distribution for the practical application in industry.

Adaptive Progressive Type-II Censored Scheme
In the life test and reliability test, the experimenters use various censoring schemes to reduce the total test time and related costs among which type-II censoring scheme is one of the most widely used. Under type-II censoring scheme, only the first m failed units in a random sample in size n (m < n) are observed. This method is easy to be carried out by people, while it has the possibility of wasting a lot of testing time. Therefore, the progressive type-II censoring as a generalized form of type-II censoring is proposed to improve the efficiency of the experiment. The progressive type-II censoring can be implemented as follows. Let n units be placed on a life test and denote X 1 , X 2 , · · · , X n as their lifetimes, respectively. It is suppose that X i , i = 1, 2, . . . , n are independent to each other and all follow pdf f (x) and cdf F(x). Before the experiment begins, the number of units observed m (m < n) as well as the progressive type-II censoring scheme R = (R 1 , R 2 , . . . , R m ) are determined where R i > 0, i = 1, 2, · · · , m and ∑ m i=1 R i + m = n. When i-th failure is observed, the surviving R i units are randomly removed from the experiment. The experiment continues according to this rule until m failures are observed and the experiment is over.
Denote by m the completely observed lifetimes by X R i:m:n , i = 1, 2, · · · , m (simplified to X i:m:n in a specific experiment), which are the observed progressive type-II right censoring order statistics.
A basic assumption of the design of this experiment is that the censoring scheme R = (R 1 , R 2 , · · · , R m ) is set before the experiment begins . This model can simulate the real-life situation where the loss  or removal of some units happens during the experiment, which makes it more reasonable than  type-II censoring. However, this method may still fail to be satisfied in reality since experimenters may want to change the censoring numbers according to the actual situation during the experiment process. Thus, adaptive progressive type-II censoring is proposed [10] to improve the flexibility and efficiency of the experiment. In the adaptive progressive type-II censoring scheme, the size of units m observed is decided before the experiment and the censoring scheme R = (R 1 , R 2 , · · · , R m ) is roughly decided, but R i may be modified according to certain needs of experimenters. The working process is demonstrated as follows: suppose there is a pre-provided expected finishing time T of the test, though we allow the total testing time to exceed T. If X m:m:n < T, the experiment is carried out in the same way as type-II progressive censoring. Otherwise, once the experiment time pass time T while we have not observed m failures, we leave as many surviving units as possible, hoping to see more failures in a short time to finish the experiment in the most efficient way. Therefore, suppose X J:m:n is the last failure observed before time T (i.e., X J:m:n < T < X J+1:m:n ) where X 0:m:n ≡ 0 and X m+1:m:n ≡ ∞. After time T, we do not drop any units until the m-th failure is observed. Then we withdraw all surviving units. That is to say, Figures 7 and 8 illustrate progressive type-II censoring and adaptive progressive type-II censoring, respectively.  Based on adaptive type-II progressively censored data, the maximum likelihood estimation and Bayes estimation for two parameters of Chen(β, θ) distribution have been studied. The rest of this paper is organized as follows. Maximum likelihood point estimation is considered and asymptotic confidence intervals using the asymptotic normality are constructed in Section 2. Then in Section 3, we discuss Bayes estimation and obtain numerical results to Bayes point estimators using Lindley approximation and Gibbs sampling. Furthermore, the latter is used to construct Bayes two-sided credible intervals of unknown parameters. In Section 4, the R program is used for simulation study to compare different estimation methods. In Section 5, a real data set is analyzed to demonstrate the feasibility of the proposed estimation methods. Finally, the thesis is summarized in the Section 6.

Maximum Likelihood Estimation
Suppose that X = X R 1:m:n , X R 2:m:n , . . . , X R m:m:n is an adaptive progressive type-II censored sample of size m from a sample of size n with censoring scheme R = (R 1 , R 2 , · · · , R m ) taken from distribution having f (x) as the pdf and F(x) as the cdf, and X J:m:n is the last observed failure before T which is prefixed best testing time. The parameter x = x R 1:m:n , x R 2:m:n , . . . , x R m:m:n (simplified as x = x 1 , x 2 , · · · , x m in later equations) is used to represent the observed values of such an adaptive type-II progressively censored sample. On this basis, the corresponding likelihood function is given by where Therefore, the likelihood function for X R 1:m:n , X R 2:m:n , ..., X R m:m:n taken from the Chen(β, θ) is correspondingly written as Further, the log-likelihood function can then be written as Then take the partial derivatives of the log-likelihood function, and obtain the likelihood equations as follows: From (8), the maximum likelihood estimate of θ can be illustrated as a function ofβ Thus, givenβ, the correspondingθ with (9) can be obtained. However, we cannot get the closed forms of (7) and (9). Therefore, numerical methods are included, such as Newton's method, to determine the expected estimates.
Further, asymptotic confidence intervals are discussed. When the sample size n is large enough, the Fisher information matrix can be involved to get approximate confidence intervals of two unknown parameters. The Fisher information matrix is derived based on likelihood functions. By formula derivation, the pivot quantities with unknown parameters subject to asymptotically normal distribution are constructed. Then the asymptotic confidence intervals are obtained through Monte Carlo simulations.
The Fisher information matrix I(β, θ) is written as and from the log-likelihood function (6), we have the second partial of l(β, θ | x) with respect to β and θ which is illustrated in the Appendix A. It is obvious that it is difficult to get the exact expression of the above expectation. Therefore, the observed Fisher information matrix without taking expectation will be used instead. Let (β,θ) represent the maximum likelihood estimate of (β, θ), the observed Fisher information matrix I(β,θ) is given as follows: Correspondingly, an approximated covariance matrix of (β,θ) is given as follows: The asymptotic normality of maximum likelihood estimation (MLE) supports that (β,θ) follow bivariate normal distribution with mean vector (β, θ) and covariance matrix I −1 (β,θ) (12) approximately. Thus, (1 − α) × 100% approximate confidence intervals for the two parameters β and θ can be written respectively, as where z α 2 is the percentile of the standard normal distribution N(0, 1) with right-tail probability α 2 .

Bayes Estimation
Different from maximum likelihood estimation, Bayes estimation analyzes the problem by combining the prior information of random samples from a certain distribution and maximum likelihood estimation. Therefore, Bayes estimation considers both the provided data and the prior probability and can infer the interested parameters more reasonably. When there is no prior probability information, non-information priors are considered.
First, it is supposed here that the unknown parameters β and θ are independent and follow the gamma prior distributions Then, the joint prior distribution for β and θ becomes According to the Bayes theory, we subsequently get the joint posterior probability distribution π(β, θ | x) as: . Thus, the marginal posterior probability distributions of β and θ can be expressed, respectively, as Generally, squared error loss function λ(ξ, ξ) of unknown parameter vector ξ = (ξ 1 , ξ 2 , · · · , ξ n ) is as follows: whereξ = (ξ 1 ,ξ 2 , · · · ,ξ n ) represents the estimator of ξ = (ξ 1 , ξ 2 , · · · , ξ n ). Then we consider the Bayes estimators of β and θ under squared error loss function and their corresponding minimum Bayes risks. It is proved that they are the expectation and the variance of the function (17) and (18), respectively. We can easily obtain Bayes estimates for the unknown parameters β and θ and the corresponding minimum Bayes risks by using the marginal posterior distribution function to generate the first and second moments of β and θ. Using µ (r) β and µ (r) θ (r = 1, 2, · · · ) to represent the r-th central moments of β and θ, then we have equations below Thus, the Bayes point estimators and the corresponding minimum Bayes risks (MBR) of β and Furthermore, since the marginal posterior probability distributions of β and θ are obtained, the Bayesian two-sided credible intervals can also be derived.
Similarly for θ, we can construct a credible interval, say (θ L , θ U ), using integral equations below: Obviously, the Bayes estimation formulas above do not have closed forms in general, so we will use approximation techniques and numerical methods to generate the Bayes estimates. The Lindley approximation and Gibbs sampling are illustrated in the following sections.
Then the ratio of integral (27) can have a linear approximation as follows: and σ ij is the (i, j)-th element of the Fisher information matrix.
In this case, ξ = (β, θ) and σ ij is the (i, j)-th element of the matrix (10). The Lindley approximation to Bayes estimators of β and θ (i.e., (20), (21))under the squared error loss function arẽ whereβ MLE andθ MLE are maximum likelihood estimates of β and θ. From (6), other expressions in (29) and (30) can be calculated which are shown in details in the Appendix B.

Gibbs Sampling
In statistics, Gibbs sampling is one of Markov chain Monte Carlo (MCMC) algorithm for approximating a set of observations from an interested multivariate probability distribution in situations where direct sampling is difficult. In this section, Gibbs sampling is involved to calculate the expectation (20) and (21) and two-sided credible intervals (22)-(25). Although it is difficult to generate samples directly from a posterior distribution, combining Metropolis-Hasting algorithm [12] with Gibbs sampling can easily generate samples. Two-sided credible intervals can be generated with this method. Several Bayesian estimation methods of two parameters of the Chen distribution have already been discussed in many references, like in [6], but the Gibbs (M-H) sampling method has not been involved yet. Therefore, in this section, we will conduct Bayes estimates through Gibbs sampling.
We can simplify the posterior distribution of parameters β, θ to π * (β, θ | x) as the Gibbs sampling density function, because the asymptotic estimate of (β, θ) through the Gibbs (M-H) sampling method can be obtained without calculating the constants.
represent the conditional density functions of the parameters, respectively, as follows: Then, sample and get a series values (β 1 , θ 1 ), (β 2 , θ 2 ), · · · , (β τ , θ τ ) and finally obtain Bayesian point estimation and 100(1 − α)% two-sided credible intervals from those new samplings. The main idea of Gibbs sampling is that instead of getting the estimate of (β, θ) simultaneously, we estimate the parameter discretely in a sequence, which can significantly simplify the process. Obviously, the conditional density functions of β and θ given in (32) and (33) cannot transfer analytically into distributions that can be sampled directly using existing algorithms. Therefore the Metropolis-Hastings algorithm with a normal proposal distribution is involved to generate β and θ. The Gibbs (M-H) sampling procedure can be explained briefly in the following steps.
Thus, a series of Gibbs (M-H) samplings of λ = (β, θ), (β 1 , θ 1 ), (β 2 , θ 2 ), . . . , (β τ , θ τ ) is obtained through former process. Now, the approximate expectation of λ = (β, θ) could be calculated as, where bt represents the burn-in time of the Gibbs sampling process when the samples can be assumed to come from the posterior distribution (16), which is a constant determined in advance. Under the squared error loss function, the Bayesian point estimate of parameter λ = (β, θ) is generated with the procedure mentioned above.

Simulation and Comparisons
In this section, the simulation study is carried out to verify the feasibility of our estimation method and compare the effects of different estimation methods.
Using the algorithm offered by [10,13], the adaptive progressive type-II censored data from Chen(β, θ) is generated with R program. Maximum likelihood estimates are computed by solving (7) and (8) in R program. The informative hyper-parameters (a 1 , a 2 , b 1 , b 2 ) in (15) are set to be (0, 1, 0, 1) and Bayes estimates under Gibbs sampling are obtained using (34) and (35). Bayes point estimates under the Lindley approximation method are obtained by calculating (29) and (30). The simulation experiment is conducted for N = 1000 times, and the mean values are taken as expected values (EV) of the unknown parameters, as it is illustrated below.
In order to compare the performances of maximum likelihood estimators for different settings of parameters, mean square error (MSE) of the estimated values φ (β, θ) is also calculated as, The asymptotic 95% confidence interval obtained by fisher observation matrix is obtained by (13) an (14). The Bayes two-sided 95% credible interval using Gibbs sampling is generated with Algorithm 1. The simulation experiment is repeated for N = 1000 times, and the corresponding interval coverages and average interval lengths are calculated using R. We further present the average 95% intervals of the parameter setting (0.5, 7) to discuss the reason for the difference in interval coverage and average interval length of interval estimation results obtained by the two methods.
(1) Tables 2 and 3 present the comparison of the maximum likelihood method and Bayes (using Gibbs-MH) method in terms of point estimation results. For both maximum likelihood estimation and Bayes estimation, under the condition that the parameters to be estimated remain unchanged when the sample size n increases, the expected values are closer to the set true values and the mean square errors are smaller. Simulation results show a better fitness performance when the expected total time on test is set with T = 2 under most censoring schemes. Although the results in simulation experiments of other parameter settings are satisfactory, the estimation effect of (0.5, 7) is poor. When the sample sizes n and expected finishing time T are fixed, the first censoring scheme (n − m, 0, · · · , 0) is more efficient than others. Bayes estimates are, in most cases, slightly more effective than maximum likelihood estimates in terms of expected values and smaller mean square errors.
(2) Tables 4 and 5 present the coverage probability and average length of two-sided 95% confidence intervals of parameters constructed by maximum likelihood method (using fisher observation matrix) and the Bayes method (using Gibbs (M-H)). For the same set of parameter, both the asymptotic confidence interval obtained by fisher observation matrix and two-sided probability interval constructed with Gibbs sampling are more accurate when the sample size n increases. That is, the interval coverage probability of the 95% asymptotic confidence interval of β and θ will be closer to 95%; and the average length of the interval is getting shorter. Meanwhile, other variables in the simulation experiment, namely the experiment end time T and different censoring schemes R, had no significant effect on the results of interval estimation. For the asymptotic confidence intervals using fisher observation matrix, under almost all parameter settings, the asymptotic confidence intervals of β are more ideal than that of θ, that is, the coverage of 95% confidence intervals of β are closer to 95%, while the coverage probability of θ are always low. In this term, the two-sided credible intervals of the Gibbs sampling method perform comparatively better.   (3) The coverage probability and average length of the interval obtained by two methods are significantly different under the set of parameters (0.5, 7). For the estimation results of β, the interval coverage probability of the asymptotic confidence interval obtained by the Fisher observation matrix is slightly higher than that obtained by Gibbs sampling, while the average interval length of the latter is slightly lower. For the estimation results of θ, the interval coverage probability of the asymptotic confidence interval obtained by the Fisher observation matrix is noticeably higher than that obtained by Gibbs sampling, but the average interval length of the latter is only half of that of the former. Therefore the average interval results of (0.5, 7) are shown in Table 6 for further illustration. There are obviously different results for the upper bound of the interval of θ under the two methods. For the asymptotic confidence interval obtained by the Fisher observation matrix (where the average length is longer), the truth value of θ = 7 is in the middle of the interval and can be well covered. For the interval obtained by Gibbs sampling, the truth value of θ = 7 is very close to the upper bound of the credible interval, so the randomness of the simulation experiment may lead to the fact that the true value is not covered in the credible interval.
(4) Lindley approximation method, as a kind of linear approximation method in the Bayesian estimation, whose estimation result is not ideal from the perspective of estimated value and mean square error when the lifetime distribution property is special and the sample structure is complex. Therefore, we only show Lindley approximation estimation results of (0.5, 7) when T = 1 in Table 7 to illustrate. In addition, Lindley approximation can only obtain point estimation, but not confidence interval which is more widely used in real life, so compared with it, maximum likelihood estimation and Bayes estimation obtained by Gibbs-MH sampling are better. Table 4. Coverage probability (CP) and average length (AL) of 95% intervals when T = 1.

Real Data Analysis
In this section, we apply our estimation methods in a real data set to verify the feasibility. The data set from Michele D. Nichols and W. J. Padgett [14] was used for illustration. A complete sample of data shows the observed fracture stress of 100 carbon fibers (shown in Table 8).
The data was previously thought to come from the exponentiated Weibull distribution with probability density function of αθx α−1 · exp(−x α ) · (1 − exp(−x α )) θ−1 (EW distribution) in [15]. We used the Kolmogorov-Smirnov (K-S) teat and compared the K-S distances (maximum vertical deviation between the probability density curve of the dataset and the given probability density curve) and p-value (in statistical hypothesis test) of four multi-parameter distributions, including the Chen distribution and the EW distribution. The results in Table 9 show that the K-S distance of the Chen distribution is smaller than that of the EW distribution, indicating that the fitting effect of the data set is better, so it is reasonable to assume that the data set comes from the Chen distribution.  Table 9. The results of Kolmogorov-Smirnov (K-S) test for alternative statistical models.

Statistical Models Probability Density Function K-S Distance P Value
Generalized Gompertz [16] θβ · e λx e − β λ (e λx −1) 1 − e − β λ (e λx −1) Then consider the estimation problem under adaptive type-II progressively censored data. We withdraw samples of size m = 60 from the complete sample of size n = 100. Censoring schemes are set to be R = (20, 0 * 58, 20) and R = (1 * 20, 0 * 20, 1 * 20), while two values of expected finish time T are (1.6,3.2). Tables 10-13 illustrate our adaptive type-II progressively censored data in details.  We consider gamma priors for Bayes estimation on both β and θ where the hyper-parameters are 0 (i.e., a 1 = a 2 = b 1 = b 2 = 0). Using the formulas mentioned in the simulation experiment, we can obtain the maximum likelihood point estimates, confidence intervals using the observed fisher matrix, Bayes point estimates, and Bayes credible intervals of unknown parameters β and θ of the real data set. The estimation results are shown in Table 14. We have conclusions as follows: (1) With the same estimation method (i.e., MLE, Gibbs-MH, Lindley), results show little difference under different censoring schemes R and best finishing times T.
(2) It can be seen from the estimated value of the parameters (β, θ) that the Chen distribution of the data follows has a bathtub-shaped failure rate function. There is no significant difference between the results of point estimation and interval estimation obtained by MLE and Gibbs sampling. The point estimation obtained by the Lindley approximation method is different from that of the former two methods. The reason can be that the setting of parameter (β, θ) allows the failure rate function of the Chen distribution takes the shape of a bathtub, which makes it more complicated to estimate its parameters. This difficulty is also illustrated in our simulation process. We find that the MSE of the point estimates, as well as interval coverage probability of the confidence interval for the setting (0.5, 7), are apparently larger than other parameter settings.

Conclusions
In this paper, the estimation problem of two unknown parameters β and θ from the Chen distribution is discussed. This distribution has a bathtub-shaped failure rate function, making it suitable for fitting many real-world data. The statistical properties of this highly applicable distribution are studied based on adaptive type-II progressively censored data, which make our research easier to be applied to practical industrial fields. This censoring scheme gives people better control over the process of experiments in terms of time and funding because it can shorten or lengthen the process.
We discuss the maximum likelihood estimators as well as asymptotic confidence intervals by including the observed Fisher information matrix for unknown parameters. On the premise that the prior distributions are gamma distributions, the theoretical results of the Bayes point estimation and Bayes two-sided credible intervals under the squared error loss function are also given. Finding that accurate estimates are not available, we use approximation techniques such as Lindley approximation and Gibbs sampling to calculate the results.
In the simulation and real data tests, the maximum likelihood estimates and the Bayes estimate under Gibbs sampling show an insignificant difference, though the latter has slightly better performances than the former. In terms of estimated value and mean square error, the result of the point estimator obtained by Lindley approximation is not as good as the former two. When the parameters are set so that the Chen distribution has a bathtub-shaped failure rate function, the accuracy of the estimation results decreases. In the future work, we may explore the optimization algorithm for this remaining problem and propose a more accurate estimation method.
Moreover, we will also study the change-point analysis of the lifetime distribution having a bathtub-shaped failure rate function, which has a wide range of applications. For example, it can help people better analyze the performance of products and predict the total life of products in the actual production.