Inferences for a Two-parameter Lifetime Distribution with Bathtub Shaped Hazard Based on Censored Data

We consider statistical inference of the unknown parameters of a two-parameter bathtub-shaped distribution (Chen, 2000) [Stat. & Prob. Letters 49 (2000) 155-161]. The inference will be conducted for Type-II censored and progressively Type-II censored data using the maximum likelihood and Bayes techniques. There are no explicit expressions for the estimators of the parameters. In the case of the maximum likelihood estimator (MLE), we propose a simple fixed point algorithm to compute the MLE and construct different confidence intervals and confidence regions of the unknown parameters. Bayes analyses of the unknown parameters are also discussed under fairly general priors for the unknown parameters. We propose to use the Markov Chain Monte Carlo (MCMC) and simulation-based technique to compute the Bayes estimates and the two-sided Bayesian probability intervals of the parameters. Also, we use the rejection sampling algorithm to produce the exact Bayes estimates. The methods developed will be applied in the analyses of two real data sets and a simulated data set. A Monte Carlo simulation is used to compare the results from the MLE and Bayes techniques.

The two parameter Weibull distribution is a very popular distribution for modeling phenomenon with monotone hazard rates.Its negatively and positively skewed density makes it an initial choice for modeling monotone hazard rates.However, Weibull distribution does not provide a reasonable fit in modeling phenomenon with non-monotone hazard rate such as the bathtub-shaped hazard rate.A distribution with a bathtub-shaped hazard rate provides an appropriate conceptual model for some electronic and mechanical products as well as the lifetime of humans.
There is extensive literature on parametric probability distributions with bathtub-shaped hazard rate function (e.g., Smith and Bain (1975), Leemis (1986), Gaver and Acar (1979), Hjorth (1980), and Mudhilkar and Srivastava (1993)).Among several extensions of the Weibull distribution, with 4 parameters, are the exponentiated generalized exponential linear distribution (Sarhan et. al, 2013 ) and the exponentiated modified Weibull extension distribution (Sarhan and Apaloo, 2013).Chen (2000) revisited a two parameter distribution with the following survival function (sf) where α > 0 is a parameter which does not affect the shape of the failure rate function and β > 0 is the shape parameter.The corresponding probability density function (pdf) is , x > 0. (3) The corresponding hazard rate function (hrf) of this distribution is h(x; α, β) = αβx β−1 e x β , x > 0. (4) The hrf can be either (1) an exponentially increasing function when β ≥ 1; or (2) of bathtub-shape when β < 1.For simplicity, we will refer to this distribution as TPBT(α, β).It has been observed in Nadarajah and Kotz (2007) and Pham and Lai (2007), that TPBT(α, β) is a special case of the distribution considered by Gurvich et al. (1997) and by Haynatzki et al. (2000).However, TPBT(α, β) has an advantage over all bathtub hazard shaped distributions since it has only two parameters.
The TPBT(α, β) distribution has been considered for a type-II censored data (Chen, 2000), type-II right censored data (Wu et al., 2004), for a progressively type-II censored data (Wu, 2008).In Chen (2000), exact onfidence interval for β and exact joint confidence regions for (α, β) were presented while Wu et al. (2004) presents statistical inference about the shape parameter.Sarhan et al. (2012) derived both the maximum likelihood estimates (MLE) and Bayes estimates (BE) of the two unknown parameters using a complete sample.They applied the simulationbased (SB) method, the Monte Carlo integration (MCI) technique and Lindley approximation (LA) method to approximate the Bayes estimates of the parameters and reported that the SB method was the best among the three methods used.
The aim of this paper is to discuss Bayes inferences of the TPBT distribution's parameters using different techniques when the data are complete, cenceored and progressively censored.We use Markov-Chain Monte-Carlo (MCMC), rejection sampling (RS) algorithms and compare them with the SB algorithm.
The rest of the paper is organized as follows.In section 2 we present the censored data scheme and the maximum likelihood estimation procedure.Confidence interval estimations are discussed in section 3. Bayes procedure is discussed in section 4. Two approximation techniques (the RS and MCMC) are presented in section 5. Progressively Type-II censored data is used in section 6.Two real and a simulated data sets are analyzed in section 7.
A large simulation study is provided in section 8 to compare the maximum likelihood technique with Bayesian techniques.Concluding remarks are made in section 9.

