An Extended Lindley Distribution with Application to Lifetime Data with Bayesian Estimation

A four-parameter extended of Lindley distribution with application to lifetime data is introduced. It is called extended Marshal-Olkin generalized Lindley distribution. Some mathematical properties such as moments, skewness, kurtosis and extreme value are derived. These properties with plots of density and hazard functions are shown the high flexibility of the mentioned distribution. The maximum likelihood estimations of proposed distribution parameters with asymptotic properties of these estimations are examined. A simulation study to investigate the performance of maximum likelihood estimations is presented. Moreover, the performance and flexibility of the new distribution are investigated by comparing with several generalizations of Lindley distribution through two real data sets. Finally, Bayesian analysis and efficiency of Gibbs sampling are provided based on the two real data sets.


Introduction
The statistical analysis and modeling of complex real data sets such as lifetime data are incapable to do with well-known standard distributions. This matter is created some grave concerns among researchers who attempt for modeling using new distributions and a lot of attention has been recently paid in researches. Nevertheless, many types researches in their pursuit have attempted for modeling complex real data by introducing new family distributions with three or four parameters. To this end, researchers have recently suggested diverse methods for generating new continuous family distributions in concerning to lifetime data analysis such that these distributions are able to fit various lifetime real data sets that have very complex with a high degree of kurtosis and skewness. Such distributions are usually obtained based on the T-X class which is defined by Alzaatreh et al. [8]. The distribution theory has widely extended using the T-X class. Given the vast amount of papers published after 2013, we can only mention a few of the most recent contributions: Kumaraswamy generated, beta generated, odd log-logistic family.
The most well-known one-parameter statistical distributions that use for modeling lifetime data are as: exponential, gamma, Weibull, Lindley and lognormal. In the recent past, since the cumulative distribution function where ξ = (θ, α, γ), G(x; θ) is the baseline cdf with a parameter vector θ and α, γ > 0 are additional shape parameters. Also, the pdf of the family is where g(x; θ) is the baseline pdf. It is clear that the EE-G is a two-parameter family without considering the parameters of the baseline distribution. Since the EE-G family is more flexible than some other families, so this family can be proposed as a serious competitor for other extended families of distributions.
In order to obtain more flexibility in modeling observed data in comparison with other generalizations that already introduced, we introduce another class of distributions using the ideas of EE-G and Marshal-Olkin families. In fact, we can add another parameter, which is called the Marshal-Olkin parameter, to EE-G family. This family is called extended Marshal-Olkin generalized-G (EMOG-G). Based on the T-X class and given a continuous baseline cdf G(x; θ) with a parameter vector θ, the cdf of the EMOG-G family is given by where β > 0 is the Marshal-Olkin parameter and ξ = (θ, α, γ, β). It is obviously that in special case, the EMOG-G family reduces to EE-G family when β = 1. For α = γ = 1, it reduces to Marshal-olkin family. If β = 1 and α = γ, then it reduces to Exp-G family. By considering α = β = γ = 1, we obtain the baseline distribution G.
Since the EMOG-G is a three-parameter family without considering the parameters of the baseline distribution, so use of the baseline distribution with only one parameter such as exponential and Lindley is more efficient than that of with two or three parameters. On the other hand, because of the Lindley distribution has hrf with increasing, decreasing, unimodal and bathtub shapes and is usually applied for the lifetime real data, so we consider the Lindley distribution as the baseline distribution in this paper. Suppose that the cdf of Lindley distribution is as follows where λ > 0 is the scale parameter and pdf corresponding to (4) is given by By considering G(x; θ) as (4) in Eq. (3) where θ = λ, we obtain a new distribution that is called extended Marshal-Olkin generalized Lindley (EMOGL) distribution. In the next section, we present the cdf, pdf and some properties of EMOGL distribution. The rest of the paper is organized as follows. In Section 2, the EMOGL distribution is introduced. Some properties of the EMOGL distribution such as the nth moment, the nth incomplete moment, moment generating function and Bonferroni and Lorenz curves are presented in Section 3. In Section 4, some asymptotic properties and extreme values of the EMOGL distribution are investigated. The maximum likelihood estimators for the unknown parameters of the EMOGL distribution are obtained in Section 5. A simulation study is presented in Section 6. In Section 7, the applications using two real data sets are reported. Bayesian inference and a Gibbs sampling procedure for the considered data sets are investigated in Section 8. Finally, some conclusions are stated in Section 9.
Some of the possible shapes of the density function (27) for the selected parameter values are illustrated in Figure 1(a). As seen in Figure 1(a), the density function can take various forms depending on the parameter values. It is evident that the EMOGL distribution is much more flexible than the Lindley distribution, i.e. the additional shape parameters cause high flexibility of the EMOGL distribution. Both unimodal and monotonically decreasing and increasing shapes appear to be possible.

