Accelerated Life Test Method for the Doubly Truncated Burr Type XII Distribution

The Burr type XII (BurrXII) distribution is very flexible for modeling and has earned much attention in the past few decades. In this study, the maximum likelihood estimation method and two Bayesian estimation procedures are investigated based on constant-stress accelerated life test (ALT) samples, which are obtained from the doubly truncated three-parameter BurrXII distribution. Because computational difficulty occurs for maximum likelihood estimation method, two Bayesian procedures are suggested to estimate model parameters and lifetime quantiles under the normal use condition. A Markov Chain Monte Carlo approach using the Metropolis–Hastings algorithm via Gibbs sampling is built to obtain Bayes estimators of the model parameters and to construct credible intervals. The proposed Bayesian estimation procedures are simple for practical use, and the obtained Bayes estimates are reliable for evaluating the reliability of lifetime products based on ALT samples. Monte Carlo simulations were conducted to evaluate the performance of these two Bayesian estimation procedures. Simulation results show that the second Bayesian estimation procedure outperforms the first Bayesian estimation procedure in terms of bias and mean squared error when users do not have sufficient knowledge to set up hyperparameters in the prior distributions. Finally, a numerical example about oil-well pumps is used for illustration.


Historical Review and Literature Review
Dr. I. W. Burr was the pioneer to propose Burr type distributions in 1942 [1]. Since then, the two-parameter Burr type XII (BurrXII) distribution has earned more attention and has been widely used for reliability inferences because of the flexibility for modeling via using its two-shape parameters. Tadikamalla [2] extended the two-parameter BurrXII distribution to be a three-parameter BurrXII distribution by introducing one scale parameter. The three-parameter BurrXII distribution is a generalized version of the BurrXII distribution and includes the gamma, bell-shaped, lognormal, and log-logistic distributions as special cases. The three-parameter BurrXII distribution is also an asymptotic limiting case of the Weibull and Pareto type I distributions.
Al-Hussaini [3] extended the theorem established by Galambos and Kotz [4] to characterize the two-parameter BurrXII distribution. Zimmer et al. [5] studied the reliability applications for the two-parameter BurrXII distribution. Jang et al. [6] used Bayesian estimation method distribution to model heavy tailed lifetime data and proposed parameter estimation methods for the doubly truncated three-parameter BurrXII distribution. Kantar and Usta [28] used the upper-truncated version of Weibull distribution for modeling wind speed data. They also proposed methods to estimate wind power density. Wang [29] developed interval inference methods for estimating general lower-truncated distribution parameters based on double type II censoring samples. He et al. [30] presented an optimization design method for implementing reliability analysis with truncated normal random variables. Dörre [31] proposed a Bayesian estimation method for lifetime inference based on doubly truncated time-restricted samples. New Bayesian estimation methods about lifetime inferences can be found in the literature. Imani and Braga-Neto [32] proposed an approximate minimum mean square error filtering algorithm and smoothing algorithm based on the auxiliary particle filter method from sequential Monte Carlo theory. For overcoming the heavy-censoring problem for Weibull mixture parameters estimation, Ducros and Pamphile [33] proposed a Bayesian bootstrap method to implement reliability or warranty analysis based on nonhomogeneous lifetime samples. Jaheen and Okasha [34] proposed expected Bayesian estimation for the BurrXII model based on type II censoring. Han [35] studied the expected Bayesian estimation and its E-posterior risk for the failure rate of exponential distribution. Afify et al. [36] proposed maximum likelihood estimation and Bayesian estimation methods to estimate the parameters of the generalized odd log-logistic exponential distribution.

Motivation and Organization
Based on our best knowledge, no literature has discussed using ALT methods to infer the reliability for the doubly truncated three-parameter BurrXII distribution. Because the likelihood function based on ALT samples from the doubly truncated three-parameter BurrXII distribution is very complicated, the explicit forms of maximum likelihood estimators for the model parameters cannot be derived. Moreover, numerical methods via gradient algorithms cannot help to obtain the maximum likelihood estimates (MLEs) of parameters due to the divergence problem. In this study, we proposed two Bayesian estimation procedures to obtain Bayes estimators of the model parameters and the credible intervals of quantiles.
The rest of this paper is organized as follows: The doubly truncated three-parameter Burr type XII distribution and ALT model are presented in Section 2. Moreover, the maximum likelihood estimation method and two Bayesian estimation procedures using the Metropolis-Hastings algorithm via Gibbs sampling to implement the MCMC approach are developed. Additionally, the credible intervals of quantiles of the doubly truncated three-parameter BurrXII distribution at the normal use condition are obtained. Monte Carlo simulations are conducted for evaluating the performance of the two proposed Bayesian estimation procedures in Section 3. One numerical example about oil-well pumps is used in Section 4 for illustrating the applications of the proposed Bayesian estimation procedure. Some concluding remarks are given in Section 5.