Maximum Likelihood Estimation
In Type-II censored data, it is assumed that n independent and identical units are put on the life test at the same time.The life test is terminated after a predetermined number of failures m results.For unit i, i = 1, 2, • • • , n, a pair of two quantities (X i , δ i ) is observed, where X i represents the testing time of unit i and δ i is a binary variable that takes the value of 1 when unit i has failed at time X i , or 0 if it is tested up to time X i without failure (X i is a censored time).
In this section we use the maximum likelihood method to estimate the two unknown parameters α and β, using Type-II censored data.Suppose (X 1 , δ 1 ), (X 2 , δ 2 ), . . ., (X n , , δ n ) is a random sample from TPBT(α, β), then the likelihood function of the observed data is where m = ∑ n i=1 δ i .The log-likelihood function becomes Taking derivatives with respect to α and β of (6), we obtain The second partial derivatives of L with respect to α and β are The information function is a two-by-two symmetric matrix For a relative maximum of L, which occurs at the MLE of the parameters α and β, the matrix I( α, β) must be positive definite.
Setting the right hand side of (7) to zero and solving for α, we get Substituting (10) into the right side of ( 8) and setting it equal to zero, we get the following non-linear equation in β The MLE of β, β, is the solution of (11) in β.A closed-form solution of (11) does not exist, so a numerical technique, e.g Newton-Raphson method, should be used to find β for any given data set.Once we get β, we can use (10) to get α = α( β).

Confidence Intervals (CIs)
We can use different techniques to approximate (1 − ϑ)100% confidence intervals of the two parameters α and β.
In the following, we describe two of such techniques.