AN EXTENDED LINDLEY DISTRIBUTION WITH APPLICATION
By definition of the survival function, we set S(x; ξ) = 1 − F (x; ξ) and then we get Now, using (27) and (8) the hrf of EMOGL is given by By choosing different values of parameters, it observes that the EMOGL distribution has various shapes of the hrf. Some of these shapes are drawn in Figure 1(b). Figure 1(b) shows that the hrf of the EMOGL distribution can have very flexible shapes, such as increasing, decreasing, bathtub followed by upside-down bathtub, and bathtub shapes for the selected values of the model parameters. This attractive flexibility makes the hrf of the EMOGL distribution useful and suitable for nonmonotone empirical hazard behaviors which are more likely to be encountered or observed in real life situations.

Mixture representations for the pdf and cdf
The cdf and pdf of the EMOGL can be written as mixture representations and such forms of cdf and pdf can be used to derive some mathematical properties, e.g. moments, moments of residual life and incomplete moments. To this purpose, first, let us remind inverse of a power series using the following Remark.
where c 0 = 1 b0 and c k = − 1 To obtain the mixture representation of the cdf of EMOGL, note that for any 0 < u < 1, By similar argument, we have where b 0 = a 0 (α) + β − a 0 (γ) and b k = a k (α) − βa k (γ), for k ≥ 1. Now, using Remark 1, we get The Eq. (10) can be interpreted as a linear combination of generalized Lindley distribution. Using this equation, the mixture representation of pdf is given by

Mathematical properties
The formulae obtained in this paper can be simply used in some mathematical and statistical software such as Mathematica, Maple and R. These software be able to deal with complex analytic expressions. To determine mathematical properties of EMOGL distribution, use of some algebraic expansions can be more efficient than computing these properties directly. In what follows, we derive the nth moment, kth central moment and moment generating function of EMOGL distribution. In addition, we provide the nth incomplete moment, mean deviations, Bonferroni and Lorenz curves and present numerical values of skewness and kurtosis using the first four ordinary moments. First of all, assume that X ∼ EM OGL(λ, α, β, γ). Using (11), we define dx.
By using generalized binomial expansion it can be shown that So, the nth moment of EMOGL distribution is given by