The Statistical Model
Let the probability density function (PDF) and cumulative distribution function (CDF) of a reliable unit lifetime, X, be f (x|θ) and F(x|θ), respectively, where θ is the vector of distribution parameters. Given two positive constants ν and µ(≤ ν) for truncation, it is easy to show that Then, g(x|θ) and G(x|θ) can be defined as the doubly truncated version of f (x|θ) and F(x|θ) over the domain µ < x < ν, respectively. When the three-parameter BurrXII distribution with θ = (α, k, c) is considered, the f (x|θ) and F(x|θ) can be represented, respectively, by and where c and k are the inner and outer shape parameters, respectively, and α is the scale parameter. The doubly truncated version of the three-parameter BurrXII distribution, BurrXII µ,ν (θ), has the following PDF and CDF: and The survival function of the BurrXII µ,ν (θ) can be obtained by Okasha and Matter [27] had established some statistical properties for the doubly truncated three-parameter BurrXII distribution. Let x p be the pth quantile of BurrXII µ,ν (θ) such that G(x p |θ) = p for 0 < p < 1. It can be show that BurrXII µ,ν (θ) is a generalized version of the BurrXII(θ) distribution. When µ = 0 and ν → ∞, G(·|θ) reduces to F(·|θ) and G(·|θ) reduces to the upper truncated BurrXII(θ) as µ → 0. When ν → ∞, G(x|θ) reduces to the lower truncated BurrXII(θ). When α = 1, BurrXII µ,ν (θ) reduces to doubly truncated two-parameter BurrXII distribution.

The ALT Model and Parameter Estimation Methods
Let the lifetimes of reliable units be obtained through using a constant-stress ALT method with one accelerated variable, denoted by S. The stress levels of S are denoted by s 1 ≤ s 2 ≤ ... ≤ s m . Assume that the outer shape and scale parameters are dependent upon the stress via the link functions k i ≡ k(s i ) = b 0 + b 1 s i and α i ≡ α(s i ) = a 0 + a 1 s i , respectively, for i = 1, 2, ..., m, and the inner parameter is free from the stress. If the sample size for each stress of the ALT is n (that is, n 1 = n 2 = ... = n m = n) and the observed lifetimes under the ALT are denoted by x = {x ij , j = 1, 2, ..., n i , i = 1, 2, ..., m}, because the constant-stress ALT design method is used, the sample {x ij , j = 1, 2, ..., n i , i = 1, 2, ..., m} is independent. Let Θ = (a 0 , a 1 , c, b 0 , b 1 ). The likelihood and log-likelihood functions can be presented by and respectively. The maximum likelihood estimates (MLEs) of a 0 , a 1 , b 0 , b 1 , and c can be obtained through solving the likelihood equations a 0 = 0, a 1 = 0, b 0 = 0, b 1 = 0, and c = 0 simultaneously, where and No explicit forms of the MLEs,â 0 ,â 1 ,b 0 ,b 1 , andĉ can be found. Moreover, it is not tractable to obtain the values ofâ 0 ,â 1 ,b 0 ,b 1 , andĉ even by using numerical methods to simultaneously solve all five likelihood equations. One typical numerical method is the Newton-Rapson method (Newton iterative method). Many software packages are available for users to implement the Newton-Rapson method. Hence, the Newton-Rapson method has been widely used by statisticians and engineers to obtain MLEs via solving nonlinear likelihood equations. However, the Newton-Rapson method is a gradient method via the iterative procedure in which the solution quality depends on the selection of initial values of parameters. It is not easy to set up good initial values for all parameters in the current case. Hence, Bayes estimates are obtained to replace the MLEs for estimating model parameters.
Let the prior distribution of Θ be π(Θ), which is defined by Therefore, the posterior distribution of Θ, given data, can be represented by The full conditional posterior distribution for each model parameter can be expressed as follows: Panahi and Sayyareh [10] and Jaheen and Okasha [34] suggested to use gamma prior distributions as informative prior distributions to obtain Bayes estimates. In this study, the follow prior distributions are utilized for the proposed Bayesian estimation procedures: and When the aforementioned gamma distributions are replaced by π a 0 (a 0 ) ∝ constant, π c (c) ∝ constant, and π b 0 (b 0 ) ∝ constant, the prior distribution π(Θ) is non-informative. The obtained Bayes estimates via using non-informative prior distributions are closed to the MLEs of the model parameters. If users have sufficient knowledge to set up the hyperparameters, β i and δ i , i = 1, 2, 3, the obtained Bayes estimates via using informative prior distributions have better quality than those obtained via using non-informative prior distributions.
The Bayes estimator of Θ, denoted byΘ B , is the estimator that minimizes the Bayes risk function, is the loss incurred due to utilizingΘ B . In this paper, the square loss function, δ(Θ,Θ B ) = (Θ B − Θ) 2 is used for the Bayesian estimation method and the Bayes estimate of Θ, denoted byΘ B = (ã 0 ,ã 1 ,b 0 ,b 1 ,c), is the mean of the posterior distribution. The posterior distribution in Equation (15) is very complicated, and the explicit forms of the Bayes estimates are not available. Therefore, a Metropolis-Hastings algorithm via Gibbs sampling is proposed based on the following steps and Algorithm 1 to implement the MCMC method to obtain the Bayes estimate of Θ.