Likelihood Intervals
A 100p% likelihood interval (LI), p ∈ (0, 1), for a parameter θ is the range of all values of θ for which the relative likelihood function, R(θ) = L(θ) L( θ) , is greater than or equal to p, for more details we refer to Kalbfleisch (1985).For simplicity, we can take the natural logarithm of R(θ), therefore the 100p% LI of θ becomes the range of all values of θ which satisfy that r(θ) ≥ ln p.Here, r(θ) is the log-relative likelihood function, given by r(θ) = L(θ) − L( θ).To determine the 100p% LI for θ, we can either use the the graph of r(θ) against θ or use the Newton-Raphson Method to get the bounds of the LI by solving r(θ) − ln p = 0.When θ is a vector of two unknown parameters, say θ = (α, β) ′ , as in the case studied here, the 100p% likelihood region is the set of parameter values (α, β) such that r(α, β) ≥ ln p, where r(α, β) is the joint log-relative likelihood function of α and β, r(α, β), given by r(α, Substituting from ( 6) into (12), we get where The 100p% likelihood contour is boundary of 100p% likelihood region, which is formed by the curve r(α, β) = ln p.The 100p% likelihood region for (α, β) is an approximate 100(1 − ϑ)% confidence region when p = e −χ 2 1 (ϑ)/2 , where χ 2 1 (ϑ) is the upper ϑ quantile of the chi square distribution with one degree of freedom.
Another way to get the 100p% LI is to use the maximum log-relative likelihood function of β, which is given by

Large-sample Intervals
The MLE of the parameters α and β are asymptotically normally distributed with means equal to the true values of α and β and variances given by the inverse of the information matrix.In particular, where Î−1 is the inverse of I( α, β), with main diagonal elements Î11 and Î22 given by Using (17), large-sample (1 − ϑ)100% confidence intervals for α and β are where Z ϑ/2 is the upper ϑ/2 quantile of the standard normal distribution.

Bayes Inferences
To make Bayes inferences about the parameters α and β, we assume that α and β are independent random variables having gamma prior distributions with parameters (a 1 , a 2 ) and (b 1 , b 2 ), respectively.Using the results in Sarhan et al. (2012), we can get the posterior joint pdf of (α, β), the marginal posterior distributions of α and β, the Bayes point estimates of α and β and the corresponding minimum Bayes risks, and the two-sided Bayesian probability intervals for α and β, by replacing ∑ n i=1 δ i x β i , respectively, and using the same . The posterior density of (α, β) is where Sarhan et al. (2012), all Bayes results have no analytic solutions in the general case.Therefore, numerical approximations and/or simulation techniques should be used.

Approximation to Bayes Estimates
Sarhan et al. ( 2012) used (i) Simulation-based method (SB); (ii) Monte-Carlo integration (MCI) method; and (iii) Lindley approximation method (LA), to make Bayes inferences about the two parameters and they concluded, based on a simulation study, that SB method provides better estimates than the other two methods.In this paper we use rejection sampling (RS) algorithm and Markov Chain Monte Carlo (MCMC) techniques to perform Bayesian inferences about the model parameters and compare it with SB method.

Rejection Sampling (RS)
The rejection sampling technique is used to produce simulated independent samples from a given density function.
For more details about the RS, we refer the readers to Albert (2009) and Gelman et al. (2014).Suppose we want to sample from the density p(θ|data), which is not an easy function to simulate from.The rejection sampling requires a positive function g(θ) defined for all θ for which p(θ|data) > 0 such that: (1) we can draw from a probability density function proportional to g; (2) it is not required that g(θ) be a density function, but must have a finite integral; (3) p(θ|data) < Cg(θ) for all θ, where C > 1 is an appropriate bound on p (θ|data)  g(θ) .The functions p(θ|data) and g(θ) are called the target and proposal functions, respectively.In our case here, the target function is the posterior density (18).The main task in the rejection sampling is to find a suitable proposal density g(θ) and constant value C that satisfy the above restriction.The following is the RS algorithm to simulate draws from (18): 1. Simulate θ = (α, β) from the proposal g(θ).
2. Simulate U from a uniform distribution on the unit interval.
The main advantage of the RS algorithm is the accepted θ has the correct posterior distribution.
The main problem with this process is that C is generally large in high-dimensional spaces and since P(accept) ∝ 1 C , many samples will get rejected.

Markov Chain Monte Carlo (MCMC)
Since early 1990s, Markov Chain Monte Carlo (MCMC) techniques have been used extensively because of their generality and flexibility along with the massive development of computing facilities.It has become one of the main computational tools in modern Bayesian statistical inference.Metropolis et al. (1953) developed a simple version of MCMC, known as the Metropolis-Hastings algorithm (MHA).A generalization of the MHA was proposed by Hastings (1970).A comprehensive theoretical exposition of this algorithm is given by Tierney (1994) and an excellent tutorial of this topic is provided by Chib and Greenberg (1995).
MCMC methods enable quantitative scientists to use highly complicated models and estimate the corresponding posterior distributions with accuracy.For extensive details of the use of MCMC methods, we refer readers to Gilks et al. (1996) and Chen et al. (2000).In this paper, we use the Metropolis-Hastings algorithm (MHA) to generate samples from the posterior distribution π(θ|data) in ( 18).
The MHA can be described by the following iterative steps: The MHA converges to its equilibrium distribution regardless of the proposal distribution q.Neverthless, in practice, choice of the proposal is important since a poor choice considerably delays the convergence towards the equilibrium distribution.An efficient proposal, for the underlying model, is the bivariate normal with mean of the MLE of θ and covariance matrix of the Fisher-information matrix evaluated at the MLE of θ (Lecam, 1986).

Progressive Type-II Censored Data
In this section, we use progressive type-II censored data from TPBT distribution to make inference about the distribution parameters using Bayes technique.In this data, it is assumed that n independent and identical units are placed on a life test simultaneously with m failures to be observed.When the first failure occurs, r 1 survived units are randomly removed from the test.At the second failure, r 2 survived units are randomly removed.The test stops at the occurrence of m-th failure and the remaining r m = n − ∑ m−1 i=1 r i − m survived units are all removed.The likelihood function in this case is The log-likelihood function is ) .
Wu ( 2008) obtained the MLE of the parameters α and β and discussed their confidence intervals and joint confidence region using a progressively type-II censored sample.Wu (2008) addressed that the exact confidence interval and the exact confidence region are obtained, but (i) numerical approximation is used to get the confidence interval of β and no confidence interval of α was discussed, and (ii) the confidence limits for β, which are needed to get the confidence region for (α, β), are obtained numerically.
In this paper we discuss the likelihood intervals for both α and β, the likelihood region for (α, β), and Bayes analysis of the two parameters α and β using progressively type-II censored sample.Using all results in the previous sections, we can get the corresponding results for maximum likelihood (such as the likelihood intervals and likelihood region) and Bayes techniques (point estimate and two-sided Bayesian intervals) using a progressively type-II censored sample by replacing , respectively.Notice that K 0 = ln P x .

