Estimation of P ( Y < X ) for a Two-parameter Bathtub Shaped Failure Rate Distribution

This paper deals with the estimation of reliability R = P[Y < X] when X and Y are two independent random variables with a two-parameter bathtub shaped failure rate distribution with the same second shape parameter. Likelihood and Bayesian methods are proposed to make inferences about R. We obtain the likelihood interval and asymptotic confidence interval for R, and we consider Bayesian point estimates of R under both absolute and squared error loss, using either gamma or uniform priors for the three unknown model parameters. An equal tail Bayesian credible interval for R is investigated. Analysis of a real data set is presented for illustrative purposes, and Monte Carlo simulations are performed to compare: (1) the performance of Bayes estimates under two different loss functions; and (2) the maximum likelihood and Bayesian methods.

The probability density function (pdf) is , x > 0.
Here α and λ are both shape parameters.Henceforth, we will denote this distribution by TPBT(λ, α) or TPBT.
Due to the convenient structure of the TPBT distribution, it can be used in analyzing many lifetime data sets.
The failure rate function takes a bathtub shape when α < 1 and is increasing if α ≥ 1, Chen (2000).The TPBT distribution is the only known two parameter distribution with a bathtub shaped failure rate function that has exact joint confidence regions for the parameters Chen (2000).Sarhan et al. (2012) discussed Bayes estimation of the two parameters of the TPBT distribution.
In a stress-strength model, the stress Y and the strength X are treated as random variables and the reliability, R, of a component during a given period is taken to be the probability that its strength exceeds the stress during the entire interval.Due to its practical importance, the estimation of R = P(Y < X) has attracted the attention of many authors.The maximum likelihood estimate (MLE) of R, under the assumption that X and Y are independent and normally distributed, is derived by Church and Harris (1970).The uniformly minimum variance unbiased estimate of R under the same assumption is obtained by Downtown (1973).Recently, Kundu and Gupta (2005) have considered estimation of P(Y < X) when X and Y are independent generalized exponential random variables.
The Bayesian approach requires integration over the posterior distribution and this is often impossible to carry out analytically.In such cases, MCMC techniques can be applied to draw samples from the posterior distribution.
The MCMC methodology provides a convenient and efficient method to sample complex, highly-dimensional distributions.
The rest of the paper is organized as follows.Section 2 presents the MLE of the unknown parameters and the reliability.The asymptotic confidence intervals of the unknown parameters and reliability are discussed in Section 3. Sections 4 discusses Bayesian inference.A complete analysis of a real data set is given in Section 5 for illustrative purposes.Section 6 gives a simulation study, and Section 7 concludes the paper.

Maximum Likelihood Estimate of R
Assume that there are two independent random variables Y and X such that X ∼ TPBT(λ, α) and Y ∼ TPBT(γ, α). Therefore, To obtain the MLE of R, we first derive the MLE of the two parameters γ and λ and then use the invariance property.Now, suppose ′ denote a vector of the three unknown parameters.The likelihood function, . (2) Thus, the log-likelihood function, The likelihood equations are From ( 5) and ( 6), we obtain the MLE of λ and γ as functions of α, say λ(α) and γ(α) respectively, as Substituting λ(α) and γ(α) into (3), we get the profile log-likelihood function for α as Therefore, the MLE of α, say α, can be obtained by maximizing (9) with respect to α.It can be shown that the α can be obtained as a solution of a non-linear equation of the form where A very simple and effective iterative procedure α ( j+1) = h(α ( j) ) can be used, where α ( j) is the jth iterate of α.The iteration procedure should be stopped when α ( j+1) − α ( j) is sufficiently small.Once we obtain α, λ and γ can be calculated from ( 7) and ( 8) respectively, and the MLE of R is R = γ γ + λ .

Confidence intervals for R
Note that, α, λ and γ are not in explicit form.Further, it is not possible to obtain the variances θ = ( α, λ , γ) ′ in an explicit form.We propose to use two approximate confidence intervals for R.

Likelihood Interval
The maximum log-relative likelihood function for R, r max (R), is and α(R) is the solution of the following equation in α for a given R A 100p% likelihood interval (LI) for R is the set of all possible real solutions, in R, of the following inequality This inequality has no closed solution in R, so the interval must be obtained numerically.The 100p% LI approximates a (1 − ϑ)100% CI when p = exp{− 1 2 χ 2 1 (ϑ)}, where χ 2 1 (ϑ) is the upper ϑ quantile of χ 2 1 .

