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

Abstract

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.

1. Introduction

1.1. 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 to infer the parameters of the two-parameter BurrXII distribution with general progressive type II censoring samples. Ismail and Khalid [7] studied the parameter estimation method using expectation-maximization algorithm for BurrXII distributions. More applications of the BurrXII distributions can be found in Thupeng [8] for modeling the maximum levels of nitrogen dioxide; in Nadar and Papadopoulos [9] for modeling the record data through using Bayesian estimation methods; in Panahi and Sayyareh [10] for predicating type II censored order statistics; in Xin et al. [11], in which the authors developed a Bayesian estimation procedure through using the Markov Chain Monte Carlo (MCMC) approach for parameter estimation based on type II censored samples from the three-parameter BurrXII distribution; and in Xin et al. [12], in which they used particle swarm optimization methods to infer the three-parameter BurrXII distribution parameters based on type I censoring samples. Chaturvedi et al. [13] proposed preliminary test estimators for the power parameters of the three-parameter BurrXII distribution and the reliability functions based on the maximum likelihood estimator, uniformly minimum variance unbiased estimator, and best invariant estimator as well as proposed empirical Bayes estimator. Tahir et al. [14] studied Bayesian and maximum likelihood estimation methods for the three-component mixture of BurrXII distributions based on doubly censored data. EL-Sagheer et al. [15] studied Bayesian and maximum likelihood estimation methods for the three-parameter BurrXII distribution with a unified hybrid censoring scheme. Byrnes et al. [16] studied the performance of using different loss functions to implement the MCMC Bayesian estimation method for estimating the parameter $\delta =P(X<Y)$, where X and Y follow two-parameter BurrXII distributions. EL-Sagheer et al. [17] proposed Bayesian inference methods for the three-parameter BurrXII distribution based on randomly censoring samples.

Because of the advance technologies on manufacturing, most modern products are highly reliable. It is difficult to collect the complete lifetime of a modern product at the normal use condition for reliability evaluations. Under this situation, an accelerated life test (ALT) could be an efficient procedure to save test time and cost because engineers can accelerate the speed of a unit’s failure by enhancing the stress used during an ALT. Constant-stress ALT, step-stress ALT, and ramp-stress ALT are three widely used ALT processes in practical applications. Compared with a step-stress ALT, the constant-stress ALT method does not need a memory assumption but requires more sample resources for life testing. Moreover, the constant-stress ALT method is easy for implementation. Ali Mouse [18] proposed an empirical Bayes estimation method for BurrXII distributions under an ALT with type II censoring. Ahmad and Islam [19] proposed optimal design methods of ALT for BurrXII distributions with type I censoring and periodic inspection. Abd-Elfattah et al. [20] studied maximum likelihood estimation methods for BurrXII distributions based on step-stress partially ALT samples and developed the asymptotic variance-covariance matrix of the estimators. Under the progressive type II censoring, Abdel-Hamid  [21] studied the maximum likelihood estimation method and provided a Fisher information matrix to implement reliability inferences for BurrXII distributions under the constant-partially accelerated life testing. Srivastava and Mittal [22] proposed optimal designs for BurrXII distributions through using a multi-objective ramp-stress ALT method. Zhao et al. [23] proposed a simple constant-stress ALT to collect type I progressively hybrid censored samples from BurrXII distribution and obtained the maximum likelihood estimators of the model parameters via using numerical methods. They also obtained approximate confidence intervals of the distribution parameters through using the normal approximation and bootstrap methods. Ismail and Al-Habardi [24] proposed interval inference methods to estimate BurrXII distribution parameters under the step-stress partially accelerated life testing. Prakash [25] proposed methods to find approximate confidence intervals through using normal approximations and bootstrap procedure under the constant-partially ALT. Ismail and Al-Habardi [26] proposed optimal plans for step-stress ALT for the BurrXII distribution with failure-censored samples.

In some occasions, truncation is a required condition for lifetime modeling in practical environments due to the ability to record data or it is even known that the failure occurs over or below a threshold or within a specific timespan. Okasha and Matter [27] used three-parameter BurrXII 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.

1.2. 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.

2. The Doubly Truncated Three-Parameter Burrxii Distribution and ALT Model

2.1. The Statistical Model

Let the probability density function (PDF) and cumulative distribution function (CDF) of a reliable unit lifetime, X, be $f\left(x\right|\theta )$ and $F\left(x\right|\theta )$, respectively, where $\theta$ is the vector of distribution parameters. Given two positive constants $\nu$ and $\mu (\le \nu )$ for truncation, it is easy to show that ${\int }_{\mu }^{\nu }f\left(x\right|\theta )dx=F\left(\nu \right|\theta )-F\left(\mu \right|\theta )$ and $ßn{t}_{\mu }^{\nu }g\left(x\right|\theta )dx=1$, where $g\left(x\right|\theta )=\frac{f\left(x\right|\theta )}{F\left(\nu \right|\theta )-F(\mu \left|\theta \right)}$ for $0<\mu \le x\le \nu$. Let $G\left(x\right|\theta )={\int }_{\mu }^{x}g\left(t\right|\theta )dt=\frac{F\left(x\right|\theta )-F(\mu \left|\theta \right)}{F\left(\nu \right|\theta )-F(\mu \left|\theta \right)}$ for $\mu <x<\nu$. Then, $g\left(x\right|\theta )$ and $G\left(x\right|\theta )$ can be defined as the doubly truncated version of $f\left(x\right|\theta )$ and $F\left(x\right|\theta )$ over the domain $\mu <x<\nu$, respectively. When the three-parameter BurrXII distribution with $\theta =(\alpha ,k,c)$ is considered, the $f\left(x\right|\theta )$ and $F\left(x\right|\theta )$ can be represented, respectively, by

