A New Compound Gamma and Lindley Distribution with Application to Failure Data

In this paper, we propose a new lifetime distribution by compounding the gamma and Lindley distributions. Construction of it can be interpreted in the viewpoint of the reliability analysis and Bayesian inference. Moreover, the distribution has decreasing and unimodal hazard rate shapes. Several properties of the distribution are obtained, involving characteristics of the (reverse) hazard rate function, quantiles, moments, extreme order statistics and some stochastic order relations. Estimating the distribution parameters is discussed by some estimation methods and their performance is evaluated by a simulation study. Also, estimation of the stress-strength parameter is investigated. Usefulness of the distribution among other models is illustrated by fitting two failure data sets and using some goodness-of-fit measures.


Introduction
Different nature of the lifetime data requires introducing new distributions with various failure rate shapes to capture the analysis of such data. Among those shapes, we mention: the decreasing and unimodal (upside-down bathtub) hazard (failure) rates. The former rates have many applications in reliability and survival analysis. Further, a decreasing failure rate describes a phenomenon where the probability of an event in a fixed time interval in the future decreases over time. A decreasing failure rate can describe a period of infant mortality where earlier failures are eliminated or corrected (Finkelstein 2008). A practical example is failure in the air conditioning systems (Proschan 1963). On the other side, as mentioned by (Lai and Xie 2006), when the main reasons of the failures of products are caused by fatigue and corrosion, the failure rates of those products exhibit unimodal shapes. Moreover, in some medical situations, for example breast cancer and infection with some new viruses, the hazard rate is unimodal shape, see (Demicheli, Bonadonna, Hrushesky, Retsky, and Valagussa 2004). and having the probability density function (pdf) f (x; θ, β) = θ 2 ((β + βθ − θ)x + 1) e −θx β(1 + θ) ; x > 0, θ > 0, β > 0. Nedjar and Zeghdoudi (Nedjar and Zeghdoudi 2016) discussed various statistical properties and simulation of the GamL distribution. Such GamL distribution has some weaknesses, where there is an incorrectness in the parameter space. The above f (x; θ, β) is not a proper pdf, since it can be negative for some values of the parameters θ and β (see (Messaadia and Zeghdoudi 2017)). For a proper f (x; θ, β), (Messaadia and Zeghdoudi 2017) modified the parameter space to be θ > 0; β > θ 1+θ . Moreover, when α = β(1 + θ) − θ, then the GamL distribution gives the two parameter Lindley distribution proposed by (Shanker, Sharma, and Shanker 2013), which has the following the pdf f (x; θ, α) = θ 2 (αx + 1)e −θx α + θ ; x > 0, θ > 0, α > 0.
Therefore, the GamL distribution is only a reparametrization for the distribution above proposed by (Shanker et al. 2013). On the other side, the GamL distribution proposed by (Messaadia and Zeghdoudi 2017) does not allow a decreasing or unimodal hazard rate shapes, which makes it non-suitable for modeling data with those hazard rate shapes. Other generalizations of the Lindley distribution based on the same mixture approach can be found in (Abouammoh, Alshangiti, and Ragab 2015) and references therein.
The aim of this paper is to propose a new gamma Lindley distribution by compounding the gamma and Lindley distributions in a different manner than used in the GamL distribution.
The new distribution has decreasing and unimodal hazard rates shapes and its construction is as follows. Let X | λ follow the gamma distribution with pdf f (x | λ) = λ α x α−1 e −λx Γ (α) , x > 0; α, λ > 0, and λ | β having the Lindley distribution (Lindley 1958) with the pdf f (λ | β) = β 2 1 + β (1 + λ) e −βλ , λ > 0; β > 0, then the marginal distribution of X is called gamma-Lindley (GaL) distribution. The pdf of X is obtained by (1 + λ) λ α e −(β+x)λ dλ, which after some algebra, we get the GaL pdf as Moreover, the cumulative distribution function (cdf) of the GaL distribution is hence, the corresponding reliability (survival) function is given by The usefulness of model (1) in reliability analysis comes in noting that X can be the lifetime of a component and λ is the scale parameter of its distribution. Suppose that, in the population, there is some variability in the scale parameters; this variability can be described by the distribution for λ. Also, comparing the GaL distribution with gamma and Lindley distributions implies the flexibility of GaL in terms of its density and hazard rate shapes as shall be seen later. A practical importance, in the context of reliability, for the GaL distribution is shown by modeling two failure data sets ( see application section) with comparing it to many recent versions of the Lindley distribution. Further, in several practical situations, objects in a certain population differ substantially from each other, hence this heterogeneity must be taken into account for accurate analysis on data of such population. Therefore, another motivation of the GaL distribution is represented by its mixture representation which is recommended for such heterogeneity. In the context of Bayesian inference, the GaL distribution given by (1) arises when the gamma f (x | λ) represents the distribution of future observations and the Lindley f (λ | β) is the posterior distribution of the parameters of f (x | λ), given the information in a sample of observed data. By making use of (3), we find that lim x→∞ e δx R(x) = ∞, for all δ > 0, as e δx is unbounded for any δ and grows fastly to ∞ as x → ∞. Hence, the GaL is a heavy tail distribution . Further, it can be noted that lim x→∞ R(t + x)/R(x) = 1, the GaL distribution has also a long right tail.
The rest of the paper is organized as follows. The shape characteristics of the pdf, hazard rate function (hrf) and reverse hazard rate function (rhrf) of the GaL distribution are discussed in Section 2. Quantiles, moments and extreme order statistics are discussed in Section 3. Also, some stochastic order relations are discussed. Estimation by the methods of maximum likelihood and least squares is presented in Section 4. Estimation of the stress-strength parameter is discussed in Section 5. In Section 6, we suggest two different algorithms for generating random data from the new distribution. Further, a simulation study is given to compare the performances of the proposed estimators under the three estimation methods. Finally, two real-life data applications are illustrated in Section 7 for evaluation of the GaL distribution among some recent models.

