Morgenstern type bivariate Lindley distribution

In this paper, a bivariate Lindley distribution using Morgenstern approach is proposed which can be used for modeling bivariate life time data. Some characteristics of the distribution like moment generating function, joint moments, Pearson correlation coefficient, survival function, hazard rate function, mean residual life function, vitality function and stress-strength parameter R = Pr(Y < X), are derived. The conditions under which the proposed distribution is an increasing (decreasing) failure rate distribution and positive (negative) quadrant dependent is discussed. Also, the method of estimating model parameters and stress-strength parameter by maximum likelihood is elucidated. Numerical illustration using simulated data is carried out to access the estimates in terms of mean squared error and relative absolute bias.


Introduction
In statistical literature, normal distribution and its associated forms have been studied extensively than any other distribution.One reason for this is the flexibility of the distribution to mathematical treatments.This makes it a preferred choice to model random phenomena.However, when the underlying process generates skewed data or the happening of the event of interest is rare, one has to necessarily use non-normal distributions.Some examples to this kind of situation include modeling medical and economic data.Thus, the construction and study of skewed distributions (univariate and multivariate) is an active field of research in statistics.A detailed review on construction of multivariate distributions using different approaches can be found in [23].Methods of constructing bivariate distributions under discrete and continuous set up are available in [16] and [17] respectively.One simple method of constructing bivariate family of distributions using marginals was proposed by Morgenstern [20].A primary advantage of this method is that the resulting form of the distribution function is less complex and is amenable to mathematical treatments.Also, this method can be used when information about marginals and their dependence structure is available.A generalization of Morgenstern method was proposed by [9], which is known as Farlie-Gumbel-Morgenstern (FGM) family of distributions.There are lots of works available in literature on Morgenstern type distributions.[7] introduced a Morgenstern type bivariate gamma distribution and studied its moments and correlation structure.[12] derived the distribution of product and quotient of variates from Morgenstern type bivariate gamma distribution.[24] derived distribution of concomitant of order statistics arising from Morgenstern family.Estimation of parameters in Morgenstern type bivariate logistic, exponential and uniform distributions using ranked set sampling have been carried out by [4], [5] and [26] respectively.
The corresponding p.d.f is The above p.d.f can be expressed as where f 1 (x) = θe −θx and f 2 (x) = θ 2 xe −θx .Thus, Lindley distribution is a mixture of exponential(θ) and gamma(θ, 2) distributions with mixing proportions θ θ+1 and 1 θ+1 respectively.From (5), it can be observed that Lindley distribution approaches exponential distribution for large values of θ.The mode of the distribution is attained at (1−θ)   θ for θ < 1 and zero for θ ≥ 1.A plot of the p.d.f at θ = 0.3 and 1.3 is given in Figure 1 and Figure 2 respectively.The hazard rate function of Lindley distribution is given by Since this is an increasing function in both x and θ, Lindley distribution is an increasing failure rate (IFR) distribution.As IFR⇒IFRA⇒NBU⇒NBUE, where IFRA, NBU and NBUE denote respectively increasing failure Properties and inferential aspects of Lindley distribution have been studied by [10].[27] introduced a three parameter generalization of Lindley distribution that include exponential and gamma distributions as special cases.[21] developed another generalization of Lindley distribution that has monotone, constant and bathtub shape hazard rate functions.A weighted two parameter Lindley distribution having increasing and bathtub shape hazard rate function is proposed by [11].[2] introduced an extended version of Lindley distribution to make it more flexible in terms of shape of hazard rate function.[8] introduced beta generalized power Lindley distribution and studied its properties.More recent generalization of Lindley distribution is given by [1].Though different forms of univariate Lindley distributions are available, not much work has been attempted under bivariate setup.A bivariate extension of generalized Lindley distribution is proposed by [27] by considering two vectors, (V 1 , V 2 ) and (W 1 , W 2 ) of independent random variables distributed according to gamma (α, θ) and gamma(α + 1, θ) respectively.Apart from this, no other bivariate extension of Lindley distribution is available in the literature.This stands as a motivation to propose an alternate yet simple method of obtaining bivariate Lindley distribution using Morgenstern approach and study some of its properties.The proposed Morgenstern Type Bivariate Lindley Distribution (MTBLD) can be used to model life time of coherent system with dependent components.Another application is that it can be used in analyzing competing risk data arising in clinical trails and epidemiological studies.Also, the joint distribution of two adjacent intervals in a Markov dependent point process can be modelled using MTBLD.The paper is organized as follows.Section 2 gives the definition of MTBLD.The corresponding moment generating function (m.g.f), joint moments and correlation coefficient are derived in Section 3. Section 4 discusses positive (negative) quadrant dependence property of MTBLD.In Section 5, the expression for stress-strength parameter of the proposed bivariate distribution is derived.Section 6 deals with obtaining different reliability measures for MTBLD.In Section 7, estimation of the parameters in MTBLD by maximum likelihood (ML) method is explained.Numerical illustration using simulated data is given in Section 8. Concluding remarks are given in Section 9.