Applications
In this section, two applications are discussed and one illustrative example.All the computations were done using Matlab 7.11.0(R2010) programming language, except for RS we used R programming language.In all calculations, the values of M and T for SB and MCMC are 20000 and N = 10000 for RS, also since we do not have any prior information, we assume a 1 = a 2 = b 1 = b 2 = 0.001.Although these values of the prior parameters imply improper priors on α and β, but the corresponding joint posterior is proper.

Application 1: Serum data
In this application, the data set refers to the serum-reversal time (days) of 148 children contaminated with HIV from vertical transmission at the university hospital of the Ribeiro Preto School of Medicine (Hospital das Clnicas da Faculdade de Medicina de Ribeiro Preto) from 1986 to 2001 (Silva, 2004) .For more details, we refer the readers to Perdon (2006).We start with the TPBT(α, β) and W(α, β) distributions for fitting and classical inferences.Then, we will focus on the Bayesian inference of the TPBT(α, β) model.
It is assumed that the lifetimes are independent and identically distributed, and also independent from the censoring mechanism.Also, we assumed that the lifetimes follow either Weibull or TPBT model.In order to assess which model is more appropriate to fit this data set, we plot the scaled TTT-transform (Figure 1.a), the parametric hazard function (Figure 1.b), the empirical and estimated survival functions of the TPBT(α, β) and W(α, β) distributions (Figure 2). Figure 1.a shows that the scaled TTT-plot for the data set has first a convex shape and then a concave shape which indicates a bathtub-shaped hrf. Figure 1.b shows that the TPBT model has a bathtub hrf while Weibul has an increasing hrf.Therefore, the TPBT model would be an appropriate model for the fitting of this data.The MLEs of the model parameters and the corresponding values of the L (the log-likelihood function), AIC (Akaike Information Criterion), K-S (Kolmogorov-Smirnov) test statistic and the corresponding p-value are given in Table 1.These results show that the AIC of the TPBT is smaller than that for Weibull which indicates that the TPBT model is more appropriate to fit this data than the Weibull.Furthermore, based on the K-S test the Weibull model is rejected at any significance level greater than 0.025 while the TPBT is not rejected.To show that there is a unique maximum value for the likelihood function of β, for every model, the plot of the profile log-likelihood function of β, L β (β), for each model is provided in Figure 3. Also, the plot of log-relative Table 2 shows the 14.7% LI and 95% asymptotic CI for the TPBT and Weibull model parameters.The LI and asymptotic CI are not very similar simply because the relative log-likelihood functions of the parameters are not symmetric.The widths of the LI and Asymptotic interval for α in the TPBT model are 17.26×10 −5 and 2.37×10 −5 and for β are 0.0504 and 0.0087 respectively.These indicate that the asymptotic interval is more precise than the LI for both parameters.We used the three Bayesian techniques RS, MCMC and SB to estimate the TPBT parameters.The point estimates and the Bayesian intervals with the corresponding widths are shown in Table 3.The acceptance rates for SB, MCMC and RS are 21.288%, 64.905% and 50.66%, respectively.Figure 4 displays the simulated draws from the joint posterior distribution along with the contours of the likelihood function from the MCMC and SB procedures and the simulated draws using the RS algorithm on the contour plot of the posterior density.These plots show: (1) a good agreement between the likelihood function and the simulated draws from the MCMC; (2) the simulated draws from the SB are not well distributed within the likelihood contours; (3) the draws from the posterior distribution using the RS are well distributed on the contours of the exact posterior density.To compare the performance of the MCMC and SB to RS, we plotted the marginal posterior density functions using the three methods in Figure 5. From Figure 5, one can see that: (1) MCMC gives better estimate of the posterior densities than the SB; (2) SB provides an approximation to the posterior density with very small spread which explains why the interval estimates using SB are narrower than those obtained from RS and MCMC; (3) the posterior distribution of α is right skewed, while that of β is approximately symmetric.3,6,7,11,18,18,18,18,21,32,36,45,47,50,55,60,63,63,67,67,75,79,82,84,84,85,86) and r = (0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0).Wu (2008) assumed that the data have the TPBT(α, β) distribution and obtained the MLE of the two parameters α and β, the 95% CI for only β (as shown in Table ??) and a 95% confidence region for (α, β).The 95% confidence region is (Wu;2008) 0.2147 < β < 0.4772, 46.2248 We obtained the MLE and the corresponding standard error (se), likelihood intervals, likelihood regions and asymptotic confidence intervals for the TPBT distribution parameters.The inverse of the information matrix using the complete data and progressively Type-II censored data, respectively, are ] Also, we computed Bayes estimates and the two sided Bayesian intervals and their widths, using SB, MCMC and RS techniques, for the complete and progressively Type-II Aarset data as shown in Table 4. Figure 6 displays the simulated draws from the joint posterior distribution, using the MCMC procedure, on the contours of the likelihood function and the 14.7% likelihood region of (α, β) and the simulated draws from the posterior density, using RS, on the contour plot of (α, β), based on the complete and progressively Type-II censored Aarset data.We did not include the plot of the simulated draws obtained from SB in Figure 6 because, as we noticed in the first example, it does not provide representative samples.
Based on the results in Table 4 and plots in Figure 6, we can conclude that: (1) Bayesian intervals are more precise than the asymptotic and likelihood intervals; (2) MCMC gives more reasonable approximations to the