$$f\left(x\right|\theta )=\frac{ck}{\alpha }{\left(\frac{x}{\alpha }\right)}^{c-1}{\left(1+{\left(\frac{x}{\alpha }\right)}^c\right)}^{-(k+1)},c,k,\alpha ,x>0,$$

(1)

and

$$F\left(x\right|\theta )=1-{\left(1+{\left(\frac{x}{\alpha }\right)}^c\right)}^{-k},c,k,\alpha ,x>0,$$

(2)

where c and k are the inner and outer shape parameters, respectively, and $\alpha$ is the scale parameter. The doubly truncated version of the three-parameter BurrXII distribution, BurrXII${}_{\mu ,\nu }\left(\theta \right)$, has the following PDF and CDF:

$$g\left(x\right|\theta )=\frac{1}{d_{\mu ,\nu }}\frac{ck}{\alpha }{\left(\frac{x}{\alpha }\right)}^{c-1}{\left(1+{\left(\frac{x}{\alpha }\right)}^c\right)}^{-(k+1)},c,k,\alpha >0,\mu \le x\le \nu ,$$

(3)

and

$$\begin{array}{ccc}\hfill G\left(x\right|\theta )& =& \frac{1}{d_{\mu ,\nu }}\left\{{\left(1+{\left(\frac{\mu }{\alpha }\right)}^c\right)}^{-k}-{\left(1+{\left(\frac{x}{\alpha }\right)}^c\right)}^{-k}\right\},c,k,\alpha >0,\mu \le x\le \nu ,\hfill \end{array}$$

(4)

respectively, where $d_{\mu ,\nu }={\left(1+{\left(\frac{\mu }{\alpha }\right)}^c\right)}^{-k}-{\left(1+{\left(\frac{\nu }{\alpha }\right)}^c\right)}^{-k}$. The survival function of the BurrXII${}_{\mu ,\nu }\left(\theta \right)$ can be obtained by

$$S_g\left(x\right|\theta )=1-G\left(x\right|\theta )=1-\frac{1}{d_{\mu ,\nu }}\left\{{\left(1+{\left(\frac{\mu }{\alpha }\right)}^c\right)}^{-k}-{\left(1+{\left(\frac{x}{\alpha }\right)}^c\right)}^{-k}\right\},c,k,\alpha >0,\mu \le x\le \nu .$$

(5)

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${}_{\mu ,\nu }\left(\theta \right)$ such that $G\left(x_p\right|\theta )=p$ for $0<p<1$. It can be show that

$$x_p=\alpha {\left({\left({\left(1+{\left(\frac{\mu }{\alpha }\right)}^c\right)}^{-k}-p\times d_{\mu ,\nu }\right)}^{-1/k}-1.0\right)}^{1/c}.$$

(6)

BurrXII${}_{\mu ,\nu }\left(\theta \right)$ is a generalized version of the BurrXII$\left(\theta \right)$ distribution. When $\mu =0$ and $\nu \to \infty$, $G(·|\theta )$ reduces to $F(·|\theta )$ and $G(·|\theta )$ reduces to the upper truncated BurrXII$\left(\theta \right)$ as $\mu \to 0$. When $\nu \to \infty$, $G\left(x\right|\theta )$ reduces to the lower truncated BurrXII$\left(\theta \right)$. When $\alpha =1$, BurrXII${}_{\mu ,\nu }\left(\theta \right)$ reduces to doubly truncated two-parameter BurrXII distribution.

2.2. 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\le s_2\le \dots \le s_m$. Assume that the outer shape and scale parameters are dependent upon the stress via the link functions $k_i\equiv k\left(s_i\right)=b_0+b_1{s}_i$ and ${\alpha }_i\equiv \alpha \left(s_i\right)=a_0+a_1{s}_i$ , respectively, for $i=1,2,\dots ,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=\dots =n_m=n$) and the observed lifetimes under the ALT are denoted by $x=\{x_{ij},j=1,2,\dots ,n_i,i=1,2,\dots ,m\}$, because the constant-stress ALT design method is used, the sample $\{x_{ij},j=1,2,\dots ,n_i,i=1,2,\dots ,m$} is independent. Let $\mathsf{\Theta }=(a_0,a_1,c,b_0,b_1)$. The likelihood and log-likelihood functions can be presented by