Initial
Step: Let i = 0 and a 1 , and c (0) be the initial states of a 0 , a 1 , b 0 , b 1 , and c, respectively.
Step 2.1: Generate a Step 2.2: Generate a ( * ) Step 2.4: Generate b ( * ) Step 2.5: Generate c ( * ) ∼ q c (c ( * ) |c (i) ) and u ∼ U(0, 1). Update c (i+1) by Step 3: The Bayes estimates can be obtained byθ a 1 , b 0 , b 1 , and c, where the first M(< N) chains are used for burn-in and all burn-in chains will be removed from the computation to obtain the Bayes estimates.
Two Bayesian estimation procedures, denoted by procedure BE-I and procedure BE-II, are given as the following to obtain reliable Bayes estimates. i + 1 ← i. In many occasions, users may not have sufficient knowledge about the model parameters and they need to use wide domains for D a 0 , D a 1 , D b 0 , D b 1 , and D c . Using uniform distributions with wide domains as transition probabilities for implementing the Metropolis-Hastings algorithm via Gibbs sampling will make the Markov chains converge slowly with less efficiency to search accurate and precise Bayes estimates. Hence, the following procedure BE-II is suggested to improve the performance of procedure BE-I.

Procedure BE-II:
Step  (6) can be used to establish the empirical distribution of the quantile estimator of x p , denoted byx p , at the normal use condition. Denote the empirical distribution ofx p byFx p (t). The γth and (1 − γ)th quantiles, t γ =F −1 x p (γ) and t 1−γ =F −1 x p (1 − γ), can be used to construct the (1 − 2γ)% creditable interval of x p . Hence, the (1 − 2γ)% creditable interval of x p can be denoted by (t γ , t 1−γ ).