AN EXTENDED LINDLEY DISTRIBUTION WITH APPLICATION
The central moments µ k = E(X − µ) k of EMOGL distribution can be derived from (13) as where µ k = E(X k ), µ = µ 1 = E(X) and k is an integer value. The mean and variance of X can be particularly obtained using Eq.s (13) and (14). Furthermore, these equations are used to derive the skewness as and the kurtosis as It is to highlight that the Eq. (13) can be easily computed numerically using the mathematical or statistical software. For this purpose, one can compute this equation for a large natural number, say N , instead of infinity in the sums. Therefore, several quantities of X such as moments, skewness and kurtosis can be computed numerically using (13). Table 1 shows numerical values of the first four ordinary moments, skewness and kurtosis of the EMOGL distribution for different values of parameters λ, α, β, γ. Also, the skewness and kurtosis plots of the EMOGL distribution for selected values of α, β and γ when λ = 2 are drawn in Figure 2. Moreover, it is easy to verify that the moment generating function for EMOGL distribution is given by In order to obtain the nth incomplete moment of the EMOGL distribution let us define where γ(λ, z) = z 0 t λ−1 e −t dt stands for the incomplete gamma function. Note that the second equality of (15) is obtained by generalized binomial expansion. Hence, using (15) the nth incomplete moment of the EMOGL distribution is derived by In what follows, we provide two measures of deviation, i.e. mean deviation about the mean (δ 1 ) and the mean deviation about the median (δ 2 ). By definition of these measures, it is easy to show that where M denotes the median of X. Therefore, it can be verified that measures δ 1 (X) and δ 2 (X) are given by Stat., Optim. Inf. Comput. Vol. 9, December 2021 Table 1. Moments, skewness, and kurtosis of EMOGL distribution for some parameter values.

Bonferroni and Lorenz curves
Bonferroni and Lorenz curves were first presented by Bonferroni [13] and Lorenz [23] to measure the inequality of the distribution for a random variable, respectively. These curves are widely used in the literature related to some fields such as reliability, economics, insurance, and etc. The Bonferroni and Lorenz indexes are defined as is the quantile function. If X ∼ EM OGL(λ, α, β, γ), then it can be shown that the Bonferroni curve of the EMOGL distribution is as follows In order to use the Lorenz curve as a measure of inequality of income, one should investigate the area between the Lorenz curve and the line L(p) = p that is called the area of concentration and such area is important in economics, reliability, insurance, and medicine.

Asymptotic properties and extreme value
In this section, we present some asymptotic properties and extreme values of EMOGL distribution. For simplicity, assume that F (x), f (x) and h(x) are cdf, pdf, and hrf of EMOGL distribution, respectively. If X 1 , ..., Xn is a random sample from EMOGL distribution, then asymptotic properties and extreme values based on the EMOGL distribution are obtained as follows.

Asymptotic properties
The asymptotic of cdf, pdf and hrf of the EMOGL distribution as x → 0 are, respectively, given by The asymptotic of cdf, pdf and hrf of the EMOGL distribution as x → ∞ are, respectively, as follows These equations show the effect of parameters on the tails of the EMOGL distribution.

Extreme value
.. + Xn)/n denotes the sample mean, then by the usual central limit theorem, √ n(X − E(X))/ Var(X) approaches the standard normal distribution as n → ∞. One may be interested in the asymptotic of the extreme values Mn = max(X 1 , ..., Xn) and mn = min(X 1 , ..., Xn). Let τ (x) = 1 λ , we obtain following equations for the cdf in (6) as and  [20], there must be norming constants an > 0, bn, cn > 0 and dn such that as n → ∞. The form of the norming constants can also be determined. For instance, using Corollary 1.6.3 in [20], one can see that bn = F −1 (1 − 1 n ) and an = λ, where F −1 (·) denotes the inverse function of F (·).

