On the Bivariate Generalized Gamma-Lindley Distribution

1. Introduction

Statistical distributions are intrinsic in terms of describing and predicting real world phenomena. They are frequently used in life data analysis and reability engineering. However, many of these classical distributions do not provide enough flexibility for analysing different types of lifetime data and are not sufficient for modelization in several other areas such the discipline of signal processing applications or the area of image processing and computer vision (for more details, see [15]). Therefore, a new class of distribution was introduced by Lindley [12] to analyse failure time data, especially in application modeling stress-strength reability (see, for example, [5]). This new class of distribution can be also used in a wide variety of fields including biology, engineering and medicine. The one-parameter Lindley distribution is defined by its probability density function (pdf) expressed as

It has been explored and generalized by numerous authors. Basically, Ghitany et al. [5] provided a treatment of the mathematical properties of equation (1). Next, Shanker et al. [22] derived a two-parameter Lindley distribution of which the one parameter Lindley distribution is a particular case, and they demonstrated its applicability in modeling waiting and survival times data. In [26], Zakerzadeh and Dolati introduced a three-parameter generalization of the one-parameter Lindley distribution and examined multiple properties of this new distribution. Furthermore, Torabi et al. [24] set forward a four parameters extension of the generalized Lindley distribution. In addition, discrete versions of the Lindley distribution were proposed in literature by Gómez-Déniz and Calderín-Ojeda [7], Gómez-Déniz et al. [8], Hussain et al. [10] and many other authors. Zeghdoudi and Nedjar [27] proposed another extension of Lindley distribution called Gamma Lindley distribution which offers a more flexible model for modeling lifetime data and they investigated its properties. Recently, Laribi et al. [13] derived a new generalized Gamma-Lindley distribution through mixing the Gamma and the three-parameters generalized Lindley distributions. Its pdf is (for x, α, θ > 0, γ ≥ 0 and β ≥ 1) indicated by

where Γ(α):=∫0∞xα−1e−xdx is the Euler Gamma function.

In this paper, we generate a new class of bivariate distributions called Bivariate Generalized Gamma-Lindley (BGGL) distribution. The joint pdf of this new distribution corresponds to

where x₁, x₂, α, θ > 0, γ ≥ 0, β ≥ 1 and Θ := (α, β, γ, θ).

The paper is organized as follows: In Section 2, we display some properties of our distribution such as joint moments, correlation and mixture construction of the BGGL distribution defined by equation (3). Additionally, for different values of parameters, we provide graphically a variety of forms of the BGGL density. In Section 3, we derive the exact densities and moment functions of the random variables Y = X₁ + X₂, Z = X₁/(x₁ + x₂), V = X₁/X₂ and W = X₁X₂, where (X₁, X₂) follows the BGGL distribution. Furthermore, we exhibit approximations for the distribution of W = X₁X₂, and discuss evidence for their robustness. The objective consists in providing simple approximations in terms of the beta and inverted beta distributions so as to facilitate the use of the known procedures for inference, prediction, etc. In Section 4, we determine the exact forms of the Rényi and Shannon entropies for the BGGL distribution. In Section 5, we investigate the estimation of the BGGL distribution parameters by the E.M algorithm. In the last section, we provide an application of the results within the field of precipitation data about rain and snow from France and we tackle another application in the field of demographic studies in USA.

2. Properties

In this section, we handle several outstanding properties of the BGGL distribution. This bivariate distribution is based essentially on a mixture of the product of two independent Gamma random variables (denoted by gamma(α, θ)) and bivariate generalized Lindley distribution (denoted by BGL(α,γ, θ)), introduced by Zakerzadeh and Dolati [26]. Both distributions are defined by the joint pdfs

and

where x₁, x₂, α, θ > 0, γ ≥ 0. Let (Y₁, Y₂) and (Z₁, Z₂) be two random vectors, such that Y₁ and Y₂ are independent random variables distributed according to gamma(α + 1, θ) and (Z₁, Z₂) is a vector distributed according to BGL(α, γ, θ). For β ≥ 1, consider the random pair (X₁, X₂) which stands for the mixture between (Y₁, Y₂) with probability (β−1)β, and (Z₁, Z₂) with probability 1β. This mixture law is called BGGL distribution and the joint pdf of (X₁, X₂) is obtained by using the total probability law:

where x₁, x₂, α, θ > 0, γ ≥ 0 and Θ := (α, β, γ, θ), which coincides with equation (3). Note that, the joint pdf (3) may be written in terms of the gamma function as:

It is easy to demonstrate that the univariate marginal density functions of equation (3) are of the form (2). When β = 1, the BGGL distribution reduces to the bivariate generalized Lindley distributions with three parameters (see [26]). For β = 1 and α = 1, the BGGL distribution reduces to the bivariate generalized Lindley distributions with two parameters (see [23]). For β = 1 and γ = 0, the random variables X₁ and X₂ become independent and the joint pdf (3) reduces to the product of two gamma density functions with the same parameters.

In addition, Figure 1 depicts a few graphs of the pdf defined by equation (3) for different values of parameters. Here, one can record the wide range of forms that result from the BGGL distribution.

Figure 1.

Plots of the pdf of BGGL distribution for some selected values of parameters (α, β, γ, θ) as: (1) (0.05, 1.2, 1.5, 0.3); (2) (0.05, 1.2, 1.5, 2) ; (3) (2, 1.2, 1.5, 0.3); (4) (2, 1.2, 1.5, 2) ; (5) (7, 1.2, 1.5, 0.3); (6) (7, 1.2, 1.5, 3).

Furthermore, simple calculations yield that for each positive integer m and n, the joint (m, n)-th moment of the variables X₁ and X₂ having BGGL distribution, can be stated as

Now, substituting appropriately for m and n in equation (5), we obtain

Besides, we can easily show that the correlation coefficient of X₁ and X₂ is indicated by

It is straightforward to notice that if α → 0, β = 1 and θ → ∞ or if α → 0, γ = 0 and β = 1, we have Corr(X₁, X₂) → 1/2; which corresponds also to the maximum value of the correlation for this family. The BGGL distribution stands for a mixture of i.i.d random variables since the BGL distribution is itself a mixture of independent random variables. Relying upon the result of Drouet-Mari and Kotz [3: Section 3.2.8], we infer that X₁ and X₂ are positively correlated by mixtures.

A bivariate distribution is said to be positively likelihood ratio dependent if the density f_Θ(x₁, x₂) satisfies

for all x₁ > x₂ and y₁ > y₂ (see Lehmann [14]). For the equation (3), the above inequality reduces to

or x₂y₂ + x₁y₁ ≥ x₁y₂ + x₂y₁, which clearly holds. Hence, the BGGL distribution is positively likelihood ratio dependent.

3. Some transformations

In this section, we derive explicit expressions for the pdfs and moments of Y = X₁ + X₂, Z = X₁/(X₁ + X₂), V = X₁/X₂ and W = X₁X₂ when the random vector (X₁, X₂) is set according to BGGL distribution. In fact, for given random variables X₁ and X₂, the distributions of the product, sum, and ratio arise in many fields as biology, economics, engineering, genetics, medicine, number theory, physics, etc. (see for example [17]). For instance, if X₁ and X₂ denote the rainfall intensity and the duration of a storm, then W = X₁X₂ will represent the amount of rainfall produced by that storm. Another expressive example is the following: if X₁ and X₂ denote the proportion of days with rain and the proportion of days with snow, then Y = X₁ + X₂ will represent the proportion of days with precipitation (rain or snow). Moreover, here is an example of the quotient V = X₁/X₂ which consists in representing X₁ as the crude divorce rate and X₂ as the crude marriage rate. Then, V will represent the divorce to marriage ratio. Furthermore, we set forward approximations for the distribution of W = X₁X₂, and discuss evidence for their robustness.

3.1. PDFS

Theorems 1–4 and Corollary 2 yield the pdfs of Y, Z, V and W when X₁ and X₂ are distributed according to equation (3).

Theorem 1.

If X₁ and X₂ are jointly distributed according to equation (3), then the pdf of Y is expressed by

6 fY(y) = θ2α+1y2α−1e−θyβ(γ + θ) Γ (2α)(1 + (βθ + βγ − θ)θy22α (2α + 1)),     y > 0,

and the pdf of Z is given by

7 fZ(z) = zα−1(1 − z)α−1Γ (2α)β(γ + θ)Γ2 (α + 1) (2α (2α + 1) (βθ + βγ − θ)z(1 − z) + α2θ),     z ∈ (0, 1).

Proof. Substituting x₁ = yz and x₂ = y(1 – z) with the Jacobian J(x₁, x₂ → y, z) = y, in the joint pdf of X₁ and X₂, we obtain the joint pdf of Y and Z as