Complete data
Progressively censored data Figure 6.The 100p% likelihood contours along with the simulated draws of (α, β) from MCMC (top row) and the simulated draws from the posterior distribution of (α, β) on the contour plot using RS (bottom row) for complete (left column) and progressively censored (right column) Aarset data.

Simulated Data
To compare our proposed technique with technique in Wu (2008), we use the same simulated data in Wu (2008).Wu (2008) generated a progressively Type-II censored sample with n = 15 from the TPBT distribution with α = 0.02 and β = 0.5, with censoring scheme r = (1,0,2,0,0,0,1,0,1,0) and got the observations x = (0.78, 3. 15, 5.15, 6.69, 7.09, 7.40, 14.28, 15.72, 15.92, 22.59).Wu (2008) calculated the 95% CI for β as given in table ??, and the 95% confidence region as 0.2299 < β < 0.7075, 8.5737 Table 5 shows the results obtained from the maximum likelihood and Bayes techniques.Based on the percentage errors, MCMC gives better estimation than both SB and RS algorithms for this simulated data set.The inverse of the information matrix is Figure 7 displays The likelihood contours, the 95% confidence region along with the simulated draws from the posterior distribution using MCMC and SB, contour plot of the parameters (α, β) of the posterior density and of the transformed parameters (ln α, β).The contour plot of (α, β) shows the skewness in the posterior density, especially towards when α gets larger.This is why we consider the transformation (ln α, β).The skewness has been reduced and the distribution becomes more symmetric as shown in the bottom right plot in Figure 7.We also plot the simulated draws from the RS on the contour plot of the log posterior density in Figure 7.As expected, most of the draws fall within the contour of the exact posterior density.