Monte Carlo Simulations
In this section, the performance of procedures BE-I and BE-II is evaluated by conducting Monte Carlo simulations through using R codes. Procedures BE-I and BE-II with N = 12, 000 and M = 2000 are used to obtain the Bayes estimates; that is, 12,000 Markov chains are generated and the first 2000 Markov chains are removed for burn-in to obtain the Bayes estimates for procedure BE-I and procedure BE-II. All the simulation procedures are done by using R codes.
The constant-stress ALT in Monte Carlo simulations is set up to have two normalized stress levels: the low stress level is denoted by s 1 = s L = 0.45 and the high stress level is denoted by s 2 = s H = 1. The parameters a 0 = 12, a 1 = −5, b 0 = 6, and b 1 = 2 are used in the link functions k i = b 0 + b 1 s i and α i = a 0 + a 1 s i to generate samples from the doubly truncated three-parameter BurrXII distribution with parameters k i , α i , and c = 2.5 for i = 1, 2, µ = 0, and ν = 20. It can be shown that the normal use condition sample follows the doubly truncated BurrXII(c, k 0 , α 0 ) with k 0 = b 0 = 6, α 0 = a 0 = 12, c = 2.5, µ = 0, and ν = 20. In the simulation study, ALT samples with sizes (n 1 , n 2 ) = (10, 10), (20,20), (30,30), and (50, 50) are generated from the doubly truncated BurrXII(c, k i , α i ), i = 1, 2. Procedures BE-I and BE-II are used to search the Bayes estimates of the model parameters. Let , and D c = {1 ≤ a 0 ≤ 10} be the domains of the model parameters to implement the Metropolis-Hastings algorithm via Gibbs sampling described in Section 2. We hope that the Bayes estimate is close to the MLEs. Hence, non-informative prior distributions are used. The Bayes estimate can be obtained based on the sample mean from the last 10,000 Markov chains after removing the leading 2000 Markov chains for burn-in.
Repeat procedures BE-I and BE-II 1000 times, respectively. The relative bias (RB) and relative square root of mean square error (RsqMSE) of parameter δ are obtained based on the following equations: and 2 , andδ can be a 0 , a 1 , b 0 , b 1 , or c. The simulation results are reported in Tables 1 and 2. In view of Tables 1 and 2, we find that procedure BE-II outperforms procedure BE-I with smaller RB and a smaller RsqMSE for almost all the cases shown in Tables 1 and 2. These results indicate that the obtained Bayes estimates via using procedure BE-II are closer to their true values than the obtained Bayes estimates via using procedure BE-I, generally. Moreover, procedure BE-II can provide more reliable estimation results with a smaller RsqMSE than procedure BE-I. The boxplots in Figures 1-5 also indicate that the obtained Bayes estimates via using procedure BE-II are more reliable than those obtained via using procedure BE-I.     Figure 4. The boxplots of 1000 Bayes estimates of b 1 , where "n.I" and "n.II" indicate that procedures BE-I and BE-II are implemented with sample size n, respectively.
n.I=10 n.II=10 n.I=20 n.II=20 n.I=30 n.II=30 n.I=50 n.II=50 If users have sufficient knowledge to set up the hyperparameters and to use informative prior distributions to implement the proposed Bayesian estimation procedures, it is easier to obtain reliable Bayes estimates of a 0 , a 1 , b 0 , b 1 , and c. An additional simulation study is conducted to verify the performance of the Bayesian estimation procedures through using informative prior distributions. Following the parameters a 0 = 12, a 1 = −5, b 0 = 6, b 1 = 2, and c = 2.5 that are used for simulation in Tables 1 and 2, constant-stress ALT samples of n 1 = n 2 = 50 were generated from BurrXII µ=0,ν=20 (c, k i , α i ), where k i = b 0 + b 1 s i and α i = a 0 + a 1 s i , i = 1, 2, to implement the Metropolis-Hastings algorithm via Gibbs sampling in Algorithm 1 with s 1 = s L = 0.45 and s 2 = s H = 1 and the transition probabilities of N(a 0 , 1), N(a 1 , 1), N(b 0 , 1), N(b 1 , 1), and N(c, 1) were used to generate the values of a 0 , a 1 , b 0 , b 1 , and c. Denote the first scenario of simulation with β 1 = δ 1 = a 0 , β 2 = δ 2 = c, and β 3 = δ 3 = b 0 by Infor-1, and denote the second scenario of simulation with β 1 = a 0 , β 2 = c, β 3 = b 0 , and δ 1 = δ 2 = δ 3 = 1 by Infor-2. The Metropolis-Hastings algorithm via Gibbs sampling in Algorithm 1 with N = 12, 000 and M = 2000 was repeated 1000 times to obtain 1000 Bayes estimates of a 0 , a 1 , b 0 , b 1 , and c. Please note that we use normal distributions as transition probabilities to generate the parameter values. Hence, the Metropolis-Hastings algorithm via Gibbs sampling and using informative prior distributions are different from procedure BE-I. We use procedure BE-III to denote the Bayesian estimation procedure that uses the Metropolis-Hastings algorithm via Gibbs sampling and by using informative prior distributions.
The RB and RsqMSEs of each Bayes estimator are evaluated based on the obtained 1000 Bayes estimates of a 0 , a 1 , b 0 , b 1 , and c through using procedure BE-III. All simulation results are reported in Table 3. From Table 3, we can find that procedure BE-III can be used to quickly obtain reliable Bayes estimates with one step based on ALT samples. Unlike procedure BE-II being time consuming for implementation due to using two steps to obtain Bayes estimates, procedure BE-III is efficient in performing Bayesian estimation for saving computation time. Because procedures BE-II and BE-III are very competitive, procedure BE-II using non-informative prior distributions is recommended to obtain reliable Bayes estimates of the model parameters when users do not have good knowledge to set up hyperparameters in practical applications. For helping users to have a clear picture to use the proposed Bayesian estimation methods, a flowchart is given in Figure 6 as a guideline. Table 3. The RB and RsqMSEs of Bayes estimates via using procedure BE-III for n 1 = n 2 = 50.