Estimation
Point estimation is the first step of statistical inference on the unknown parameters of the underlying population. In order to find point estimations there are different methods which the most well-known of them is the maximum likelihood estimation (MLE) method. In the present paper, we obtain the MLEs of the unknown parameters of the EMOGL distribution. To this end, let X 1 , X 2 , ..., Xn be a random sample from EM OGL(λ, α, β, γ) with observed values as x 1 , x 2 , ..., xn. According to the sample, it is easy to see that the log-likelihood function of (α, γ, λ, β) is as By differentiating from the log-likelihood function with respect to the parameters α, γ, λ and β and after simple algebra, we have We denote the MLEs of (α, γ, λ, β), by ( α, γ, λ, β), where these estimations are derived by solving simultaneously the equations ∂ n ∂α = 0; ∂ n ∂γ = 0; ∂ n ∂λ = 0; ∂ n ∂β = 0. Since Eq.s (19)- (22) are non-linear with respect to parameters, so we can 858 AN EXTENDED LINDLEY DISTRIBUTION WITH APPLICATION not obtain the explicit expression of ( α, γ, λ, β). However, by using some numerical iterative techniques we can solve these equations and compute the global maxima of the log-likelihood. Moreover, using Eq.s (19)- (22), we can present the Fisher information matrix as follows ∂i∂j for any i, j = α, γ, λ, β. According to importance of the Fisher information matrix for interval estimating, we follow asymptotical properties of the MLEs. Let θ = (α, γ, λ, β) T , then under standard regularity conditions (see [21, pp. 461-463] , where θ and K(θ) stand for the MLE of θ and the expected Fisher information matrix, respectively. The asymptotic behavior will be established provided that K(θ) is substituted by the observed Fisher information matrix which is multiplied by 1/n, that is I(θ)/n, approximated by θ, i.e. I( θ)/n.

Simulation study
The performance of the maximum likelihood method is evaluated for estimating the EMOGL parameters using a Monte Carlo simulation study. The mean squared error (MSE) and the bias of the parameter estimates are calculated. We are generated N = 10, 000 samples of sizes n = 50, 55, ..., 500 from the EMOGL distribution with α = 0.5, β = 2, γ = 1.5, λ = 2.5. We are applied the following Algorithm 1 for generating random data from EMOGL distribution. Algorithm 1.
• Step 3. Solve numerically the non-linear equation and compute values of x i for i = 1, · · · , n. • Step 4. Repeat Steps 1 to 3 for N times.

859
To compute the MSE and bias of the MLEs, let ( α, β, γ, λ) be the MLEs of the new model parameters and (s α , s β , s γ , s λ ) be the standard errors of the MLEs. The estimated biases and MSEs are given by for ξ = α, β, γ, λ. Figure 3 displays the numerical results for the above measures. From Figure 3, the following results are deduced: The estimated biases decrease when the sample size n increases, The estimated MSEs tend to zero as n increases, These results reveal the consistency property of the MLEs.

AN EXTENDED LINDLEY DISTRIBUTION WITH APPLICATION
• Odd Burr-Power Lindley [7], OBu − P L(α, β, γ, λ), with distribution function • Extended generalized Lindley [28], EGL(α, γ, λ), has distribution function as  For both data sets, the MLEs of unknown parameters, Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), Cramer Von Mises and Anderson-Darling statistics (W * and A * ) are computed for the EMOGL and all the above distributions. Moreover, the Kolmogorov-Smirnov (K-S) statistic with its corresponding p-value and the maximized loglikelihood function (l( ξ, x)) are considered. All of the computations were carried out using the software R.

First application
The first data are the exceedances of flood peaks (in m 3 /s) of the Wheaton River near Carcross in Yukon Territory, Canada. The data consist of 72 exceedances for the years 1958-1984 rounded to one decimal place and exist in Table 2. These data were analyzed in [2].   Tables 3 and  4, respectively. One can easily observe that the smallest values of AIC, BIC, A * , W * and −l statistics and the largest p-value belong to the EMOGL distribution. Therefore, the EMOGL distribution has better performance than the other considered competitors based on these criteria. In addition, the likelihood ratio (LR) test is applied to compare the EMOGL distribution with its sub-models. For example, a testing hypothesis of H 0 : β = 1 versus H 1 : β = 1 is equivalent to compare the EMOGL distribution with EGL distribution. To investigate this testing hypothesis, the LR statistic can be computed by the following equation    where α * , γ * , and λ * are the MLEs of α, γ, and λ, obtained under H 0 , respectively. Under H 0 and by considering the regularity conditions, the LR test statistic converges to χ 2 r in distribution, where r is equal to difference between the number of parameters estimated under H 0 and the number of parameters estimated in general, that is, under H 0 : β = 1, we have r = 1. Table 5 gives values of the LR statistics with their corresponding p-values. From Table 5, we can see that the obtained p-values are too small so, we can reject all the null hypotheses. In other words, based on the LR criterion we conclude that fitting of the EMOGL has better performance than the considered sub-models for the first data set. Also, the fitted pdfs, cdfs and P-P plots of the considered models for the sake of visual comparison are provided in Figures 4 and 5, respectively. Figure 4 suggests that the EMOGL fits the skewed data very well. Figures 5 shows that the plotted points for the EMOGL distribution best capture the diagonal line in the probability plots. Therefore, the EMOGL distribution can be considered as an appropriate model for fitting the first data set.

