The Exponentiated Poisson-Exponential Distribution: A Distribution with Increasing, Decreasing and Bathtub Failure Rate

In this paper we propose a new three-parameters lifetime distribution with increasing, decreasing and bathtub failure rate dependingonitsparameters.Theproperties oftheproposed distributionarediscussed,includingexplicitalgebraicformulaefor its reliability and failure rate functions, quantiles and moments. A simulation study is performed in order to verify the behavior of the maximum likelihood estimates. The methodology is illustrated on a real data set.


INTRODUCTION
The exponential distribution (ED) provides an elegant, simple and close form solution to many problems in lifetime testing and reliability studies. However, the ED does not provide a reasonable parametric fit for some practical applications where the underlying failure rates are nonconstant, presenting monotone shapes. In recent years, in order to overcome such problem, new classes of models were introduced based on variations of the ED. For instance, Gupta and Kundu [1] proposed a generalized ED, which can accommodate data with increasing and decreasing failure rate function. Kus [2] proposed another modification of the ED with decreasing failure rate function, while Barreto-Souza and Cribari-Neto [3] generalizes the distribution proposed by [2] by including a power parameter in his distribution.
Recently, Cancho et al. [4] introduced the Poisson-Exponential (PE) distribution whose cumulative distribution function is given by where > 0 and > 0. This distribution family based on the ED has increasing failure rate function.
As described by Marshall and Olkin [5], an exponentiated distribution can be easily constructed. It is based on the observation that by raising any baseline cumulative distribution function (cdf) F baseline (y) to an arbitrary power > 0, a new cdf F(y) = [F baseline (y)] > 0 is obtained, but now with the additional parameter , which can be refereed as a resilience parameter and F(y) is a resilience parameter family. Although it is not our case, the term resilience easily emerges if we let be an integer. In this case, F(y) can be seem as the cdf of a parallel system with independent components, which is less likely to failure as the number of components increases, leading to a resilient structure. Following this idea, several authors have considered extensions from usual survival distributions. For instance, Mudholkar et al. [6] considered the exponentiated Weibull distribution as a generalization of the Weibull distribution, [1] introduced the exponentiated ED as a generalization of the usual ED, Nadarajah and Kotz [7] proposed exponentiated type distributions extending the Fréchet, Gamma, Gumbel and Weibull distributions and Ristic and Nadarajah [8] introduced the exponentiated exponential-Poisson distribution as a generalization of the exponential-Poisson distribution due to [2]. Hence, the cdf of the exponentiated PE (EPE) distribution is defined from (1) as where > 0, > 0 and > 0.
The paper is organized as follows: In Section 2 we present the probability density function, survival function and hazard function of the new distribution, the EPE distribution, and present graphics of such functions for some parameters values. Some properties of the EPE distribution are also presented in Section 2. In Section 3 we present the inferential procedure for the model parameters. Section 4 contains the results of a Monte Carlo simulation on the finite sample behavior of maximum likelihood estimates. In Section 5 the methodology is illustrated on a real data set. Some final comments in Section 6 conclude the paper.

THE EPE DISTRIBUTION
The probability density function (pdf) of the EPE distribution corresponding to the cdf (2) is given by where > 0, > 0 and > 0. When = 1 the model (3) reduces to the PE distribution proposed by [4], and when = 1 and letting → 0 + the model (3) reduces to the ED with parameter .
If u and v are real numbers with |u| > |v| and > 0, we have the series representation If is an integer, the sum in (4) stops at − 1. Hence now e − e − y > e − , we expand (e − e − y − e − ) −1 as in Eq. (4) and then the pdf in (3) can be reduced to where f PE (y; (k + 1), ) is the pdf of a random variable with PE distribution with parameters (k + 1) and . We then have the EPE pdf can be written as a linear combination of PE pdfs. Figure 1 (left panel) shows the pdf (3) for some values of parameters.