8 fY,Z(y, z) = θ2α+1zα−1(1 − z)α−1 y2α−1e−θyβ(γ + θ) Γ2(α + 1)[ (βθ + βγ − θ) θy2z (1 − z) + α2 ],

where z ∈ (0,1) and y > 0. Now, integrating this equation appropriately, we obtain marginal densities of Y and Z.

Remark 1.

Simple calculations reveal that equation (8) can be rewritten as

fY,Z(y, z) = θβ(γ + θ) fg(y; 2α, θ) fb(z; α, α) + (βθ + βγ − θ)β (γ + θ) fg(y; 2α + 2, θ) fb(z; α + 1, α + 1),

where f_g (y; a, b) and f_b (z; a, b) denote the pdfs of the gamma and the beta distributions with the parameters a and b, respectively, for all 0 < z < 1 and y > 0.

Corollary 2.

If X₁ and X₂ are jointly distributed according to equation (3), then the pdf of the variable V = X₁/X₂ is expressed by

9 fV(v) = vα−1Γ(2α)β(γ + θ)(1 + v)2αΓ2(α + 1)(2α (2α + 1) (βθ + βγ − θ)v(1 + v)2 + α2θ),v > 0.

Proof. We can state that

X1X1 + X2 = X1/X2X1/X2 + 1

which, in terms of Z = X₁/(X₁ + X₂) and V = X₁/X₂, can be written as Z = V/(1 + V). Therefore, making the transformation Z = V/(1 + V) with the Jacobian J(z → v) = (1 + v)^-2 in the pdf of Z given in (7), we get the pdf of V.

In the following, we derive the exact pdf of the product W = X₁X₂. The calculations involve the modified Bessel function of the second kind defined by:

Kv(x) = π (I−v(x) − Iv(x))2 sin(vπ),

where I_ν(·) denotes the modified Bessel function of the first kind of order ν defined by

Iv(x) = xv2vΓ (v + 1) ∑k=0∞1(v + 1)k k!(x24)k.

Furthermore, K₀(·) is interpreted as the limit K0(x)=limν→0Kν(x). We also need the following important lemma:

Lemma 3

([18: Equation 2.3.16.1]). For p > 0 and q > 0,

∫0∞xα−1e−px−q/xdx = 2(qp)α/2 Kα (2pq) .

Theorem 4.

If X₁ and X₂ are jointly distributed according to equation (3), then the pdf of the variable W = X₁X₂ is provided by

10 fW(w) = 2θ2α+1wα−1 ((βθ + βγ − θ) θw + α2)β(γ + θ)Γ2(α + 1)K​0 (2θw),w > 0.

Proof. From (3), the joint pdf of (X₁, W) = (X₁, X₁X₂) with the Jacobian J(x1,  x2→x1,  w)=x1−1, becomes

fX1, W (x1, w) = θ2α+1wα−1 ((βθ + βγ − θ) θw + α2)β(γ + θ)Γ2(α + 1)x1−1e−θ(x1 + wx1).

Hence, integrating this equation with respect to x₁ we find that

11 fW(w)=θ2α+1wα−1((βθ+βγ−θ)θw+α2)β(γ+θ)Γ2(α+1)∫0∞x1−1e−θ(x1+wx1)dx1.

Direct application of Lemma 3 implies that

12 ∫0∞x1−1e−θ(x1+wx1)dx1=2K0(2θw).

The result of the theorem follows by combining equations (11) and (12). □

By using the mathematical software MATLAB, we can present numerically in Figure 2, the graphs of pdfs of Y and V defined by (6) and (9), respectively, for different values of parameters. The graph of W will be given later.

Figure 2.

PDF of Y for the choices of α, β, θ and fixed parameter γ = 2 (a) and PDF of V for the choices of α, β and fixed parameters θ = 4, γ = 3 (b).

3.3. Approximations of the product

In this section, we derive the approximate distribution of the product W = X₁X₂. It is clear that the random variable P := 1/(1 + W) has support in the interval (0, 1). For this reason, in what follows we shall account for the approximate distribution by the beta distribution with parameters a, b > 0 (say P~dβ(a,  b)) and pdf

19 fP(y)=Γ(a+b)Γ(a)Γ(b)ya−1(1−y)b−1.