$$\begin{array}{ccc}\hfill L\left(\mathsf{\Theta }\right|x)& =& \prod _{i=1}^m\prod _{j=1}^n\frac{\frac{c{k}_i}{{\alpha }_i}{\left(\frac{x_{ij}}{{\alpha }_i}\right)}^{c-1}{\left(1+{\left(\frac{x_{ij}}{{\alpha }_i}\right)}^c\right)}^{-k_i-1}}{{\left(1+{\left(\frac{\mu }{{\alpha }_i}\right)}^c\right)}^{-k_i}-{\left(1+{\left(\frac{\nu }{{\alpha }_i}\right)}^c\right)}^{-k_i}},\hfill \end{array}$$

(7)

and

$$\begin{array}{ccc}\hfill \ell \left(\mathsf{\Theta }\right|x)& =& m\times nln\left(c\right)+n{\displaystyle \sum _{i=1}^m}ln\left(k_i\right)+(c-1){\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}ln\left(x_{ij}\right)\hfill \\ & & -cn{\displaystyle \sum _{i=1}^m}ln\left({\alpha }_i\right)-{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}(k_i+1)ln\left[1+{\left(\frac{x_{ij}}{{\alpha }_i}\right)}^c\right]\hfill \\ & & -n{\displaystyle \sum _{i=1}^m}ln\left\{{\left(1+{\left(\frac{\mu }{{\alpha }_i}\right)}^c\right)}^{-k_i}-{\left(1+{\left(\frac{\nu }{{\alpha }_i}\right)}^c\right)}^{-k_i}\right\},\hfill \end{array}$$

(8)

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 ${\ell }_{a_0}=0$, ${\ell }_{a_1}=0$, ${\ell }_{b_0}=0$, ${\ell }_{b_1}=0$, and ${\ell }_c=0$ simultaneously, where

$$\begin{array}{ccc}\hfill {\ell }_{a_0}=\frac{\partial \ell \left(\mathsf{\Theta }\right|x)}{\partial a_0}& =& -cn{\displaystyle \sum _{i=1}^m}\frac{1}{{\alpha }_i}+c{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}\frac{(k_i+1)x_{ij}^c}{{\alpha }_i({\alpha }_i^c+x_{ij}^c)}\hfill \\ & & -cn{\displaystyle \sum _{i=1}^m}{k}_i\frac{{[{\alpha }_i^c+{\mu }^c]}^{-k_i-1}{\mu }^c-{[{\alpha }_i^c+{\nu }^c]}^{-k_i-1}{\nu }^c}{{\alpha }_i[{({\alpha }_i^c+{\mu }^c)}^{-k_i}-{({\alpha }_i^c+{\nu }^c)}^{-k_i}]},\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill {\ell }_{a_1}=\frac{\partial \ell \left(\mathsf{\Theta }\right|x)}{\partial a_1}& =& -cn{\displaystyle \sum _{i=1}^m}\frac{s_i}{{\alpha }_i}+c(k_i+1){\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}\frac{s_i{x}_{ij}^c}{{\alpha }_i^{c+1}}\hfill \\ & & -cn{\displaystyle \sum _{i=1}^m}{s}_i{k}_i\frac{{[{\alpha }_i^c+{\mu }^c]}^{-k_i-1}{\mu }^c-{[{\alpha }_i^c+{\nu }^c]}^{-k_i-1}{\nu }^c}{{\alpha }_i[{({\alpha }_i^c+{\mu }^c)}^{-k_i}-{({\alpha }_i^c+{\nu }^c)}^{-k_i}]},\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill {\ell }_{b_0}=\frac{\partial \ell \left(\mathsf{\Theta }\right|x)}{\partial b_0}& =& n{\displaystyle \sum _{i=1}^m}\frac{1}{k_i}-{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}ln\left[1+{\left(\frac{x_{ij}}{{\alpha }_i}\right)}^c\right]\hfill \\ & & +n{\displaystyle \sum _{i=1}^m}\frac{{({\alpha }_i^c+{\mu }^c)}^{-k_i}ln\left[1+{\left(\frac{\mu }{{\alpha }_i}\right)}^c\right]-{({\alpha }_i^c+{\nu }^c)}^{-k_i}ln\left[1+{\left(\frac{\nu }{{\alpha }_i}\right)}^c\right]}{{({\alpha }_i^c+{\mu }^c)}^{-k_i}-{({\alpha }_i^c+{\nu }^c)}^{-k_i}},\hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill {\ell }_{b_1}=\frac{\partial \ell \left(\mathsf{\Theta }\right|x)}{\partial b_1}& =& n{\displaystyle \sum _{i=1}^m}\frac{s_i}{k_i}-{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}{s}_{i}ln\left[1+{\left(\frac{x_{ij}}{{\alpha }_i}\right)}^c\right]\hfill \\ & & +n{\displaystyle \sum _{i=1}^m}{s}_i\frac{{({\alpha }_i^c+{\mu }^c)}^{-k_i}ln\left[1+{\left(\frac{\mu }{{\alpha }_i}\right)}^c\right]-{({\alpha }_i^c+{\nu }^c)}^{-k_i}ln\left[1+{\left(\frac{\nu }{{\alpha }_i}\right)}^c\right]}{{({\alpha }_i^c+{\mu }^c)}^{-k_i}-{({\alpha }_i^c+{\nu }^c)}^{-k_i}}\hfill \end{array}$$

(12)

and