Failure Rate
From (3) and (6) it is easy to verify that the failure rate function is given by The EPE distribution accommodates increasing, decreasing and bathtub failure rates depending on the parameter values. Define the function ] .
Considering the inequality in equation above implies that if ≥ 1 then ′ (y) > 0, ∀y > 0. Thus, from Glaser's theorem (see [9]), it follows that the failure rate function is increasing.
If ≤ 1 and < 1 we have ( −1) Hence, from Glaser [9], we conclude that the failure rate function is decreasing. We also observed that the failure rate function may be bathtub. The distinct types of failure rate shapes are illustrated in Figure 2 for some different parameter combinations of the EPE distribution.

Some Properties
Many of the most important features and characteristics of a distribution can be studied through its moments, such as mean, variance, tending, dispersion, skewness and kurtosis. We shall now provide a general expression for the rth ordinary moment of the EPE distribution. Cancho et al. [4] showed that the rth moment of the random variable Y with PE distribution with parameters and can be expressed as where is the generalized hypergeometric function, with a = [a 1 , … , a p ], p is the number of operands of a, b = [b 1 , … , b q ] and q is the number of operands of b.
We then formulate the following result. Proposition 2.1. Let Y be a random variable with EPE distribution and parameters , and . If > 0 non-integer, using (5) and (7) we obtain the raw rth ordinary moment of Y can be written as Proof. The proof is obtained by direct integration and it is then omitted.
Order statistics are among the most important tools in nonparametric statistics and inference. Let Y 1 , … , Y n be a random sample taken from the EPE distribution and let Y 1∶n , … , Y n∶n denote the corresponding order statistics. The probability density function f i∶n (y) of the ith order statistic Y i∶n is given by According to Barakat and Abdelkader [10], the rth moment of the ith order statistic Y i∶n can be represented as Consider the binomial series expansion given by where (r) k = r(r − 1)...(r − k + 1) is a Pochhammer symbol. If |x| < 1 the series converges and We then state the following.

Proposition 2.2. For the random variable Y with EPE distribution, we have that the rth moment of the ith order statistic is given by
Proof. From (6) and (8), and using (9), (10) and the Taylor series of the exponential function This completes the proof.
Given that a component survives up to time t ≥ 0, i.e., there was no failure prior to time t, the residual life is the period beyond t until the time of failure and defined by the conditional random variable Y − t|Y > t. In reliability, it is well-known that the mean residual life function and the ratio of two consecutive moments of residual life characterize the distribution uniquely (see [11]). Thus, we obtain the rth order moment of the residual life via the general formula Proposition 2.3. For the random variable Y with EPE distribution, we have that the rth order moment of the residual life is given by Proof. From (11), using (6), (9), (10) and the Taylor series of the exponential function e −x , and proceeding a similar way in the proof of Proposition 2.2 follow the result.
Then, we obtain the mean residual life of the EPE distribution as On the other hand, we analogously consider the reversed residual life and some of its properties. The reversed residual life can be defined as the conditional random variable t − Y|Y ≤ t, which indicates the time elapsed from the failure of a component given that its life is less than or equal to t. This random variable may also be called the inactivity time or time since failure (for more details, you may see [12] and [13]). Also, in reliability, the mean reversed residual life function and the ratio of two consecutive moments of reversed residual life determine the distribution uniquely.
The rth order moment of the reversed residual life can be obtained via the well-known formula where WhittakerM ( , , z) is the Whittaker M function, which can be defined in terms of the hypergeometric function as follows: Proof. From (12) and using F ( y ) given by (2), we have that dy. Now using (9) and proceeding a similar way in the proofs of Propositions 2.2 and 2.3 follow the result.
Thus, the mean of the reversed residual life of the EPE distribution is given by ] .
The Bonferroni and Lorenz curves (see [14]) and Gini index have many applications not only in economics to study income and poverty, but also in other fields like reliability, insurance and medicine. The Bonferroni curve B F [F(y)] is given by or equivalently by where p = F(y) and F −1 (t) = inf {y ∶ F(y) ≥ t}.

Proposition 2.5. From the relationship between the Bonferroni curve and the mean residual lifetime given by Theorem 2.1 of Pundir et al. [15], the Bonferroni curve of the distribution function F of EPE distribution is given by
Pdf_Folio:280 Proof. From (13) and using f ( y ) given by (3) and F ( y ) given by (2), we have that Now using (9) and making binomial expansion a similar way in the proofs of Propositions 2.2 2.4 follow the result.