Shape characteristics
In this section, we discuss the shape characteristics of the pdf, hrf and rhrf of the GaL distribution.

Shape of pdf
We can see from (1) that and lim x→∞ f (x) = 0. Figure 1 shows the pdf of the GaL distribution for some selected choices of α and β. From Figure 1, we see that the pdf of GaL distribution is decreasing for α ≤ 1 and unimodal for α > 1. The shape characteristics of the pdf of GaL distribution are discussed theoretically in the following theorem.
Theorem 1. The pdf of GaL distribution given by (1) is decreasing for α ≤ 1 and unimodal for α > 1.
Proof. The logarithm of (1) is where C is a real valued constant. We have If α ≤ 1, we easily see that So in this case, f (x) is decreasing for all x.

Shapes of hazard rate and reversed hazard rate functions
The hrf and rhrf corresponding to (1) and (2) are given, respectively, by and The behavior of h(x) when x → 0 and x → ∞ , respectively, are given by Also, we can find that lim x→0 r(x) = ∞ and lim x→∞ r(x) = 0. Figures 2 and 3, show the hrf h(x) and rhrf r(x) of the GaL distribution for some choices of α and β, respectively. For establishing the main result regarding the shape of the hazard function h(x) in (4), we use the following lemma due to (Glaser 1980).
Lemma 1. Suppose f (t), for t > 0, is the density function of a positive real-valued continuous random variable, f (t) is the derivative of f (t), and η(t) = −f (t)/f (t). Then, if there exists a t 0 such that η (t) > 0 ∀ t ∈ (0, t 0 ), η (t 0 ) = 0 and η (t) < 0 ∀ t ∈ (t 0 , ∞), the hazard function corresponding to f (t) is either an upside down or a decreasing function of t.
Proof. See (Glaser 1980). The following theorem shows that there are two shapes for the hazard rate function of the GaL distribution, depending on the values of the parameter α and one shape for the reversed hazard rate function of the GaL distribution.
Theorem 2. a) The hazard rate function of the GaL distribution in (4) is decreasing for α ≤ 1 and unimodal for α > 1.
b) The reversed hazard rate function of the GaL distribution in (5) is decreasing.
Proof. a) For the pdf (1), we have It follows that for all α and β. Therefore, the reversed hazard rate function (5) is decreasing.

Quantiles and moments
The pth quantile x p of the GaL distribution defined by F (x p ) = p, is the root of the equation Note that x p can be used to generate GaL random variates. Moreover, the median of the GaL distribution is obtained by using the equation above for p = 1 2 . Now, we discuss moments of the GaL distribution. If X ∼ GaL(α, β), then we find and generally which implies all moments of the distribution are infinite. Thus, the GaL distribution has no mean and this means in practice that if we take a sample X 1 , X 2 , ..., X n from the GaL distribution, then the average X does not tend to a particular value.
On the other side, the negative moments are useful in several applications, such as life testing problems and estimation purposes. Therefore, we discuss the negative moments for this distribution.
The r th negative moment (about the origin) of the GaL distribution is given by in particular, for r = 1 For α ≤ 1, the negative moments undefined.

Extreme order statistics and their limiting distributions
Let X 1:n , · · · , X n:n be the order statistics of a random sample of size n from the GaL(α, β) distribution with distribution function F (x) given by (2). The cdf of the minimum order statistic X 1:n is given by The cdf of the maximum order statistic X n:n is given by The minimum (maximum) order statistics represents the lifetime of a series (parallel) system in reliability studies.
In the following theorem, we provide the limiting distributions of X 1:n and X n:n for the GaL(α, β) model. Theorem 3. Let X 1:n , · · · , X n:n be the order statistics of a random sample of size n from the GaL(α, β) distribution, then

