The McDonald Lindley-Poisson Distribution

We propose the McDonald Lindley-Poisson distribution and derive some of its mathematical properties including explicit expressions for moments, generating and quantile functions, mean deviations, order statistics and their moments. Its model parameters are estimated by maximum likelihood. A simulation study investigates the performance of the estimates. The new distribution represents a more flexible model for lifetime data analysis than other existing models as proved empirically by means of two real data sets.


Introduction
The Poisson distribution has been used to generate several flexible continuous distributions by compounding methods for modeling survival data. Many generalizations based on the Poisson distribution are investigated in recent years. For example, the Conway-Maxwell-Poisson discussed in Minka et al. (1) and Shmueli et al. (2), exponential Poisson proposed by Kus (3), Weibull Poisson studied by Hemmati et al. (4) and exponentiated Burr XII Poisson proposed and studied by Silva et al. (5). Further, the Lindley distribution was introduced by Lindley (6) to illustrate the difference between fiducial and posterior distributions. Ghitany et al. (7) investigated the properties of the latter distribution and showed that it is a better model than the exponential distribution. The Lindley has also some advantages since to its hazard rate function (hrf) can exhibit bathtub shapes, and then it becomes more versatile and flexible for compounding with other distributions. Sankaran (8) introduced the Poisson-Lindley distribution by assuming that the Poisson parameter follows a Lindley distribution. Mahmoudi and Zakerzadeh (9) generalized the Poisson-Lindley distribution and showed that their generalization has more flexibility in analyzing count data. Bhati (11). Modeling real life data is very useful in many applied areas and considerable efforts have been done to construct new distributions for survival data. However, there still remain many problems involving real data, which are not contemplated by existing distributions. Researchers in recent years proposed new distributions by extending the exponential, gamma and Weibull, among others. Some of these generalizations have been used extensively for modeling data in several areas such as the actuarial, engineering, biological studies, economics, finance and medical sciences. However, extra works are needed for extended forms of these distributions in many applied areas of lifetime analysis. For an arbitrary parent cumulative distribution function (cdf) G(x; τ ) with parameter vector τ and g(x; τ ) = dG(x;τ ) dx , the probability density function (pdf) of the McDonald class ("Mc-G" for short) proposed by Alexander et al. (12) is where a > 0, b > 0 and c > 0 are additional shape parameters which introduce skewness and provide greater flexibility of the tails of f (x). The class (1) includes two important family models: the beta generalized family (Eugene et al. (13)) when c = 1 and Kumaraswamy generalized family (Cordeiro and Castro,(14)) when a = 1. The Mc-G distribution with baseline G(x; τ ) is simply the beta generalized distribution with baseline G(x; τ ) c . This simple transformation may facilitate the computation of several of its properties. For G(x; τ ) = x, we obtain the McDonald distribution and the beta and Kumaraswamy distributions follow as special cases when c = 1 and a = 1, respectively. The cdf corresponding to (1) is where is the quantile function (qf) of G and V has a beta distribution with parameters a and b. One important characteristic of the Mc-G class is its ability to fit skewed data that cannot be properly fitted by existing distributions. In this paper, we study the McDonald Lindley-Poisson (McLP) distribution, which extends the Kumaraswamy Lindley-Poisson explored and studied by Pararai et al. (11). The motivation to study the McLP distribution is because of the extensive use of Lindley distribution in survival analysis and the fact that existing generalizations can be improved. Some important contributions of this paper are: The proposed distribution overcomes a limitation of its baseline, whose hazard function presents only monotonic increasing and upside-down bathtub shapes. The McLP hrf admits the five main characteristics: monotonically increasing, monotonically decreasing, bathtub, upside-down bathtub and sigmoid shapes.
An advantage of fitting a wider model to real data is that we can verify, using the likelihood ratio (LR) statistics, whether its sub-models with fewer parameters, can be more proper to the data. The McLP distribution extends eight lifetime distributions, including the beta LP (BLP) and exponentiated LP (ELP) models. The BLP distribution becomes a very competitive model to some well-known fourparameter lifetime distributions such as the Burr XII Poisson, beta inverse Weibull, beta exponentiated Lindley and ELP.
Although the proposed model has five parameters, it can provide better fits, based on Anderson-Darling and Cramér-von Mises statistics, than ones with less number of parameters.
The paper is organized as follows. In Section 2, we define the McLP distribution and highlight some special cases. In Section 3, we demonstrate that the McLP density function is a linear combination of LP densities. Some mathematical properties of the McLP distribution are also obtained in this section including ordinary The McDonald Lindley-Poisson Distribution and incomplete moments, mean deviations and generating and quantile functions, order statistics and their moments. The estimation of the model parameters by maximum likelihood is addressed in Section 4. In Section 5, we perform a simulation study to verify the adequacy of the estimates. Two applications to real data illustrate the flexibility of the new distribution in Section 6. Section 7 offers some concluding remarks.