Also, the Lorenz curve of F that follows the EPE distribution can be obtained via the expression L F [F(y)] = B F [F(y)]F(y).
The scaled total time and cumulative total time on test transform of a distribution function F are given by respectively (see [15]). If F(y) is the EPE distribution function specified by (2), then using formula (9) and the Taylor series of the exponential function e −x , Finally, the Gini index can be obtained from the relationship G = 1 − C F .

INFERENCE
In this section we discuss inference and interval estimation for the model parameters.

Estimation by Maximum Likelihood
Let y = ( y 1 , … , y n ) be a random sample of the EPE distribution with unknown parameter vector = ( , , ) . The logarithm of likelihood function for is given by The maximum likelihood estimates (MLEs) of = ( , , ) can be directly obtained by maximizing the log-likelihood function (14), or alternatively, by finding the solution for the following three nonlinear equations: The Eq. (17) can be solved exactly for , i.e.,̂= conditional upon the values of̂and, where,̂and̂are the MLEs for , and , respectively. Therefore, it is sufficient to solve the Eqs. (15) and (16) iteratively in order to find the MLEs for and .

Interval Estimation
For EPE distribution the Fisher information matrix is given by where , = The above expressions depend on some expectations that can be easily computed using numerical integration.
Large sample inference for the model parameters can be based, in principle, on the MLEs and their estimated standard errors. Following Cox and Hinkley [16], under suitable regularity conditions, it can be shown that as n → ∞. This asymptotic behavior remains valid if Fisher information is evaluated at the MLEs,, say −1 () . Alternative estimates can be obtained from the inverse observed information matrix −1 () since it is a consistent estimator of −1 ( ) . We can use the approximate trivariate normal N 3 ( , −1 ()) distribution of̂to construct approximate confidence regions for some parameters and survival function of EPE model.
For testing goodness of fit of the EPE distribution and for comparing this distribution with some of its special cases, we consider the likelihood ratio statistics (LRS). For instance, for comparing the EPE distribution with the PE distribution proposed by [4] for a given data set we consider the hypothesis test given by H 0 ∶ = 1 versus H 1 ∶ ≠ 1. In this case, the LRS is given by Λ = 2 {ℓ (̂,̂,̂) − ℓ (̃,̃, 1 ) }, where,̂and̂are the MLEs without restriction and̃and̃are the estimates under H 0 . Following [16], under H 0 we assume an asymptotic chi-square distribution with one degree of freedom for Λ.

SIMULATION STUDY
This section presents the results of a simulation study that has been performed in order to verify the behavior of the standard errors, bias and root mean squared errors of the maximum likelihood estimators of , and . We generated 1, 000 samples of size n = 30, 50, 100, 200 and 500 from the EPE distribution for each one of the three sets of values of = ( , , ): (0.5, 2.0, 0.5), which corresponds to a decreasing failure rate function; (0.5, 2.0, 2.0), from which follows that the failure rate is increasing; and (2.0, 2.0, 0.5), whose failure rate function has a bathtub shape. The results are condensate in Table 1, which shows the averages of the 1, 000 MLEs, Av () , together with their standard errors, SE () , bias, B () , and root mean squared errors, RMSE () . Besides, Table 1 shows the coverage probability of the 95% confidence intervals for parameters of the EPE distribution, CP ( ). Table 1 indicates the following results: convergence has been achieved in all cases, i.e., for all three sets of parameter values; the biases are always smaller than one empirical SE. These results suggest the MLEs have performed consistently. The standard errors of the MLEs decrease when sample size increases. The empirical coverage probabilities are close to the nominal coverage level, 0.95. In particular, as sample size increases.