The idea of approximating distributions including complicated formulae with the beta distribution is frequently tackled in the statistics literature, based on the interesting work of Dasgupta [1] (see also [4,11,16] and [21]). Note that recently, Ghorbel [6] developed this single beta distribution of the form (19) for the bivariate beta distribution. Subsequently, this transformation was opted for owing to its popularity, simplicity and its association with the derived distribution. The basic reason for doing this doesn’t lie in the fact that the calculation of equation (10) cannot be handled; it rather resides in providing a simple approximation in terms of the beta distribution so that one can use the known procedures for inference, prediction, etc. Furthermore, the presented approximation seems to be useful especially to practitioners who not only avoid the use of the modified Bessel function of the second kind but also regard the beta distribution as widely accessible in standard statistical packages.

The choice of the beta parameters a and b is enacted using the method of moments. The first two moments of P can be written as

E(P)=aa+bandE(P2)=a(a+1)(a+b)(a+b+1).

Solving simultaneously the two expressions above, it is straightforward to find that

20 a=E(P)E(P)−E(P2)E(P2)−E2(P)

and

21 b=(1−E(P))E(P)−E(P2)E(P2)−E2(P).

where the two moments E(P) and E(P²) can be computed numerically by evaluating the two integrals

I1=∫0∞fW(w)1+wdwandI2=∫0∞fW(w)(1+w)2dw

for the given parameters α, β, γ and θ, where f_W (w) is presented by (10). To demonstrate the closeness of the proposed approximation, we use equations (20) and (21) with eight selected values of the parameters (α, β, γ, θ). Both integrals as well as the corresponding estimated values of the parameters a and b are depicted in Table 2.

Table 2.

Estimates of (I₁, I₂, a, b) for selected (α, β, γ, θ).

α	β	γ	θ	I1	I2	a	b
0.5	1.5	0.5	2.2	0.8662	0.7807	2.4314	0.3757
1.5	1.5	1.5	6	0.9087	0.8339	8.2760	0.8318
4	1.5	4	1.5	0.1269	0.0234	1.8010	12.3888
5	1.5	3	3	0.2700	0.0886	3.1271	8.4546
4	1.5	0.5	3.2	0.4106	0.1960	3.2063	4.6032
2	1.5	2	2	0.4666	0.2650	1.9894	2.2745
3.5	1.5	0.5	7	0.7849	0.6309	8.1810	2.2418
1	1.5	2	2	0.6835	0.5200	2.1181	0.9808

The corresponding graphics are outlined in Figure 3, demonstrating comparison between exact and approximate densities of P = 1/(1 + W). The exact pdf corresponds to the solid curve and approximate pdf stands for the broken curve. It is evident that the approximate density is quite close to the exact density. Note that the exact pdf of P is stated by (with 0 < p < 1)

fP(P)=2θ2α+1(1−p)α−1(θ(βθ+βγ−θ)(1−p)p−1+α2)β(γ+θ)Γ2(α+1)pα+1K0(2θ1−pp),

where K₀(·) is interpreted as the limit (as ν → 0) of the modified Bessel function of the second kind K_ν(x). Similar findings were recorded when this exercise was repeated for many other combinations of (α, β, γ, θ). Finally, It is worth noting that the numerical results are obtained by the use of besselk function in mathematical software MATLAB.

Figure 3.

The exact pdf of P (solid curve) and its approximated function (broken curve) for β = 1.5 and (1): (α, γ, θ) = (0.5, 0.5, 2.2) ; (2): (α, γ, θ) = (1.5, 1.5, 6) ; (3): (α, γ, θ) = (4, 4, 1.5) ; (4): (α, γ, θ) = (5, 3, 3); (5): (α, γ, θ) = (4, 0.5, 3.2); (6): (α, γ, θ) = (2, 2, 2); (7): (α, γ, θ) = (3.5, 0.5, 7); (8): (α, γ, θ) = (1, 2, 2).

Afterwards, based on the performed relationship in the first approximation, we present another direct approximation of W by the inverted beta distribution of the two parameters p > 0 and q > 0 say β′ (p, q), defined by the pdf

22 Γ(a+b)Γ(a)Γ(b)xp−1(1+x)−(p+q),x>0.

From this, we use the method of moments to determine the parameters p and q. The first two moments of W become

E(W)=qp−1andE(W2)=q(q+1)(p−1)(p−2).

After some algebraic manipulation, it is straightforward to verify that

23 p=1+E(W2)−E(W)E(W2)−E2(W)

and