The McLP distribution
Suppose that Z has a truncated Poisson distribution with parameter θ > 0 and probability mass function We consider independent and identically distributed (iid) random variables {W i } Z i=1 having Lindley density function where β > 0 is the scale parameter.
Assuming that W 's and Z are independent random variables, X = min{W 1 , . . . , W Z } defines the LP distribution. The pdf of X takes the form (for x > 0) g(x; β, θ) = The cdf corresponding to (3) is The LP model is well-motivated for biological studies. For example, consider the time to relapse of cancer under the first-activation scheme. Suppose that the number of carcinogenic cells for an individual left active after the initial treatment follows a truncated Poisson distribution and let W i be the time spent for the ith carcinogenic cell to produce a detectable cancer mass (for i ≥ 1). If {W i } i≥1 is a sequence of iid Lindley random variables independent of Z, then the time to relapse of cancer of a susceptible individual can be modeled by the LP distribution. The five-parameter McLP distribution is defined from (1) by taking G(·) and g(·) to be the cdf and pdf of the LP distribution. Its density function is (for where ϕ(x) = θ 1 − (1 + βx β+1 ) e −βx , θ > 0, β > 0 is a scale parameter, a > 0, b > 0 and c > 0 are shape parameters. Henceforth, the McLP distribution is denoted by the random variable X ∼McLP(a, b, c, β, θ). According to Gui et al. (15), the LP distribution has only monotonically increasing or upside-down bathtub hazard rates. On the other hand, as we shall see later, the hrf of the McLP distribution allows for monotonically increasing, monotonically decreasing, bathtub, upside-down bathtub and sigmoid shapes. The density function (4) includes as special models some distributions. In fact, the LP distribution is a basic exemplar when a = b = c = 1. The BLP distribution is obtained when c = 1. For a = 1, we have the Kumaraswamy LP (KwLP) distribution. The McLP distribution becomes the ELP when a = b = 1. In addition, we obtain the exponentiated Lindley distribution when θ → 0 + . The cdf and hrf of X are, respectively, Some possible shapes of the density function (4) and hrf (6) for selected parameter values are displayed in Figures 1 and 2, respectively.
The random variable X with cdf (5) has a simple stochastic representation X = Q G (V 1/c ), where x = Q G (u) = G −1 (u) denotes the qf corresponding to G and V has a beta distribution with parameters a and b. The McLP distribution can be simulated by inverting Equation (5) as follows: If V is a beta random variable with parameter a and b, then follows the McLP(a, b, c, β, θ) distribution. This method is useful because of the existence of strong generators for beta random variables.