Proof. (i)
In the GaL(α, β) model, we have, by using L'Hospital rule, Therefore, by Theorem 8.3.6 of (Arnold, Balakrishnan, and Nagaraja 1992), the minimal domain of attraction of the GaL(α, β) distribution is the standard Weibull distribution, proving Part(i). Also, we obtain that Therefore, by Theorem 8.3.2 of (Arnold et al. 1992), the maximal domain of attraction of the GaL(α, β) distribution is the standard inverted exponential distribution, which proves Part(ii). Also, we obtain that a n = 0, b n = F −1 (1 − n −1 ).

Stochastic orders
Stochastic ordering of positive continuous random variables is an important tool to judge the comparative behavior of such variables. For this purpose, we shall recall some basic definitions.
A random variable X 1 is said to be smaller than a random variable X 2 in the It is well known that the likelihood ratio order implies hazard rate order which in turn implies stochastic order, see (Shaked and Shanthikumar 1994) for additional details.

Estimation
In this section, we consider estimation of the unknown parameters of the GaL(α, β) distribution by maximum likelihood, least squares and weighted least squares methods.

Maximum likelihood estimation
Let x 1 , x 2 , · · · , x n be the observed values of a random sample X 1 , X 2 , · · · , X n from the GaL(α, β) distribution. Then, the log-likelihood function is given by .
The log-likelihood function of the two parameters is ln L(α, β) = n ln α + 2n ln β − n ln(1 + β) + n i=1 ln(1 + α + β + x i ) It follows that the maximum likelihood estimators (MLEs) of α and β, say α and β, are the simultaneous solutions of the following equations: For interval estimation of (α, β) and tests of hypothesis, one requires the Fisher information matrix: The elements of this matrix for the log-likelihood function (7) are It can be shown that The expressions above follows by using the series see the website "https://www.wolframalpha.com".

Least squares and weighted least squares estimators
In this section, we present the regression based method estimators of the unknown parameters of the GaL distribution. The method of ordinary least squares and the method of weighted least squares were originally proposed by (Swain, Venkatraman, and Wilson 1988) to estimate the parameters of Beta distributions.
Suppose X 1 , ..., X n is a random sample of size n from a distribution function F (.) and X 1:n < ... < X n:n be the ordered statistics of the sample. The least squares estimators (LSEs) can be obtained by minimizing with respect to the unknown parameters of F (.). Therefore in the case of GaL distribution, the least squares estimators of α and β, sayα LSE ,β LSE , respectively, can be obtained by minimizing with respect to α and β.
While, the weighted least squares estimators (WLSEs) of the unknown parameters can be obtained by minimizing n j=1 w j [F (X j:n ) − j n + 1 ] 2 , with respect to the unknown parameters, where w j = 1 V ar[F (X j:n )] = (n + 1) 2 (n + 2) j(n − j + 1) .
Therefore, in case of GaL distribution, the weighted least squares of α and β, sayα W LSE and β W LSE , respectively, can be obtained by minimizing n j=1 w j 1 1 + β X α i:n [X i:n (1 + β) + (1 + α + β) β] (β + X i:n ) 1+α − i n + 1 2 with respect to α and β. Numerical solutions of those estimators will be obtained by a simulation study in Section 6.

Estimation of the stress-strength parameter R = P (X > Y )
In reliability, the stress-strength model describes the life of a component which has a random strength X subjected to a random stress Y . The component fails at the instant that the stress applied to it exceeds the strength, and the component will function satisfactorily whenever X > Y . In this section, we consider the problem of estimating R = P (X > Y ), under the assumption that X ∼ GaL(α 1 , β 1 ), Y ∼ GaL(α 2 , β 2 ), and X and Y are independently distributed. Then it can be easily seen that y α 1 +α 2 −1 [y(1 + β 1 ) + (1 + α 1 + β 1 )β 1 ][1 + α 2 + β 2 + y] (β 1 + y) α 1 +1 (β 2 + y) α 2 +2 dy.
To compute the maximum likelihood estimator (MLE) of R, let us first obtain the MLEs of α 1 , α 2 , β 1 and β 2 . Suppose x 1 , x 2 , · · · , x n be the observed values of a random sample of size n from GaL(α 1 , β 1 ) and y 1 , y 2 , · · · , y m be the observed values of a random sample of size m from GaL(α 2 , β 2 ). Therefore, the log-likelihood function of α 1 , α 2 , β 1 and β 2 is given by ln L(α 1 , α 2 , β 1 , β 2 ) = n ln α 1 + 2n ln β 1 − n ln(1 + β 1 ) + It follows that the MLEs of α 1 , α 2 , β 1 and β 2 , say α 1 , α 2 , β 1 and β 2 , are the simultaneous solutions of the following equations: Once we obtain α 1 , α 2 , β 1 and β 2 , then, we compute the MLE of R as Here the maximum likelihood approach does not give an explicit estimator for the MLEs of the parameters and hence the MLE of R. In practice, one has to use numerical methods to find the MLEs, such methods are well implemented in MATLAB and R packages .