Definition
Let X and Y be two random variables each having Lindley distribution with respective parameters θ 1 and θ 2 .Let F X , G Y denote the corresponding c.d.f's and f X , g Y be the corresponding p.d.f's.Using (1), (2), (3) and (4), the c.d.f and p.d.f of MTBLD are obtained as and A plot of the density function for different choices of parameters is given in Figure 3 and Figure 4.The conditional density of Y given X = x for Morgenstern family is defined as Using ( 9), the conditional density of Y given X = x in MTBLD is obtained as ] .(10) In a similar manner the conditional density of X given Y = y can also be obtained.Survival function of Morgenstern family is of the form Using ( 3) and ( 11), survival function of MTBLD is found to be

Moment Generating Function
In this section, m.g.f, joint moments and correlation coefficient of MTBLD are derived.The m.g.f of Morgenstern family of distributions is defined as where M X (t 1 ) and M Y (t 2 ) are marginal m.g.f's.For Lindley distribution, m.g.f is given by Using ( 13) and ( 14), m.g.f of MTBLD is obtained as where From m.g.f, one can obtain the joint moments as Using ( 16), (r, s) th moment of MTBLD is given by Stat., Optim.Inf.Comput.Vol. 4, June 2016 V.S. VAIDYANATHAN AND SHARON VARGHESE.A 137 From (17), the covariance and correlation (ρ) between X and Y is obtained as and Note that ρ = 0 ⇒ α=0 ⇒ X and Y are independent.Based on the conditional density of Y given X = x in (10), the conditional expectation of Y given X = x in MTBLD is obtained as It can be seen that the above conditional expectation is non-linear in X.In a similar manner it can be shown that conditional expectation of X given Y = y is also non-linear in Y .

Positive Quadrant Dependence
In this section, positive quadrant dependence property of MTBLD is discussed.Positive quadrant dependence is a form of dependence between random variables introduced by [18].Two random variables X and Y are said to be positive quadrant dependent (PQD) if If inequality in ( 21) is reversed, then X and Y are said to be negative quadrant dependent (NQD).In the following theorem, it is established that MTBLD is PQD(NQD) for positive(negative) values of α.
Theorem 1 MTBLD is PQD (NQD) for positive (negative) value of α. Proof: , which is always non-negative for all values of x and y, since c.d.f and survival function takes values from zero to one.Therefore, for positive values of α, αβ(x, y) ≥ 0 ∀x, y.This implies the condition given in (21).Hence, MTBLD is PQD for positive values of α.Similarly, for negative values of α, αβ(x, y) ≤ 0 ∀x, y.Therefore, inequality in ( 21) is reversed, hence MTBLD is NQD for negative values of α.
Thus MTBLD possesses both positive and negative quadrant dependence, while the bivariate Lindley model proposed by [27] is only PQD.

Stress-Strength Parameter
In this section, stress-strength parameter R = P r(Y < X) of MTBLD is derived.Stress-strength parameter plays an important role in studies involving product reliability.In this context, R is considered as a measure of reliability of the system and it gives the probability of strength (X) exceeding stress (Y ).Assuming that strength (X) and stress (Y ) are jointly distributed according to MTBLD with dependence parameter α, R is obtained as ) ) ) Estimate of R can be obtained by substituting estimates of θ 1 and θ 2 in (22) for some specified value of dependence parameter α.

Reliability Measures
In this section, reliability measures like hazard rate, mean residual life and vitality function in the context of MTLBD are derived.

Hazard Rate Function
In statistical literature, bivariate hazard rate function is defined in different ways.One due to Basu [3] is given by Using the above definition, the hazard rate function of MTBLD is obtained as A primary limitation of Basu's definition is that it is defined from ℜ 2 → ℜ i.e. h(x, y) is not a vector quantity.To overcome this limitation, [14] defined bivariate hazard rate function in vector form as follows: where S(.) denotes the bivariate survival function.From ( 12), we get Substituting the above expressions in (25) gives the vector hazrad rate function of MTBLD.This function is an increasing (decreasing) function for positive (negative) values of α as proved in the following theorem.

Theorem 2 MTBLD is IHR (DHR) for positive (negative) values of α.
Proof: To prove MTBLD is IHR for positive values of α, it is sufficient to show that ( 26) and ( 27) are increasing functions in x and y respectively.Consider where h X is univariate hazard rate function which is of the form given in (6).Now for This implies that the term is a positive increasing function in x since F X is an increasing function in x.Also, h X is a positive increasing function in x.Thus − ∂lnS(x,y) ∂x is an increasing function in x.Similarly it can be shown that − ∂lnS(x,y) ∂y is an increasing function in y.Thus for α positive, MTBLD is IHR.In a similar manner, it can be proved that MTBLD is DHR for negative values of α.

Mean Residual Life
Mean residual life (m.r.l) denotes the average remaining life of a unit after it has survived for a specified time t.For vector valued random variables, [25] defined m.r.l as where and The expressions for m 1 (x, y) and m 2 (x, y) in MTBLD is obtained as ) , (29) ) .