$$\begin{array}{ccc}\hfill {\ell }_c=\frac{\partial \ell \left(\mathsf{\Theta }\right|x)}{\partial c}& =& \frac{mn}{c}+{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}ln\left(x_{ij}\right)-n{\displaystyle \sum _{i=1}^n}ln\left({\alpha }_i\right)\hfill \\ & & -{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}\frac{k_i+1}{1+{\left(\frac{x_{ij}}{{\alpha }_i}\right)}^c}{\left(\frac{x_{ij}}{{\alpha }_i}\right)}^{c}ln\left(\frac{x_{ij}}{{\alpha }_i}\right)+n{\displaystyle \sum _{i=1}^m}{k}_i.\hfill \end{array}$$

(13)

No explicit forms of the MLEs, ${\widehat{a}}_0$, ${\widehat{a}}_1$, ${\widehat{b}}_0$, ${\widehat{b}}_1$, and $\widehat{c}$ can be found. Moreover, it is not tractable to obtain the values of ${\widehat{a}}_0$, ${\widehat{a}}_1$, ${\widehat{b}}_0$, ${\widehat{b}}_1$, and $\widehat{c}$ 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 $\mathsf{\Theta }$ be $\pi \left(\mathsf{\Theta }\right)$, which is defined by

$$\pi \left(\mathsf{\Theta }\right)={\pi }_{a_0}\left(a_0\right)\times {\pi }_{a_1}\left(a_1\right)\times {\pi }_c\left(c\right)\times {\pi }_{b_0}\left(b_0\right)\times {\pi }_{b_1}\left(b_1\right).$$

(14)

Therefore, the posterior distribution of $\mathsf{\Theta }$, given data, can be represented by

$$\pi \left(\mathsf{\Theta }\right|x)\propto \pi (\mathsf{\Theta })\times L(\mathsf{\Theta }\left|x\right).$$

(15)

The full conditional posterior distribution for each model parameter can be expressed as follows:

$$\begin{array}{ccc}\hfill {\pi }_{a_0}\left(a_0\right|a_1,b_0,b_1,c,x)& \propto & {\pi }_{a_0}\left(a_0\right)\times L\left(\mathsf{\Theta }\right|x),\hfill \\ \hfill {\pi }_{a_1}\left(a_1\right|a_0,b_0,b_1,c,x)& \propto & {\pi }_{a_1}\left(a_1\right)\times L\left(\mathsf{\Theta }\right|x),\hfill \\ \hfill {\pi }_c\left(c\right|a_0,a_1,b_0,b_1,x)& \propto & {\pi }_c\left(c\right)\times L\left(\mathsf{\Theta }\right|x),\hfill \\ \hfill {\pi }_{b_0}\left(b_0\right|a_0,a_1,b_1,c,x)& \propto & {\pi }_{b_0}\left(b_0\right)\times L\left(\mathsf{\Theta }\right|x),\hfill \end{array}$$

and

$${\pi }_{b_1}\left(b_1\right|a_0,a_1,b_0,c,x)\propto {\pi }_{b_1}\left(b_1\right)\times L\left(\mathsf{\Theta }\right|x).$$

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:

$$\begin{array}{ccc}\hfill {\pi }_{a_0}\left(a_0\right)& =& \frac{{\beta }_1^{{\delta }_1}}{\mathsf{\Gamma }\left({\delta }_1\right)}{a}_0^{{\delta }_1-1}exp\{-{\beta }_1{a}_0\},a_0,{\beta }_1\ge 0,{\delta }_1>0,\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill {\pi }_{a_1}\left(a_1\right)& \propto & \mathrm{constant}\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill {\pi }_c\left(c\right)& =& \frac{{\beta }_2^{{\delta }_2}}{\mathsf{\Gamma }\left({\delta }_2\right)}{c}^{{\delta }_2-1}exp\{-{\beta }_{2}c\},c,{\beta }_2\ge 0,{\delta }_2>0,\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill {\pi }_{b_0}\left(b_0\right)& =& \frac{{\beta }_3^{{\delta }_3}}{\mathsf{\Gamma }\left({\delta }_3\right)}{b}_0^{{\delta }_3-1}exp\{-{\beta }_3{b}_0\},b_0,{\beta }_3\ge 0,{\delta }_3>0,\hfill \end{array}$$

(19)

and

$${\pi }_{b_1}\left(b_1\right)\propto \mathrm{constant}.$$

(20)