Generation algorithms and Monte Carlo simulation study
In this section, we propose two different algorithms for generating the random data x 1 , x 2 , · · · , x n from the GaL distribution. Further, a simulation study is given to compare the performances of different estimators using the different estimation methods.

Algorithms
Two proposed algorithms for generating the random data x 1 , x 2 , · · · , x n from the GaL distribution are as follows.
(i) The first algorithm is based on generating random data from the Lindley distribution and conditional gamma distribution.
(ii) The second algorithm is based on generating random data from the inverse CDF of the GaL distribution.

Monte Carlo simulation study
Here, a simulation study is given to compare the performances of maximum likelihood, least squares and weighted least squares estimators of the unknown parameters α and β via Monte Carlo simulation. For a given n and (α, β), we have generated the sample x 1 , x 2 , · · · , x n from the GaL(α, β) model and then obtain the estimates using the preceding estimation methods. We used Algorithm 1 to generate data from the GaL distribution. The simulation experiment was repeated N = 10, 000 times each with sample sizes n = 50, 100, 200 and (α, β) = (0.5, 0.5), (0.5, 1.5), (1.0, 0.5), (1.0, 1.5), (2.0, 0.5), (2.0, 1.5). Note that the selected values of (α, β) give decreasing and unimodal shapes of density as displayed in Fig. 1. Two quantities were examined in this Monte Carlo study: (a) Average bias of the MLE θ of the parameter θ = α, β: Mean square error (MSE) of the MLE θ of the parameter θ = α, β: The results of this study are reported in Table 1. The following conclusions can be noted: i) As expected, the biases and MSEs of the MLEs decrease as n increases. Therefore, the MLEs are asymptotically unbiased and consistent.
ii) Table 1 also shows that for most of cases considered, the biases and MSEs of the LSEs and WLSEs decrease as n increases. iii) With respect to the bias, the LSEs work better than the MLEs and WLSEs. iv) With respect to the MSEs, the MLEs work better than the LSEs and WLSEs for α < 1 but for α ≥ 1, the LSEs and WLSEs work better than the MLEs.
For each model, we estimate the unknown parameters using the maximum likelihood approach. All computations are performed using the statistical software R (see the Appendix A).
Tables 2 and 3, for both data sets and each fitted model, list the MLEs of the parameters and their standard errors (in parentheses) from the above fitted models and the values of the statistics: − log L = − log L( α, β), Kolmogorov-Smirnov (K-S) statistics with their p-values, Akaike information criterion (AIC) and Bayesian information criterion (BIC), which are defined, respectively, by −2 log L + 2q and −2 log L + q log(n), where ( α, β) are the MLEs vector, q is the number of parameters estimated and n is the sample size. The best distribution corresponds to lower − log L, K-S statistic, AIC and BIC values, and large p-values associated with K-S. The values of AIC, BIC, − log L, K-S statistics with their p-values in Tables 2 and 3, indicate that the GaL distribution is a strong competitor to the other distributions commonly used in literature for fitting lifetime data, moreover being the best fitting considering the previous goodness of fit statistics.
In order to assess if the new model is appropriate, Figures 5 and 6 display the histograms of two data sets and the fitted density functions, and plots of the empirical and estimated cumulative distribution functions of these fitted distributions. We can conclude that the GaL distribution is a very suitable model to fit both data sets.

Conclusions
In this paper, we propose a new gamma Lindley distribution by compounding the Lindley and gamma distributions in a different manner than used in the GamL distribution (Zeghdoudi and Nedjar 2015). The new distribution has decreasing and unimodal hazard rates shapes which adapt many applications in reliability and survival analysis. We study some of its mathematical properties such as characteristics of the pdf, hrf, rhrf, quantiles, moments, extreme order statistics and their limiting distributions, and stochastic ordering. The model parameters are estimated by the maximum likelihood, least squares and weighted least squares methods. A simulation study is given to compare the performances of the proposed estimators under the mentioned estimation methods. Further, Estimation of the stress-strength parameter is discussed. The potentiality of the new model among recent versions of the Lindley and gamma distributions is illustrated by means of two failure data sets, and using some goodness-of-fit statistics.