Vitality Function
For a system with two components, [22] defined bivariate vitality function as where Also, v i (x, y) is related to m i (x, y) by Here v 1 (x, y) measures the expected life time of first component as the sum of present age x and the average life time remaining to it, given that second component has survived beyond age y. v 2 (x, y) can be interpreted in similar way.Using Equations ( 29), ( 30) in (32) we obtain v 1 (x, y) and v 2 (x, y) of MTBLD as ) , (33) ) .( From ( 33) and (34), vitality function of MTBLD can be obtained using (31).
where C denotes constant term independent of θ 1 and θ 2 .Differentiating lnL partially with respect to the parameters θ 1 , θ 2 , α and equating them to zero, we get the following log likelihood equations. ) ) It can be seen that the above equations are non-linear with respect to the parameters and hence obtaining closed form expressions for the estimators is not possible.However, when the dependence parameter α is fixed at some specified constant, using (36) and (37), the ML estimates (MLEs) of θ 1 and θ 2 can be obtained by using the method given in [13] as fixed point solutions of the above equations.The method of obtaining estimates is explained below.Define where ] −1 . (41) is the MLE of (θ 1 , θ 2 ), then ( θ 1 , θ 2 ) will be a fixed point solution of (39).The MLE of (θ 1 , θ 2 ) can be found by implementing the following iterative procedure.
• Using initial value, obtain (θ 1j+1 , θ 2j+1 ) as a solution of The above process is continued till the difference between successive values of (θ 1 , θ 2 ) is less than some specified threshold limit.The solutions arrived at the final iteration are taken as estimates for the unknown parameters.One may also use two-dimensional Newton-Raphson method or any root finding algorithm to obtain solution to the system of non-linear equations given in (36)-(38).
Taking negative expectations of the second order partial derivatives and mixed derivative with respect to θ 1 and θ 2 , the Fisher information matrix F is given by where From the log-likelihood function given in (35), second order partial derivatives are obtained as Since the second order partial derivatives have complex expressions, finding expected values of the same is difficult.Hence, the sample (observed) Fisher information matrix can be used instead which is given by where

Simulation Study
In order to access the performance of the proposed estimation procedure, a simulation study is carried out by generating random samples from MTBLD.The results are evaluated in terms of mean squared error (MSE) and relative absolute bias (RAB).The following procedure is adopted to generate samples (X , Y i ), i = 1, 2, ..., n from MTBLD with parameters θ 1 ,θ 2 and α.
1. Generate n bivariate observations (U i , V i ), i = 1, 2, ..., n from Farlie-Gumbel-Morgenstern (FGM) copula for given value of α .2. For the specified values of θ 1 and θ 2 , find X i and Y i such that where W (.) is the Lambert's W function.For more details about Lambert's W function refer [6].
In the simulation study, the following choices of parameter values namely, θ 1 = 0.4, 1, 1.9, θ 2 = 0.7, 1, 1.4 and α = 0.9, −0.4 are considered.The number of Monte Carlo (MC) runs is taken to be 1000.For each MC run, random samples of sizes n = 50 and n = 100 are simulated from MTBLD with the different parameter values.Packages 'copula' and 'gsl' available in R are used for generating bivariate samples from FGM copula and evaluating Lambert's W function, respectively.The MLEs of parameters are determined using fixed point solution method discussed in section 7 with initial values θ 1 = 0.2 and θ 2 = 0.4.Also for each MC run, the MLE of stressstrength parameter R is calculated using the MLE of θ 1 and θ 2 in (22).Let Θ k = ( θ 1k , θ 2k , R k ) be the MLE of Θ, Θ = (θ 1 , θ 2 , R) based on k th MC run, k = 1, 2, ..., 1000.Then average MLEs and respective MSEs and RABs are computed as where r denotes the number of MC runs.The results are presented in Table 1 and Table 2.
From Table 1 and Table 2, it is observed that as sample size increases the MSEs and RABs of the estimates decreases and tend towards zero.Also, the average MLEs are closer to their respective true parameter values.Thus, the proposed method produces estimates that are consistent.Since the expected strength E(X) of the system is inversely related to θ 1 through E(X) = θ1+2 θ1(θ1+1) , as θ 1 increases, E(X) decreases, resulting in a decrease in the system reliability.Similarly as θ 2 increases, expected stress E(Y ) decreases, resulting in an increase in R. The same pattern is observed from the estimates of R given in Table 1

Conclusion
A Morgenstern type bivariate Lindley distribution is proposed in this paper and joint moments, correlation and certain reliability characteristics of the same is obtained.It is shown that the proposed model satisfies positive (negative) quadrant dependence property and is IFR (DFR) distribution for positive (negative) values of dependence parameter.An explicit expression for the stress-strength parameter is also derived for the proposed distribution.A method of obtaining MLE of the parameters using fixed point solution is also proposed.Numerical illustration through simulation study reveal that the proposed method results in estimates that are consistent.