24 q=E(W)E(W2)−E(W)E(W2)−E2(W),

where both moments E(W) and E(W²) follow from (16) by taking n = 1 and n = 2, respectively.

In order to corroborate the robustness as well as goodness of fit of the second approximation, we use (23) and (24) with given selected values of the parameters (α, β, γ, θ). The corresponding estimated values of the parameters p and q are summarized in Table 3. Afterwards, the notion of robustness is assessed by comparing the exact and the approximated pdfs of W as exhibited in (10) and (22). These comparisons are illustrated in Figure 4.

Table 3.

Estimates of (a, b) for selected (α, β, γ, θ)

α	β	γ	θ	p	q
1	1.5	0.5	1.7	1.9208	0.7820
0.5	1.5	0.5	2	1.3898	0.1153
2	1.5	0.5	2.5	1.9962	0.9918
0.5	1.5	1	1	2.1309	1.7906
2.5	1.5	2.5	2	2.7739	4.6442
4	1.5	1.5	4	2.4129	1.8223
2	1.5	0.5	2	2.3193	2.0888
1.5	1.5	0.7	2	2.0428	1.1144

Figure 4.

The exact pdf of W (solid curve) and its approximated function (broken curve) for β = 1.5 and (1): (α, γ, θ) = (1, 0.5, 1.7); (2): (α, γ, θ) = (0.5, 0.5, 2); (3): (α, γ, θ) = (2, 0.5, 2.5) ; (4): (α, γ, θ) = (0.5, 1, 1) ; (5): (α, γ, θ) = (2.5, 2.5, 2) ; (6): (α, γ, θ) = (4, 1.5, 4); (7): (α, γ, θ) = (2, 0.5, 2) ; (8): (α, γ, θ) = (1.5, 0.7, 2).

4. Entropies

The entropy is a crucial concept in such areas such as classical thermodynamics, where it was first recognized, statistical mechanics and physics, information theory and many other fields. It stands for an excellent tool to quantify the amount of uncertainty involved in the value of a random variable or the outcome of a random process. Several measures of entropy have been examined and compared in literature. The simplest known entropy is the Shannon entropy (see [20]). Let (X,  ℬ,  P) be a probability space. Consider a pdf f associated with P, dominated by σ-finite measure μ on X. Thus, the Shannon entropy is defined by

One of the main extensions of the Shannon entropy was established by Rényi (see [19]). This generalized entropy measure is expressed in terms of

where

Note that the Shannon entropy is the particular case of the Rényi entropy for μ ↑ 1. The calculations of these entropies involve the hypergeometric function identified as

where (i)_k := i(i + 1) … (i + k – 1) denotes the ascending factorial.

Before deriving these entropies, we need the following lemma:

5. Estimation by Expectation-Maximization (E.M) Algorithm

The Expectation-Maximization (E.M) algorithm (see [2]) is a prominent way of finding maximum-likelihood estimates (MLEs) for model parameters when data are incomplete. It is an itetarive method to approximate the maximum likelihood function. Each iteration of the E.M algorithm involves two steps. In the first one (E-Step), we compute the conditional expectation denoted by ℚ(Θ|Θ(j)), where Θ^(j) = (α^(j), β^(j), γ^(j), θ^(j)) is the estimated parameters vector of Θ = (α, β, γ, θ) at the j^th iteration. In the second step of the E.M alghorithm (M-Step), we maximize ℚ(Θ|Θ(j)) with respect to the vector of parameters Θ.

In this section, we discuss the problem of computing MLEs of the unknown parameters of the BGGL distribution for a given random sample of size n, of the form Λ = {(x₁₁, x₁₂),..., (x_n1, x_n2 )}· Using the mixture density representation of the BGGL distribution given by (4), the incomplete likelihood function can be written as

where Θ₁ = (α, θ) and Θ₂ = (α, θ, γ).

To complete the data of this function, we associate a discrete random vector Δ, where Δ follows a bivariate Bernoulli distribution with parameters p=β−1β and 1−p=1β. This implies that

where δ_i ∈ {0, 1}². The complete samples are as follows:

Resting on the complete observations, the complete likelihood function is

Thus, the log-likelihood function becomes

Now, we compute the conditional expectation ℚ(Θ|Θ(j)) of the complete log likelihood function, resting upon the observed data and the current value Θ^(j) = (α^(j), β^(j), γ^(j), θ^(j)) at the j^th iteration. We obtain