Simulation Studies
To evaluate the performance of maximum likelihood and Bayes procedures (using the SB, MCMC and RS techniques) based on the sample and the censored sizes, a large simulation study using Monte Carlo method is carried out according to the following scheme: 1. Specify the sample size n.
2. Specify the censored size n − m as a percentage of n.
3. Generate a random sample with sizes n and m from TPBT(α, β).
5. Compute the Fisher information matrix at the MLE's of α and β.
6. Compute Bayes estimates of α and β using SB, MCMC, and RS methods.
7. Compute the squared deviation of the point estimate of each parameter using each of the four procedures from the corresponding true value.
8. Repeat Steps 3-7 1000 times.9. Calculate the average of the point estimates of every parameter and the mean squared error (MSE) associated with each estimate for the two parameters using the four techniques.12. Steps 1-11 are carried out when α = 0.1 and β = 0.5.
Table 6 presents the average estimates (first row) and the corresponding MSE (second row) corresponding to every parameter at different set of values of (n, m).
Based on the simulation results, one can conclude that: (1) the MSE decreases when n increases for all cases at every level of censoring, (2) the MSE increases when the percentage of censorship increases mainly for large n, (3) there is no significant difference in the values of the MSE corresponding to the three techniques.However, based on the MSE, (i) the MCMC generally provides better approximations than SB; (ii) RS produces 75.62% better estimations than MCMC.

Conclusion
This paper extends the work of Sarhan et al. (2012) which considered a two parameter bathtub shaped distribution that was revisited by Chen (2000).We provide statistical inference of the two parameters using maximum likelihood and Bayesian methods.For Bayesian, the simulation based, Markov Chain Monte Carlo and rejection sampling techniques are applied.The specific interest is in the MCMC and RS methods.The estimation techniques were applied to two real data sets and a simulated data set.In addition, the estimation methods were compared by a Monte Carlo simulation.
For the two real data sets, the Bayesian estimates of the parameters from the MCMC and RS methods were not very similar to those obtained by the SB method.The widths of the Bayesian intervals for the parameters constructed by the MCMC and RS methods were larger.In each parameter case, the widths of the MCMC and RS intervals are about 5 times larger than the SB interval widths.The main reason was that the updates from the SB algorithm are not representative.
Several phenomena can be observed from the results from the Monte Carlo simulation.First the mean squared error decreases as the sample size increases for fixed level of the number of units censored.For a fixed value of the sample size beyond some threshold value, the mean squared error increases with increasing level of units censored.
The RS method gives better approximations of the point estimates of the parameters as the corresponding mean squared errors are generally smaller than those from the other methods.
Figure 1.a) The Scaled TTT-Transform and b) the parametric hazard functions.

Figure 2 .Figure 3 .
Figure 2. The empirical and parametric survival functions of the serum.

Figure 4 .
Figure 4.The 100p% likelihood contours along with the simulated draws of (α, β) from MCMC (top left) and SB (top right) and the simulated draws from the posterior density on contour plot of (α, β) using RS (bottom) for serum data.

Figure 5 .
Figure 5.The marginal posterior density function using the three methods for serum data.

Figure 7 .
Figure7.The 100p% likelihood contours, 95% confidence region along with the simulated draws of (α, β) using MCMC (top left) and SB (top right), contour plot of (α, β) (bottom left) and of (log α, β) (bottom right) along with simulated draws from the posterior density using RS for the progressively censored simulated data.

Table 1 .
MLEs of the parameters, the corresponding L, AIC, K-S and p-value.

Table 3 .
Bayesian results for the TPBT model.

Table 4 .
MLE and Bayes results for the Aarset data.

Table 5 .
Point estimates and the corresponding percentage errors and interval estimates of the parameters for the simulated data.

Table 6 .
The mean estimated values and corresponding MSE using MLE and Bayes technique methods.