When the aforementioned gamma distributions are replaced by ${\pi }_{a_0}\left(a_0\right)\propto \mathrm{constant}$, ${\pi }_c\left(c\right)\propto \mathrm{constant}$, and ${\pi }_{b_0}\left(b_0\right)\propto \mathrm{constant}$, the prior distribution $\pi \left(\mathsf{\Theta }\right)$ 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, ${\beta }_i$ and ${\delta }_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 $\mathsf{\Theta }$, denoted by ${\tilde{\mathsf{\Theta }}}_B$, is the estimator that minimizes the Bayes risk function, $R(\mathsf{\Theta },{\tilde{\mathsf{\Theta }}}_B)=E_{\pi }\{\delta (\mathsf{\Theta },{\tilde{\mathsf{\Theta }}}_B)\}$, where $\delta (\mathsf{\Theta },{\tilde{\mathsf{\Theta }}}_B)$ is the loss incurred due to utilizing ${\tilde{\mathsf{\Theta }}}_B$. In this paper, the square loss function, $\delta (\mathsf{\Theta },{\tilde{\mathsf{\Theta }}}_B)={({\tilde{\mathsf{\Theta }}}_B-\mathsf{\Theta })}^2$ is used for the Bayesian estimation method and the Bayes estimate of $\mathsf{\Theta }$, denoted by ${\tilde{\mathsf{\Theta }}}_B=({\tilde{a}}_0,{\tilde{a}}_1,{\tilde{b}}_0,{\tilde{b}}_1,\tilde{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 $\mathsf{\Theta }$.

Initial Step: 

Let $i=0$ and $a_0^{\left(0\right)}$, $a_1^{\left(0\right)}$, $b_0^{\left(0\right)}$, $b_1^{\left(0\right)}$, and $c^{\left(0\right)}$ be the initial states of $a_0$, $a_1$, $b_0$, $b_1$, and c, respectively.

Step 1: 

Propose the transition probabilities $q_j\left(j^{(*)}\right|j^{\left(i\right)})$ from $j^{\left(i\right)}$ to $j^{(*)}$ for $j=a_0,a_1,b_0,b_1,c$.

Step 2: 

Implement Step 2.1 to Step 2.5 N times for $i=0,1,2,\dots ,N$, where N is a huge number.

Step 2.1: 

Generate $a_0^{(*)}\sim q_{a_0}(a_0^{(*)}|a_0^{\left(i\right)})$ and $u\sim U(0,1)$, where $U(0,1)$ is the uniform distribution over the domain of (0,1). Update $a_0^{(i+1)}$ by

$$a_0^{(i+1)}=\left\{\begin{array}{cc}{a}_0^{(*)},\hfill & u\le min\left\{1,\frac{\pi \left(a_0^{(*)}\right|a_1^{\left(i\right)},b_0^{\left(i\right)},b_1^{\left(i\right)},c^{\left(i\right)};x)q_{a_0}\left(a_0^{\left(i\right)}\right|a_0^{(*)})}{\pi \left(a_0^{\left(i\right)}\right|a_1^{\left(i\right)},b_0^{\left(i\right)},b_1^{\left(i\right)},c^{\left(i\right)};x)q_{a_0}\left(a_0^{(*)}\right|a_0^{\left(i\right)})}\right\},\hfill \\ a_0^{\left(i\right)},\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(21)

Step 2.2: 

Generate $a_1^{(*)}\sim q_{a_1}(a_1^{(*)}|a_1^{\left(i\right)})$ and $u\sim U(0,1)$. Update $a_1^{(i+1)}$ by

$$a_1^{(i+1)}=\left\{\begin{array}{cc}{a}_1^{(*)},\hfill & u\le min\left\{1,\frac{\pi \left(a_1^{(*)}\right|a_0^{(i+1)},b_0^{\left(i\right)},b_1^{\left(i\right)},c^{\left(i\right)};x)q_{a_1}\left(a_1^{\left(i\right)}\right|a_1^{(*)})}{\pi \left(a_1^{\left(i\right)}\right|a_0^{(i+1)},b_0^{\left(i\right)},b_1^{\left(i\right)},c^{\left(i\right)};x)q_{a_1}\left(a_1^{(*)}\right|a_1^{\left(i\right)})}\right\},\hfill \\ a_1^{\left(i\right)},\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(22)

Step 2.3: 

Generate $b_0^{(*)}\sim q_{b_0}(b_0^{(*)}|b_0^{\left(i\right)})$ and $u\sim U(0,1)$. Update $b_0^{(i+1)}$ by

$$b_0^{(i+1)}=\left\{\begin{array}{cc}{b}_0^{(*)},\hfill & u\le min\left\{1,\frac{\pi \left(b_0^{(*)}\right|a_0^{(i+1)},a_1^{(i+1)},b_1^{\left(i\right)},c^{\left(i\right)};x)q_{b_0}\left(b_0^{\left(i\right)}\right|b_0^{(*)})}{\pi \left(b_0^{\left(i\right)}\right|a_0^{(i+1)},a_1^{(i+1)},b_1^{\left(i\right)},c^{\left(i\right)};x)q_{b_0}\left(b_0^{(*)}\right|b_0^{\left(i\right)})}\right\},\hfill \\ b_0^{\left(i\right)},\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(23)

Step 2.4: 

Generate $b_1^{(*)}\sim q_{b_1}(b_1^{(*)}|b_1^{\left(i\right)})$ and $u\sim U(0,1)$. Update $b_1^{(i+1)}$ by

$$b_1^{(i+1)}=\left\{\begin{array}{cc}{b}_1^{(*)},\hfill & u\le min\left\{1,\frac{\pi \left(b_1^{(*)}\right|a_0^{(i+1)},a_1^{(i+1)},b_0^{(i+1)},c^{\left(i\right)};x)q_{b_1}\left(b_1^{\left(i\right)}\right|b_1^{(*)})}{\pi \left(b_1^{\left(i\right)}\right|a_0^{(i+1)},a_1^{(i+1)},b_0^{(i+1)},c^{\left(i\right)};x)q_{b_1}\left(b_1^{(*)}\right|b_1^{\left(i\right)})}\right\},\hfill \\ b_1^{\left(i\right)},\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(24)

Step 2.5: 

Generate $c^{(*)}\sim q_c(c^{(*)}|c^{\left(i\right)})$ and $u\sim U(0,1)$. Update $c^{(i+1)}$ by

$$c^{(i+1)}=\left\{\begin{array}{cc}{c}^{(*)},\hfill & u\le min\left\{1,\frac{\pi \left(c^{(*)}\right|a_0^{(i+1)},a_1^{(i+1)},b_0^{(i+1)},b_1^{(i+1)};x)q_c\left(c^{\left(i\right)}\right|c^{(*)})}{\pi \left(c^{\left(i\right)}\right|a_0^{(i+1)},a_1^{(i+1)},b_0^{(i+1)},b_1^{(i+1)};x)q_c\left(c^{(*)}\right|c^{\left(i\right)})}\right\},\hfill \\ c^{\left(i\right)},\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(25)

Step 3: 

The Bayes estimates can be obtained by ${\tilde{\theta }}_j=\frac{1}{N-M}{\sum }_{j=M+1}^N{\theta }_j^{\left(i\right)}$, ${\theta }_j=a_0,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.

Algorithm 1 The Metropolis-Hastings algorithm via Gibbs sampling.

1: Initial Step. 2: Step 1. 3: Step 2 with i=0. 4: loop: 5: if$i<N$then do Step 2.1-Step 2.5 6:     $i+1\leftarrow i$. 7:     goto loop. 8:     close; 9:   Remove the first N chains for burn-in. 10: Obtain Bayes estimates via using Step 3.

Procedure BE-I: 

Let the domain of the model parameters be $D_{a_0}$, $D_{a_1}$, $D_{b_0}$, $D_{b_1}$, and $D_c$. Uniform distributions over the domains of $D_{a_0}$, $D_{a_1}$, $D_{b_0}$, $D_{b_1}$, and $D_c$, respectively, are used to be the transition probabilities to implement the Metropolis–Hastings algorithm via Gibbs sampling to obtain Bayes estimates.

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 1: 

Obtain 100 sets of Bayes estimates through using procedure BE-I, and denote them by $\{{\tilde{a}}_{0,j},j=1,2,\dots ,100\}$, $\{{\tilde{a}}_{1,j},j=1,2,\dots ,100\}$, $\{{\tilde{b}}_{0,j},j=1,2,\dots ,100\}$, $\{{\tilde{b}}_{1,j},j=1,2,\dots ,100\}$, and $\{{\tilde{c}}_j,j=1,2,\dots ,100\}$. Find the mean of each set of Bayes estimates with 5% of them trimmed from each end. Denote the trimmed mean by ${\overline{a}}_0$, ${\overline{a}}_1$, ${\overline{b}}_0$, ${\overline{b}}_1$, and $\overline{c}$, respectively.

Step 2: 

Implement the Metropolis–Hastings algorithm via Gibbs sampling with the normal distributions, $N({\overline{a}}_0,1)$, $N({\overline{a}}_1,1)$, $N({\overline{b}}_0,1)$, $N({\overline{b}}_1,1)$, and $N(\overline{c},1)$, as the transition probabilities to obtain Bayes estimates. That is, the MCMC method is implemented based on the knowledge that is obtained from Step 1.

The Markov chains of $\{{\tilde{a}}_{0,j},j=1,2,\dots ,N-M\}$, $\{{\tilde{a}}_{1,j},j=1,2,\dots ,N-M\}$, $\{{\tilde{b}}_{0,j},j=1,2,\dots ,N-M\}$, $\{{\tilde{b}}_{1,j},j=1,2,\dots ,N-M\}$, and $\{{\tilde{c}}_j,j=1,2,\dots ,N-M\}$ and Equation (6) can be used to establish the empirical distribution of the quantile estimator of $x_p$, denoted by ${\tilde{x}}_p$, at the normal use condition. Denote the empirical distribution of ${\tilde{x}}_p$ by ${\widehat{F}}_{{\tilde{x}}_p}\left(t\right)$. The $\gamma$th and $(1-\gamma )$th quantiles, $t_{\gamma }={\widehat{F}}_{{\tilde{x}}_p}^{-1}\left(\gamma \right)$ and $t_{1-\gamma }={\widehat{F}}_{{\tilde{x}}_p}^{-1}(1-\gamma )$, can be used to construct the $(1-2\gamma )$% creditable interval of $x_p$. Hence, the $(1-2\gamma )$% creditable interval of $x_p$ can be denoted by $(t_{\gamma },t_{1-\gamma })$.

3. 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 ${\alpha }_i=a_0+a_1{s}_i$ to generate samples from the doubly truncated three-parameter BurrXII distribution with parameters $k_i$, ${\alpha }_i$, and $c=2.5$ for $i=1,2$, $\mu =0$, and $\nu =20$. It can be shown that the normal use condition sample follows the doubly truncated BurrXII$(c,k_0,{\alpha }_0)$ with $k_0=b_0=6$, ${\alpha }_0=a_0=12$, $c=2.5$, $\mu =0$, and $\nu =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,{\alpha }_i)$, $i=1,2$. Procedures BE-I and BE-II are used to search the Bayes estimates of the model parameters.

Let $D_{a_0}=\{1\le a_0\le 20\}$, $D_{a_1}=\{-10\le a_1\le 0\}$, $D_{b_0}=\{0\le b_0\le 10\}$, $D_{b_1}=\{1\le b_1\le 5\}$, and $D_c=\{1\le a_0\le 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 $\delta$ are obtained based on the following equations:

$$\mathrm{RB}=\frac{\overline{\tilde{\delta }}-\delta }{\delta }$$

(26)

and

$$\mathrm{RsqMSE}=\frac{\sqrt{MSE}}{\left|\delta \right|},$$

(27)

where $\overline{\tilde{\delta }}={\sum }_{i=1}^{1000}{\tilde{\delta }}_i/1000$, $MSE=\frac{1}{1000}{\sum }_{i=1}^{1000}{({\tilde{\delta }}_i-\delta )}^2$, and $\tilde{\delta }$ can be $a_0$, $a_1$, $b_0$, $b_1$, or c. The simulation results are reported in

Table 1. The relative bias (RB) of Bayes estimates for $n_1=n_2=$ 10, 20, 30, and 50.

$({\mathit{n}}_1,{\mathit{n}}_2)$	Procedure	${\mathit{a}}_0$	${\mathit{a}}_1$	${\mathit{b}}_0$	${\mathit{b}}_1$	c
(10, 10)	BE-I	−0.1415	0.1666	−0.0540	0.5020	0.9080
	BE-II	−0.1399	0.0460	−0.1136	0.0829	0.8557
(20, 20)	BE-I	−0.0710	0.1347	−0.0710	0.4849	0.7307
	BE-II	−0.1060	0.0443	−0.0959	0.1443	0.5903
(30, 30)	BE-I	−0.0247	0.0924	−0.0665	0.4818	0.5595
	BE-II	−0.0309	0.0545	−0.0378	0.2453	0.3749
(50, 50)	BE-I	0.0614	−0.0983	−0.1629	0.2560	0.2364
	BE-II	0.0564	−0.0849	−0.1117	0.4096	−0.0068

and

Table 2. The relative square root of mean square errors (RsqMSEs) of Bayes estimates for $n_1=n_2=$ 10, 20, 30, and 50.

$({\mathit{n}}_1,{\mathit{n}}_2)$	Procedure	${\mathit{a}}_0$	${\mathit{a}}_1$	${\mathit{b}}_0$	${\mathit{b}}_1$	c
(10, 10)	BE-I	0.3578	0.3238	0.2427	0.5961	1.0592
	BE-II	0.1422	0.0995	0.1300	0.1372	0.8575
(20, 20)	BE-I	0.3594	0.2971	0.2377	0.5674	0.9505
	BE-II	0.1073	0.0756	0.1126	0.1691	0.5916
(30, 30)	BE-I	0.3907	0.2952	0.2339	0.5639	0.8834
	BE-II	0.0353	0.0707	0.0586	0.2541	0.3790
(50, 50)	BE-I	0.3921	0.3458	0.2777	0.4370	0.7075
	BE-II	0.0601	0.0963	0.1165	0.4134	0.1063

. In view of Table 1 and Table 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 Table 1 and Table 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 Figure 1, Figure 2, Figure 3, Figure 4 and Figure 5 also indicate that the obtained Bayes estimates via using procedure BE-II are more reliable than those obtained via using procedure BE-I.

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 Table 1 and Table 2, constant-stress ALT samples of $n_1=n_2=50$ were generated from BurrXII${}_{\mu =0,\nu =20}(c,k_i,{\alpha }_i)$, where $k_i=b_0+b_1{s}_i$ and ${\alpha }_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 ${\beta }_1={\delta }_1=a_0$, ${\beta }_2={\delta }_2=c$, and ${\beta }_3={\delta }_3=b_0$ by Infor-1, and denote the second scenario of simulation with ${\beta }_1=a_0$, ${\beta }_2=c$, ${\beta }_3=b_0$, and ${\delta }_1={\delta }_2={\delta }_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. The RB and RsqMSEs of Bayes estimates via using procedure BE-III for $n_1=n_2=50$.

	Methods	${\mathit{a}}_0$	${\mathit{a}}_1$	${\mathit{b}}_0$	${\mathit{b}}_1$	c
RB	Infor-1	−0.1076	−0.0824	−0.2633	−0.1361	0.0195
	Infor-2	−0.1347	−0.0992	−0.3210	−0.1667	0.0186
RsqMSE	Infor-1	0.1113	−0.0910	0.2689	0.1549	0.0593
	Infor-2	0.1394	−0.1094	0.3283	0.1884	0.0552

. 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.

4. An Example

Xin et al. (2018) proposed a Bayesian estimation method to infer the reliability of oil-well pumps using the BurrXII$\left(\theta \right)$. 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. Xin et al. (2018) used BurrXII$\left(\theta \right)$ to model the lifetimes of oil-well pumps with a type II censoring scheme. They showed that BurrXII$(c=1.982,k=5.313,\alpha =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${}_{\mu ,\nu }(c=1.982,k_1,{\alpha }_1)$ and BurrXII${}_{\mu ,\nu }(c=1.982,k_2,{\alpha }_2)$ with $\mu =1/6$ and $\nu =\infty$, respectively, and displayed in

Table 4. The regenerated accelerated life test (ALT) samples of the oil-well pump example.

Low Stress
3.009, 1.434, 3.471, 3.937, 1.605, 2.015, 1.832, 1.501, 1.324, 0.825,
2.055, 2.847, 1.033, 1.612, 2.002, 2.020, 1.603, 1.080, 1.373, 1.849,
0.456, 0.903, 0.990, 1.089, 1.520, 1.151, 3.046, 0.457, 1.966, 0.841,
2.255, 2.542, 2.181, 1.637, 1.252, 0.907, 1.296, 1.304, 2.701, 0.556,
1.552, 3.132, 0.656, 1.097, 0.544, 2.814, 1.759, 1.041, 2.544, 1.853
High Stress
0.498, 1.871, 1.554, 0.679, 1.656, 1.225, 2.027, 1.458, 0.968, 0.667,
0.263, 1.436, 0.664, 2.435, 1.438, 0.638, 1.069, 1.042, 1.293, 0.386,
1.057, 2.197, 0.657, 1.352, 1.115, 0.587, 1.405, 0.635, 1.715, 1.592,
1.886, 0.850, 0.547, 0.783, 0.405, 1.675, 2.150, 0.743, 1.299, 0.766,
0.515, 1.281, 1.738, 2.615, 0.205, 1.058, 0.415, 0.223, 0.594, 1.687

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 $D_{a_0}=\{1\le a_0\le 10\}$, $D_{a_1}=\{-5\le a_1\le 0\}$, $D_{b_0}=\{1\le b_0\le 10\}$, $D_{b_1}=\{1\le b_1\le 5\}$, and $D_c=\{1\le c\le 10\}$. Because BurrXII${}_{\mu ,\nu }\left(\theta \right)$ contains two shape parameters, model identification could be a problem for modeling. One method to overcome this problem is to cut the domain of c into two disjoint subdomains, for example, $D_c^I=\{1\le c\le 5\}$ and $D_c^{II}=\{5<c\le 10\}$, and to search the Bayes estimates of the model parameters for the domains of $D_{a_0}$, $D_{a_1}$, $D_{b_0}$, $D_{b_1}$, and $D_c^I$ (labeled by domain I) and of $D_{a_0}$, $D_{a_1}$, $D_{b_0}$, $D_{b_1}$, and $D_c^{II}$ (labeled by domain II), respectively. Moreover, we also search the Bayes estimates of the model parameters based on the domains of $D_{a_0}$, $D_{a_1}$, $D_{b_0}$, $D_{b_1}$, and $D_c$ (labeled by domain III). Then, we screen out the best set from the two sets of Bayes estimates that are obtained from domain I and domain II to implement procedure BE-II. Non-informative prior distributions are used to implement the proposed Metropolis–Hastings algorithm via Gibbs sampling with $N=12,000$ and $M=2000$.

The Bayes estimates based on using domain I are ${\tilde{a}}_0=6.408$, ${\tilde{a}}_1=-2.246$, ${\tilde{b}}_0=6.397$, ${\tilde{b}}_1=2.567$, and $\tilde{c}=2.084$; the Bayes estimates via using domain II are ${\tilde{a}}_0=1.966$, ${\tilde{a}}_1=-0.913$, ${\tilde{b}}_0=1.001$, ${\tilde{b}}_1=0.001$, and $\tilde{c}=7.450$; and the Bayes estimates based on using domain III are ${\tilde{a}}_0=4.125$, ${\tilde{a}}_1=-1.899$, ${\tilde{b}}_0=4.164$, ${\tilde{b}}_1=1.511$, and $\tilde{c}=3.435$. The Markov chains obtained based on using domain I, domain II, and domain III are given in Figure 8, Figure 9, Figure 10, Figure 11 and Figure 12. Compared with the Markov chains that are obtained using domain I and III, lower update rates for the Markov chains of ${\tilde{b}}_0$ and ${\tilde{b}}_1$, which are obtained using domain II, are found in Figure 10 and Figure 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 $\tilde{\alpha }={\tilde{a}}_0=6.408$, $\tilde{k}={\tilde{b}}_0=6.397$, and $\tilde{c}=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 for ${\tilde{x}}_p$, which can be obtained using Equation (6) and the 10,000 Markov chains of ${\tilde{a}}_0$, ${\tilde{a}}_1$, ${\tilde{b}}_0$, ${\tilde{b}}_1$, and $\tilde{c}$. The histogram of 10,000 Markov chains for ${\tilde{x}}_{0.5}$ is given in Figure 13, and the Bayes estimate of $x_{0.5}$ is ${\tilde{x}}_{0.5}=2.8$. The 95% credible interval is (1.028, 4.546), which covers the true $x_{0.5}=2.115$.

5. 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.