APPLICATION
In this section we consider a real data set of 520 failure times (10, 000's of hours) of industrial devices. Initially, in order to identify the shape of a lifetime data failure rate function we shall consider a graphical method based on the TTT plot (see [17]). In its empirical version, the TTT plot is given by , where r = 1, … , n and Y i∶n denote the order statistics of the sample. It has been shown that the failure rate function is decreasing (increasing) if the TTT plot is convex (concave). Even though the TTT plot is only a sufficient condition, not a necessary one for indicating the failure rate function shape, it is used here as a crude indicative of its shape. Figure 3 (upper-left panel) shows the TTT plot for the considered data, which is first convex and then concave. It indicates a bathtub-shaped failure rate function, which can be properly accommodated by a EPE distribution.
Then, the EPE distribution was fitted to the data. Table 2 presents the MLEs and the corresponding 95% confidence intervals of the EPE distribution parameters, which were based on the Fisher information matrix at the MLEs, given by Table 1 The averages of the 1, 000 maximum likelihood estimates (MLEs), Av () , their standard errors, SE () , bias, B () , root mean squared errors, RMSE () , and the coverage probabilities of the 95% confidence intervals for parameters of the exponentiated Poisson-Exponential (EPE) distribution, CP( ). The parameter vector is given by = ( , , ). ] .
As discussed in Section 3.2, for comparison of nested models, which is the case when comparing the EPE distribution with the PE distribution and the PE distribution with the ED, we can compute the maximum values of the unrestricted and restricted log-likelihoods to obtain the LRS. Thus, the LRS used for testing H 0 ∶ = 1 results in Λ = 112.211 with p-value < 0.0001, indicating that the EPE model presents a much better fitting than the PE model to the data set under consideration. The LRS used for testing H 0 ∶ = 0 (for details on this comparison, see results in [4]) is Λ ′ = 1.741 with p-value = 0.187, which gives favorable evidence to the restricted model, i.e., the ED, over the full model, i.e., the PE model. The last column of Table 2 presents the log-likelihood values for the three fitted models. We also compare the EPE distribution, PE distribution and ED by inspection of the Akaike's information criterion (AIC), −2ℓ () + 2k, and Schwarz's Bayesian information criterion (BIC), −2ℓ () + k log(n), where k is the number of parameters in the model and n is the sample size. The favorite model is the one with the smaller value on each criterion. The estimated statistics AIC and BIC for the EPE distribution are equal to −472.468 and −459.706, respectively. While, the estimated statistics AIC and BIC for the PE distribution are equal to −362.256 and −353.749, respectively; and the estimated statistics AIC and BIC for the ED are equal to −362.516 and −358.262, respectively. The EPE distribution overcomes the corresponding PE distribution and ED in both considered criterion. These results are corroborated by the plot in the lower-left panel of Figure 3, which compares the estimated survival functions obtained via ED, PE and EPE model fitting on the empirical Kaplan Meier survival function. Figure 3 also shows the estimated densities (upper-right panel) and the hazard function obtained via EPE model fitting (lower-right panel).
Besides, for the sake of illustration, in order to compare our EPE distribution with other lifetime distributions capable of modeling a bathtubshaped failure rate function, we fitted the exponentiated Nadarajah-Haghighi (ENH) distribution developed by Lemonte [18]

CONCLUDING REMARKS
In this paper we propose a new lifetime distribution. The EPE distribution is a generalization of the PE distribution proposed by [4], which accommodates increasing, decreasing and bathtub failure rate functions. It arises on a latent complementary competing risks scenario, where the lifetime associated with a particular risk is not observable but only the maximum lifetime value among all risks. We provide a mathematical treatment of this new distribution including explicit algebraic formulaes for its survival function and failure rate function, quantiles and moments. We discussed maximum likelihood estimation, which presented an adequate numerical stability in the simulation study performed. The EPE distribution allows a straightforwardly nested hypothesis testing procedure for comparison with its PE distribution particular case. The practical relevance and applicability of the EPE distribution were demonstrated in a real data set.

CONFLICT OF INTEREST
No conflict of interest has been reported by authors.

AUTHORS' CONTRIBUTIONS
(i) Vicente G. Cancho wrote the initial drafts of the paper and R codes; (ii) Paulo H. Ferreira helped in writing and improving the manuscript and the R programs, and also did the simulation study and analysis; and (iii) Francisco Louzada helped in the paper writing and supervised overall work.