An Example
Xin et al. (2018) proposed a Bayesian estimation method to infer the reliability of oil-well pumps using the BurrXII(θ). For sucker-rod oil pumping systems, the most important lifting equipment is the oil-well pump, but it is weak and wears out over time due to cyclic loading, liquid corrosion, or sand wear during the operation. Fatigue fracture or wear leakage could cause a critical failure for the oil-well pump.  used BurrXII(θ) to model the lifetimes of oil-well pumps with a type II censoring scheme. They showed that BurrXII(c = 1.982, k = 5.313, α = 5.694) can be a good distribution to model the lifetimes of oil-well pumps in years. Assume that the constant-stress ALT method with two stresses, s 1 = 0.45 and s 2 = 1, and parameters a 0 = 5.694, a 1 = −2.5, b 0 = 5.313, and b 1 = 2 is used to save the test time and cost on testing the oil-well pumps that do not fail within 2 months or 1/6 year. Two ALT samples, each with size 50, were generated from BurrXII µ,ν (c = 1.982, k 1 , α 1 ) and BurrXII µ,ν (c = 1.982, k 2 , α 2 ) with µ = 1/6 and ν = ∞, respectively, and displayed in Table 4 for s 1 = 0.45 and s 2 = 1. The plots of hazard rates for these two stress levels are given in Figure 7. In view of Figure 7, it can be noticed that the high-stress level results in higher hazard rate than the low-stress level does. Let the domains of a 0 , a 1 , b 0 , b 1 , and c be Compared with the Markov chains that are obtained using domain I and III, lower update rates for the Markov chains ofb 0 andb 1 , which are obtained using domain II, are found in Figures 10 and 11. This fact indicates that the Markov chains based on using domain II could not generate good Bayes estimates. Hence, we prefer to use the Markov chains obtained based on using domain I to search the Bayes estimates of the model parameters. The Bayes estimatesα =ã 0 = 6.408,k =b 0 = 6.397, andc = 2.084 can be used to infer the life quality of oil-well pumps at the normal use condition.
Assume that we would like to infer the median lifetime of oil-well pumps; the empirical distribution can be established based on the Markov chain forx p , which can be obtained using Equation (6) and the 10,000 Markov chains ofã 0 ,ã 1 ,b 0 ,b 1 , andc. The histogram of 10,000 Markov chains forx 0.5 is given in Figure 13, and the Bayes estimate of x 0.5 isx 0.5 = 2.8. The 95% credible interval is (1.028, 4.546), which covers the true x 0.5 = 2.115.

Concluding Remarks
In this study, the doubly truncated three-parameter BurrXII distribution is used to model the lifetimes of reliable units and the lifetimes of units are collected under the constant-stress ALT method to save the test time and sample resources. Because the maximum likelihood estimators of the model parameters are difficult to be obtained using gradient algorithms to simultaneously solve all likelihood equations, two Bayesian estimation procedures are proposed to obtain the Bayes estimates of the model parameters through using the Metropolis-Hastings algorithm via Gibbs sampling for generating Markov chains. The obtained Markov chains are used to establish the empirical distribution of the lifetime quantile at the normal use condition. Moreover, the Bayes estimate and creditable interval of the pth lifetime quantile at the normal use condition are obtained.
Intensive simulation studies were conducted to verify the performance of two proposed Bayesian estimation procedures. We found that procedure BE-II outperforms procedure BE-I. Hence, procedure BE-II with non-informative prior distributions is recommended to obtain reliable Bayes estimates of the model parameters when users do not have sufficient knowledge to set up hyperparameters. For an accelerated life testing that has two stress levels, low and high stress levels, at least 50 units for each stress level of the ALT are required to use procedure BE-II to obtain reliable Bayes estimates of the model parameters and to establish the empirical distribution of the estimator of a quantile at the normal use condition. The proposed procedure BE-II can release the model identification problem caused by using two shape parameters in the doubly truncated three-parameter BurrXII distribution.
A numerical example about the lifetimes of oil-well pumps is used for illustrating the applications of the proposed Bayesian estimation methods. The doubly truncated three-parameter BurrXII distribution is a generalized version of the BurrXII distribution. How to set up hyperparameters in the prior distribution to obtain reliable Bayes estimators is an important issue. Machine learning technique-based numerical methods could be competitive with the proposed methods. Extending the proposed Bayesian estimation procedures to a generalized version of other lifetime distributions under the constant-stress or step-stress ALT methods can be another good topic for study. These two topics are interesting and will be studied in the future. Funding: This study is supported by the grant of Ministry of Science and Technology, Taiwan MOST 108-2221-E-032-018-MY2.

Conflicts of Interest:
The authors declare no conflict of interest.