Second application
The second data set represents the failure times of 50 components, [26, p. 220]. The data set is reported in Table 6. The MLEs of the unknown parameters and the goodness-of-fit test statistics for the second data set are presented in Tables 7  and 8, respectively. It is easy to see that the EMOGL distribution has smaller values of the AIC, BIC, A * , W * and −l statistics than that of the other considered distributions. Moreover, the largest p-value belongs to the EMOGL distribution in Table 8. Hence, the EMOGL distribution outperforms the other considered competitors on the basis of these criteria. These conclusions can also be drawn visually in Figures 6 and 7. Furthermore, the values of the LR statistics with their corresponding p-values are reported in Table 9. From Table 9, we observe that the computed p-values are too small in comparison with their counterparts so, we can reject all the null hypotheses. That is, accordingly the LR criterion we deduce that fitting of the EMOGL is more efficient than the considered sub-models for the second data set.

Bayesian estimation
Bayesian inference procedure has been taken into consideration by many statistical researchers, especially researchers in the field of survival analysis and reliability engineering. In this section, a complete sample data is analyzed through Bayesian point of view. We assume that the parameters α, β, γ and λ of EM OGL distribution have independent prior distributions as where a,b,c,d,e, f , g and h are positive. Hence, the joint prior density function is formulated as follow:   In the Bayesian estimation, according to that we do not know the actual value of the parameter, we may be adversely affected by loss when we choose an estimator. This loss can be measured by a function of the parameter and corresponding estimator.
Hence, we get joint posterior density of parameters α, β, γ and λ for complete sample data by combining the likelihood function and joint prior density (24). Therefore, the joint posterior density function is given by π * (α, β, γ, λ|x) = Kϕ(α, β, γ, λ)L(x, ξ) and K is given as It is clear from the equation (26) that there is no closed form for the Bayesian estimators under the Five loss functions described in Table 10, so we suggest using an M CM C procedure based on 10000 replicates to compute Bayesian estimators. The corresponding Bayesian point and interval estimation and posterior risk are provided in Tables 11 and 12 for flood peaks data set. Table 12 provides 95% credible and HP D intervals for each parameter of the EM OGL distribution. The posterior samples extracted by using Gibbs sampling technique. Moreover, we provide the posterior summary plots in Figures 8 and  9. These plots confirm that the sampling process is of prime quality and convergence is occurred.

Conclusion
In present paper, a new distribution which is called extended Marshal-Olkin generalized Lindley (EMOGL) was introduced. This distribution was included some important distributions as special cases such as extended generalized Lindley, Marshall-    Figure 11. Plots of Bayesian analysis and performance of Gibbs sampling for the Failure times data set. Autocorrelation plots of each parameter of EM OGL distribution. distribution were provided. These properties were shown the EMOGL distribution has considerable flexibility. Moreover, the MLEs of EMOGL parameters and asymptotic distribution based on the MLEs were presented using Fisher information matrix. Furthermore, a simulation study was carried out to evaluate the efficiency of obtained MLEs and it was illustrated the good performance of these estimations. In addition, the EMOGL distribution was compared with some other extended distributions of Lindley such as generalized Lindley, Beta Lindley, exponentiated power Lindley and odd log-logistic power Lindley using two real data sets. Using these real data sets, it was shown numerically that the EMOGL has better performance and high flexibility in comparing to other extended distributions of Lindley which were considered. Finally, the Bayesian estimation and Gibbs sampling procedure for the considered data sets were discussed.