Quantile function
The qf of the McLP distribution, say Q(u) = F −1 (u), can be expressed in terms of the beta qf. By inverting (5), we can write where W −1 (·) is the negative branch of the Lambert function (see Corless  can be found the following expansion for the beta qf where the quantities A i (i = 1, 2, 3, 4) are and

Linear representation
In this section, we derive a useful linear representation for the McLP density function. If |z| < 1 and d > 0 is a real non-integer, the power series holds where the binomial coefficient is defined for any real. Based on Equation (8) and by expanding the binomial in (1), we have where the weights are which leads to a linear representation where θ j = (j + 1)θ > 0, g(x; β, θ j ) denotes the LP density and the weights are .
By integrating (10), we have where G(x; β, θ j ) denotes the cdf of the LP distribution.

The McDonald Lindley-Poisson Distribution
Some mathematical properties of the McLP distribution can be determined directly from (10) using the properties of the LP distribution. These properties can be found analytically in softwares such as MATHEMATICA, MATLAB and MAPLE since these computing softwares have the ability to deal with most complex analytical expressions.

Moments
In this section, we obtain the moments of the McLP distribution. Some of the most important features and characteristics of a distribution can be studied through moments (e.g., tendency, dispersion, skewness and kurtosis). The nth ordinary moment of X comes from Equation (10) as By using the power series for the exponential function, we obtain By using the binomial expansion, (12) can be rewritten as The last integral can be computed using the MATHEMATICA software. Then, dt is the well-know hypergeometric function (Gradshteyn and Ryzhik (17)).
An expression for the moment generating function (mgf) of X can be obtained from Equation (10) using the LP density function. We can express it as
Setting u = G(x) in (11) gives where the integral T k (z) can be expressed in terms of Q(u) = G −1 (x) as

The McDonald Lindley-Poisson Distribution
An alternative representation for m 1 (z) can be derived from (10) as The Lorenz and Bonferroni curves are important applications of the mean deviations in fields like economics, reliability, demography insurance and medicine. They are defined for a given probability π by B(π) = m 1 (q)/(πµ 1 ) and L(π) = m 1 (q)/µ 1 respectively, where µ 1 = E(X) and q = Q(π) is given by (7). The Bonferroni and Lorenz curves for the McLP distribution as functions of π are readily calculated from (13) for n = 1.

Order statistics
Order statistics make their appearance in many areas of statistical theory and practice. Order statistics play an important role in quality control and reliability, where some predictors are often based on their moments. We derive an explicit expression for the density function of the ith order statistic X i:n , say f i:n (x). Suppose X 1 , X 2 , . . . , X n independent and identically distributed random variables from the McLP distribution. Let X i:n denote the ith order statistic. From Equations (4) and (5), the pdf of X i:n can be expressed as an infinite linear combination of LP density functions f i:n (x) = ∞ k,r,t=0 n−i j=0 p k,r,t g(x; β, θ(t + 1)), (14) where and d j+i−1,r are defined recursively by Clearly, the cdf of X i:n can be expressed as n−i j=0 p k,r,t G(x; β, θ(t + 1)).
Equation (14) reveals that the pdf of X i:n can be represented as a finite mixture of Lindley Poisson densities. Thus, some mathematical properties of X i:n can be obtained from (14). For example, the moments and mgf of X i:n are given by where Z ∼ LP(β, θ(t + 1)).

Estimation
We calculate the maximum likelihood estimates (MLEs) of the parameters of the McLP distribution from complete samples only. Let x 1 , . . . , x n be a random sample of size n from the McLP(a, b, c, β, θ) distribution. The log-likelihood function for the vector of parameters θ = (a, b, c, β, θ) T can be expressed as The elements of the score vector are given by

The McDonald Lindley-Poisson Distribution
Note that since ϕ(x) = θ 1 − 1 + β x 1+β e −βx , we have The maximum likelihood estimates, θ of θ = (a, b, c, β, θ) are obtained by solving the nonlinear equations ∂ ∂a = 0, ∂ ∂b = 0, ∂ ∂c = 0, ∂ ∂β = 0, ∂ ∂θ = 0. These equations are not in closed form and the values of the parameters a, b, c, β and θ, must be found by using iterative methods. We maximize the likelihood function using NLmixed procedure in SAS as well as the function nlm in R. These functions were applied and executed for a wide range of initial values. This process often results or leads to more than one maximum, however, in these cases, we take the MLEs corresponding to the largest value of the maxima. In a few cases, no maximum was identified for the selected initial values. In these cases, a new initial value was tried in order to obtain a maximum. The issues of existence and uniqueness of the MLEs are of theoretical interest and have been studied by several authors for different distributions (Seregin (18); Santos Silva and Tenreyro (19); Xia et al. (20); Zhou (21)). At this point we are not able to address the theoretical aspects (existence, uniqueness) of the MLE of the parameters of the McLP distribution.

A simulation study
To examine the performance of the MLEs for the McLP distribution, we perform a simulation study: 1. generate r samples of size n from (4) using (7).
2. compute the MLEs for the r samples, say a i , b i , c i , β i , θ i for i = 1, . . . , r.
3. compute the standard errors of the MLEs for r samples, say s a , s b , s c , s β , s θ for i = 1, . . . , r. The standard errors are computed by inverting the observed information matrices.

compute the average biases (bias) and root mean squared errors (RMSE) by
for = a, b, c, β, θ.
We repeat these steps for r = 1, 000 values and n = 50, 100, . . . , 500 with a = 4.0, b = 5.0, c = 1.5, β = 2.3 and θ = 7.0, so computing bias (n) and RMSE (n). Table 1 presents the bias and RMSE values. The biases of the estimates decrease when n → ∞. The reported results hold only for the choice (a, b, c, β, θ) = (4.0, 5.0, 1.5, 2.3, 7.0). However, the results are similar for a wide range of this parameter vector. We also observe that as the sample size n increases, the RMSEs likewise decreases.

Application
In this section, we present two applications of the McLP model to the data obtained from Bader and Priest (22) (data set I) and Bjerkedal (23) (data set II). The data set I consist of 56 strength data measured in GPA, the single carbon fibers and impregnated 1000-carbon fiber tows. Single fibers were tested under tension at gauge length 1 mm. The second real data set represents the survival times (in days) of 72 guinea pigs infected with virulent tubercle bacilli. In Table 2 we provide some descriptive statistics for both sets of data. The total time on test (TTT) plot proposed by Aarset (24) is a means to check the shape of the observed hazard function. This is drawn by plotting T (i/n) = i r=1 y (r) + (n − i) y (r) / i r=1 y (r) , where, i = 1, 2, . . . , n and y (r) (r = 1, 2, . . . , n) is the order statistics of the sample, against i/n. For constant hazard plot is a straight diagonal while for decreasing (increasing) hazards it is convex (concave) respectively. The TTT plots for the data sets ( Figure 3) indicate that all the data sets have increasing hazard rate.  Table 3. We have considered known model selection criteria namely the AIC (Akaike Information Criterion), BIC (Bayesian Information Criterion), CAIC (Consistent Akaike Information Criterion) and HQIC (Hannan-Quinn Information Criterion) and Anderson-Darling (A * ) and Cramér-von Mises (W * ) goodness-of-fit statistics to compare the fitted models ( Table 4). The statistics A * and W * are described by Chen and Balakrishnan (25). In general, the smaller the values of these statistics, the better the fit to the data.  Let θ = (a, b, c, β, θ) be a parametric vector, the LR test statistic is given by LR = −2[ (θ * ; x) − (θ; x)] whereθ * is the restricted MLEs under the null hypothesis H 0 andθ is the unrestricted MLEs estimates under the alternative hypothesis H 1 . Under the null hypothesis, the LR criterion follows Chi-square distribution. The null hypothesis can not be accepted for p-value less than 0.05. A comparison of the McLP distribution with some of its nested models using LR statistics is performed in Table 5.
The Tables 4 and 5

Conclusions
The Lindley-Poisson distribution is commonly used to model the lifetime of a system. However, it does not exhibit a bathtub-shaped failure rate function and thus it can not be used to model the complete lifetime of a system. We define a new lifetime model, called the McDonald Lindley-Poisson distribution, which extends the Kumaraswamy Lindley-Poisson distribution proposed by Pararai et al. (11), whose failure rate function can be increasing, decreasing and bathtub. The McDonald Lindley-Poisson density function can be expressed as a linear combination of Lindley-Poisson densities, which allows to obtain several of its structural properties. We provide a mathematical treatment of the distribution including explicit expressions for the density function, generating function, ordinary and incomplete moments, generating function, quantile function, mean deviations and order statistics and their moments. The parameter estimation is approached by maximum likelihood and the observed information matrix is derived. The usefulness of the new model is illustrated in applications to real data using formal goodness-of-fit tests. By means of two real data sets The McDonald Lindley-Poisson Distribution application, we show that the proposed distribution can give a better fits them other competitive models.