Asymptotic Confidence Interval
The MLE of the vector of unknown parameters θ = (α, λ, γ) ′ is asymptotically normally distributed with mean of true θ and variance-covariance matrix I −1 , and is the inverse of the expected information matrix I(θ).I(θ) is consistently estimated by Î = I( θ), which is given in the Appendix.Let Î j j be the ( j, j)-entry of Î, j = 1, 2, 3, then the asymptotic (1 − ϑ)100% confidence intervals for θ j , j = 1, 2, 3, are where Z ϑ/2 is the upper ϑ/2 quantile of the standard normal distribution.
One can show that, R is asymptotically normally distributed with mean of R and variance Therefore, an asymptotic (1

Bayesian Inference
In this section, we discuss Bayesian methods for making inferences about R = γ λ+γ .We assume that the prior pdfs of the elements of θ = (α, λ, γ) ′ are independent, and that they are either gamma distributed π(u) ∝ α a−1 e −b u , u > 0 for fixed values of a, b > 0, or that they are they are uniformly distributed on (0, B), for some large B.
Using the likelihood function (2) and independent uniform priors, the joint posterior pdf of θ = (α, λ, γ) ′ is proportional to the likelihood function, while using independent gamma priors the joint posterior pdf θ is In either case the joint posterior density has a complicated form and it is unlikely that closed form inferences for the parameters α, β, γ, or for the reliability R, are possible.As indicated by Gilks et al (1996), any feature of the posterior distribution is a legitimate candidate for Bayesian inference, for example, moments, quantiles, or highest posterior density regions.Such quantities can often be expressed in terms of posterior expectations of functions of θ, and while the integrations in such expectations are typically difficult to be evaluate analytically, Markov Chain Monte Carlo (MCMC) sampling methods can often be used to approximate the integrals.In any event, Bayesian inferences are always based on the posterior distribution, and where analytical forms are not not available, inference is most often accomplished using samples from the posterior distribution.Good general references are Aykroyd (1998) and Gelman et al (2009).
We sampled from the posterior distribution p(θ|data) using the Metropolis-Hastings algorithm (Hastings, 1970), summarized as follows.
1. Choose an intitial value
Here q is the transition probability matrix of a Markov chain whose support is the same as that of the likelihood function, and q(θ * |θ) is the probability of a transition from θ to θ * .Tierney (1994) provides an extensive discussion on the use of the Metropolis-Hastings and related algorithms for sampling posterior distributions, including the choice of the proposal distribution.Under quite general conditions the Metropolis-Hastings algorithm generates a sample path {θ (t) , t = 0, 1, . . .T } from an ergodic Markov Chain whose marginal distributions converge to the equilibrium distribution of the chain -in this case the posterior distribution from which we wish to sample -as T → ∞.While there are many possible choices for the proposal distribution, in practice the choice of the proposal is important since a poor choice can delay the convergence towards the equilibrium distribution.
We generated proposals using an independence chain with proposals chosen as i.i.d.samples from N 3 ( θ, Î−1 ).The Bernstein-von Mises's theorem (Lecam, 1986) states that this is the limit distribution of the posterior p(θ|data), as n → ∞.The idea of using an independence chain in which the proposed point is independent of the past and current states was suggested by Hastings (1970), among other proposal strategies.

Application
In this section we present a data analysis of the strength data reported by Badar and Priest [4].The data represent the strength, measured in GPa, for single carbon fibers and impregnated 1000-carbon fiber tows.Single fibers were tested under tension at gauge lengths of 1; 10; 20; and 50mm.Impregnated tows of 1000 fibers were tested at gauge lengths of 20; 50; 150 and 300 mm.For illustrative purposes, we consider the single fibers of 20 mm (Data Set I) and 10 mm (Data Set II) in gauge length, with sample sizes n = 69 and m = 63, respectively.The data are presented below for convenience.We fit the TPBT model to the two data sets separately, and present the estimated parameters of both generalized exponential (GE) and TPBT distributions, log-likelihood values L, Kolmogorov-Smirnov (K-S) distances and corresponding p-values in Table 1.Also, the Akaike information criterion (AIC= −2L + 2k, k is the number of model parameters) is computed for every model using the two data sets.Based on either the P-value or the AIC, the TPBT model fits quite well (and better than the GE model) to both the data sets.Empirical and fitted survivor functions are shown in Figure 1.Assuming a common parameter α for the TPBT model, we obtained the MLEs α = 1.2015, λ = 0.042308 and γ = 0.012396 leading to the MLE of R as R = 0.22661.The inverse of the observed information matrix of ( α, λ, γ) is Using ( 12) and ( 15), we obtained ÎR = 1.2638 × 10 −3 .Therefore, using ( 11) and ( 13), we obtained the 14.7% LI (which approximates a 95% CI) and the asymptotic 95% CI for R as (0.17984, 0.28367) and (0.15693, 0.29628) respectively.The intervals are close, but not identical, reflecting the asymmetry of the maximum log-relative likelihood function r max (R), which is plotted in Figure 2 together with the endpoints of the 14.7% LI for R.   Figure 3a shows values of α, λ and γ sampled from the joint posterior distribution assuming gamma(.001,.001)priors, after deleting an initial transient of 10,000 points and subsampling each 200'th point.The empirical bivariate marginal distributions are also included in Figures 3b, c and d along with p = 0.01, 01.47 and 0.50 profile likelihood contours.There is a good agreement between the shape of the contours and the sample from the posterior indicating little dependence on the prior.
Figure 4 shows the 1000 sampled points from the marginal posterior distributions of α, λ, γ and R calculated using gamma priors, and histograms of the empirical marginal posterior pdfs.The dashed lines on the left hand panels are positioned at the upper and lower limits of the 14.7% likelihood intervals.In order to assess the sensitivity of Bayesian inferences for R to the choice of prior we ran the Metropolis-Hastings algorithm twice: with independent U(0,100) priors, and with independent gamma(0.001,0.001)priors.The uniform distribution was meant to play the role of a non-informative prior, while the gamma prior is relatively informative, being heavily skewed to the right with a single high mode near 0. We ran the sampler for 100,000 iterations and discarded the first 50,000 points as a transient.The remaining 50000 points were sub-sampled, retaining every 50th point, and the resulting sequences of 1000 sampled values of R are shown in Figure 5, together with plots of the associated simulated empirical posterior pdf.These figures suggest that the choice of prior has little influence on inferences for R, and this is further evidenced in Table 2, which shows the estimated posterior means and medians.
Equal tail probability 0.95 credible intervals for R were calculated as the interval from the 2.5'th to the 97.5'th percentile of the empirical distribution of R. The intervals, (0.177, 0.276) using the gamma prior, and (0.180, 0.277) using the uniform prior, again show little sensitivity to the choice of prior.

Simulation Study
A simulation study was carried out to assess the sampling properties of the MLE, the posterior mean and the posterior median.
At each of several parameter values and several sample sizes, 1000 simulation batches were independently generated.For each simulation batch, X 1 , . . ., X n and Y 1 , . . ., Y m were independently drawn from the TPBT(λ, α) and TPBT(γ, α), respectively.Three values were used for each parameter.These were approximately the MLE and the MLE plus or minus 2 standard errors, based on the calculations from the data analysis in Section 5.
For each simulation batch, the MLE and asymptotic 95% confidence interval (13) for R were evaluated.In addition, for each batch the Metropolis-Hastings algorithm was used to generate a sequence of 10000 samples from the posterior distribution, using independent U(0,100) priors.The data analysis of Section 4 showed little sensitivity to the choice of prior, and we used the uniform as a non-informative choice.The first 5000 iterates of each sequence were discarded and the remaining 5000 were assumed to be samples from the stationary posterior distribution of (α, λ, γ) ′ .These were transformed to give samples from the posterior distribution of R.
Based on the 1000 simulation batches, estimates of the bias and root mean squared error were made for the MLE, the posterior mean and the posterior median.The empirical coverage probabilities were estimated for the asymptotic 95% confidence interval and the equal tail probability .95credible interval, the latter being the interval between the 2.5'th and 97.5'th percentiles of the empirical posterior distribution.The results are presented in Table 3.
The Bernstein-von Mises' theorem guarantees that for any choice of prior, the MLE, the posterior mean and the posterior median will be close for large sample sizes.In the present case there are only small differences among the MLE, the posterior mean and median with respect to bias and mean squared error, even at small to moderate sample sizes.The coverage property of the frequentist intervals is close to nominal, even when m = n = 20.

Conclusion
For the situation when two independent random variables follow the TPBT distribution with equal second shape parameter, the reliability has a particularly simple form, although it still must be estimated numerically.We have examined the use of both likelihood and Bayesian inference for the reliability in this case.Uniform and gamma priors were proposed, and MCMC methods were used to examine the posterior density.Point estimates via maximum likelihood, posterior mean and posterior median all displayed excellent sampling properties.Asymptotic confidence intervals gave coverage close to the nominal level.
There is no guarantee that Bayesian inferences will have frequentist validity.This is the case in Table 3 which shows that the probability .95Bayesian credible intervals for R have frequentist coverage below 95%.There have been systematic attempts to develop families of prior distributions, so-called probability matching priors, having frequentist validity (Datta, 1996).
A typical adjunct to a Bayesian analysis involving samples generated with MCMC algorithms is the assessment of the stationarity of the sampled points.A number of methods for this purpose are discussed, for example, in Gelman et al (2009).Using several such procedures, there was no indication from the data analysis that more than 10000 points were needed for convergence, using either uniform or gamma priors.Due to the size of the simulation study, we made no assessment of convergence of the sampled points for the individual simulated data sets.

Copyrights
Copyright for this article is retained by the author(s), with first publication rights granted to the journal.
This is an open-access article distributed under the terms and conditions of the Creative Commons Attribution license (http://creativecommons.org/licenses/by/3.0/).
Empirical and fitted survival functions for the two data sets using GE and TPBT models.

Figure 2 .
Figure 2. Max Log-relative likelihood function of R along with the endpoints of 14.7% LI for R.

Figure 5 .
Figure 5. Simulated points and empirical posterior pdfs of R using gamma (top row) and uniform (bottom row)priors.

.
The local estimate of the variance-covariance matrix for the MLE of the model parameters is the inverse of the observed

Table 1 :
MLE of the parameters, K-S, P-value, L, and AIC for the two data sets.