A Bivariate Gamma Distribution Whose Marginals are Finite Mixtures of Gamma Distributions

In this article a new bivariate distribution, whose both the marginals are finite mixtures of gamma distributions, has been defined. Several of its properties such moments, correlation coefficients, measure of skewness, moment generating function, Rényi and Shannon entropies have been derived. Simulation study has been conducted to evaluate the performance of maximum likelihood method.


Introduction
The univariate gamma distribution is one of the most commonly used statistical distributions to analyze skewed data in many disciplines and has been studied extensively in scientific literature. The chi-square distribution, which is of utmost importance in statistical inference, is a special case of gamma distribution. Probability distributions such as exponential and Erlang are also special cases of the gamma distribution. Several univariate generalizations and variants of gamma distribution have also been developed and applied in various areas.
The univariate gamma distribution has been generalized to the bivariate case in many different ways and many forms of bivariate gamma distribution are available. Several techniques to generate bivariate distributions have also been proposed in the scientific literature, e.g., see Balakrishnan and Lai [1], Mardia [12], and Zhang and Singh [30].
Bivariate gamma distributions have found useful applications in many areas. They have been used for representing joint probabilistic properties of multivariate hydrological events such as floods and storms or in the modeling of rainfall at two nearby rain gauges, data obtained from rainmaking experiments, the dependence between annual streamflow and aerial precipitation, wind gust modeling (Smith and Adelfang [25], Smith, Adelfang, and Tubbs [26]), and the dependence between rainfall and runoff (see Nadarajah and Gupta [16], Nadarajah [14,15] and references therein). For an interesting review of bivariate gamma distributions for hydrological application, the reader is referred to Yue, Ouarda, and Bobée [29] and Zhang and Singh [30].
Nadarajah [13] has listed a number of bivariate gamma distributions such as McKay's bivariate gamma distribution, Dussauchoy and Berland's bivariate gamma distribution, Cherians bivariate gamma distribution, Arnold and Strauss' bivariate gamma distribution, Becker and Roux's bivariate gamma distribution, and Smith and Adelfang's bivariate gamma distribution.

The bivariate gamma distribution
The random variables X 1 and X 2 are said to have a bivariate gamma distribution with parameters α, β and k, denoted by (X 1 , X 2 ) ∼ BGa(α, β, k), if their joint density is given by where x 1 > 0, x 2 > 0, α > 0, β > 0, k ∈ N 0 and C(α, β, k) is the normalizing constant. By integrating the joint density of X 1 and X 2 over its support set, the normalizing constant is derived as Now, expanding (x 1 + x 2 ) k using binomial theorem and integrating x 1 and x 2 , we obtain

952
A BIVARIATE GAMMA DISTRIBUTION Finally, using Lemma A.1, we get and An alternative way to compute (2) is to substitute s = x 1 + x 2 and r = x 1 /(x 1 + x 2 ) and integrate s and r by using gamma and beta integrals. Since this approach works for all k > 0, we will use it to compute Shannon entropy.
Let us now briefly discuss the shape of (1). The first order derivatives of ln f (x 1 , x 2 ; α, β, k) with respect to x 1 and x 2 are and respectively. Setting (5) and (6) to zero, the only stationary point of (1) is obtained as (5) and (6), we get and Further, from (7), (8) and (9), we get

Now, observe that
Stat., Optim. Inf. Comput. Vol. 8, December 2020 2 > 0, f x1x1 (a, a) < 0 and f x2x2 (a, a) < 0 and therefore (a, a) is a maximum point. 2 < 0, and therefore (a, a) is a saddle point. Figure 1 illustrates the shape of the pdf (1) for selected values of α and β and k. It can easily be observed that (X 1 , X 2 ) and (X 2 , X 1 ) are identically distributed and hence X 1 and X 2 are exchangeable.
A distribution is said to be negatively likelihood ratio dependent if the density f ( for all x 1 > x * 1 and x 2 > x * 2 (Lehmann [9], Tong [28]). One can check that the bivariate distribution defined by the density (1) is negatively likelihood ratio dependent.
By integrating x 2 in (1) the marginal density of X 1 is obtained as Substituting x 2 /x 1 = z in (10), the marginal density of X 1 is rewritten as Now, Writing (1 + z) k using binomial theorem and integrating z in (11), the marginal density of X 1 is derived as Likewise, the marginal density of X 2 is obtained as Thus, the marginal density of X i is a finite mixture of gamma densities. Figure 2 shows some plots of the marginal density of X 1 for β = 2, k = 0, 1, . . . , 20 and some values of α. (11), one gets Now, writing  in (14) and integrating u, the density f X1 (x 1 ), in series involving generalized Laguerre polynomials, is derived as where L (α) j (·) is the generalized Laguerre polynomial (see Appendix for the definition). From the joint pdf (1) and the marginal density of X given in (12), the conditional pdf of X 2 given X 1 = x 1 is given by Also, the conditional pdf of X 1 given X 2 = x 2 is given by

Moments
By definition where the last line has been obtained by using beta and gamma integrals. Finally, simplifying the above expression, we get Further, substituting appropriately in the above expression, one gets Further, variances, covariances, correlation and several higher central moments are derived as

Moment Generating Function
By definition, the joint mgf of X 1 and X 2 is given by

958
A BIVARIATE GAMMA DISTRIBUTION in (19) and integrating r, we get where the last line has been obtained by using the integral representation of the Gauss hypergeometric function given in (A.1). Finally, substituting for C(α, β, k) and simplifying, we get ) .
which is the mgf of a gamma random variable with shape parameter 2α + k and scale parameter β.

Entropies
In this section, exact forms of Rényi and Shannon entropies are derived for the bivariate gamma distribution defined in this article. Let (X , B, P) be a probability space. Consider a pdf f associated with P, dominated by σ−finite measure µ on X . Denote by H SH (f ) the well-known Shannon entropy introduced in Shannon [24]. It is define by One of the main extensions of the Shannon entropy was defined by Rényi [21]. This generalized entropy measure is given by where The additional parameter η is used to describe complex behavior in probability models and the associated process under study. Rényi entropy is monotonically decreasing in η, while Shannon entropy (21) is obtained from (22) for η ↑ 1. For details see Nadarajah and Zografos [17], Zografos and Nadarajah [32] and Zografos [31].

Theorem 5.1
For the bivariate gamma distribution defined by the pdf (1), the Rényi and the Shannon entropies are given by and

Proof
For η > 0 and η ̸ = 1, using the joint density of X 1 and X 2 given by (1), we have where the last line has been obtained by substituting s = x 1 + x 2 and r = x 1 /(x 1 + x 2 ). Finally, evaluating above integrals by using gamma and beta integrals and simplifying the resulting expression, we get Now, taking logarithm of G(η) and using (22) we get H R (η, f ). The Shannon entropy is obtained from H R (η, f ) by taking η ↑ 1 and using L'Hopital's rule.

Sum, Quotient and Product
In this section we derive the distributions of X 1 + X 2 , X 1 /(X 1 + X 2 ), X 1 X 2 , and X 1 /X 2 when X 1 and X 2 follow a bivariate gamma distribution defined in (1).

Theorem 6.1
Let (X 1 , X 2 ) ∼ BGa(α, β, k), and define R = X 1 /(X 1 + X 2 ) and S = X 1 + X 2 . Then, R and S are independent, the distribution of R is beta with both the parameters α and the distribution of S is gamma with shape parameter 2α + k and scale parameter β.
Proof Substituting x 1 = rs and x 2 = s(1 − r) with the Jacobian J(x 1 , x 2 → r, s) = s, in the joint density of X 1 and X 2 , we obtain the joint density of R and S as where 0 < r < 1 and s > 0. Now, from (23), the desired result is obtained.
Corollary 6.1.1 Both X 1 /X 2 and X 2 /X 1 have inverted beta distribution with parameters α and α.

Theorem 6.2
Let (X 1 , X 2 ) ∼ BGa(α, β, k), and define P = X 1 X 2 . Then, the density of P is given by Proof Transforming X 1 = X and P = X 1 X 2 with the Jacobian J(x 1 , x 2 → p) = 1/x in the joint density of X 1 and X 2 and integrating x, we obtain the density of P as Now, using the integral (Gradshteyn and Ryzhik [4, Eq. 3.471.9]), where K ν is the modified Bessel function of the second kind, we obtain the desired result.
Next two theorems deal with bivariate distributions of (X 1 /Y, X 2 /Y ) and (X 1 /U,

Theorem 6.3
Let (X 1 , X 2 ) ∼ BGa(α, β, k), and Y ∼ Ga(ν, β) be independent. Then, the joint density of Z 1 = X 1 /Y and Z 2 = X 2 /Y is given by Proof Transforming X 1 = Z 1 Y and X 2 = Z 2 Y with the Jacobian J(x 1 , x 2 , → z 1 , z 2 ) = y 2 in the join density of (X 1 , X 2 ) and Y , the joint density of (Z 1 , Z 2 ) and Y is obtained as Now, integrating y by using gamma integral, we get the desired result.
For k = 0, the variables X 1 , X 2 and Y are independent gamma with scale parameter β and therefore (X 1 /Y, X 2 /Y ) has a Dirichlet type 2 distribution.

Theorem 6.4
Let (X 1 , X 2 ) ∼ BGa(α, β, k), and U ∼ B(a, b) be independent. Then, the joint density of Z 1 = X 1 /U and Z 2 = X 2 /U is given by Proof Transforming X 1 = Z 1 U and X 2 = Z 2 U with the Jacobian J(x 1 , x 2 , → z 1 , z 2 ) = u 2 in the join density of (X 1 , X 2 ) and U , the joint density of (Z 1 , Z 2 ) and U is obtained as where z 1 > 0, z 2 > 0, and 0 < u < 1. Now, integrating u by using integral representation of confluent hypergeometric function given in (A.2), we get the desired result.

Estimation
Let (X 11 , X 12 ), · · · , (X n1 , X n2 ) be a random sample from BGa(α, β, k). The log-likelihood function, denoted by l(α, β), is given by Using the duplication formula for digamma function, namely, Further, For a given observation vector (x 1 , x 2 ), the Fisher information matrix for the bivariate distribution given by the density (1) is defined as and Thus, by solving numerically (25) and (26), the MLEs of α and β can be obtained.

Simulation
In this section a simulation study is conducted to evaluate the performance of maximum likelihood method. Samples of size n = 30, 50, 200, 500 from Equation (1)
The generating function of the generalized Laguerre polynomial is Finally, we define the gamma, beta type 1 and beta type 2 distributions. These definitions can be found in Johnson, Kotz and Balakrishnan [7].

Definition A.1
A random variable X is said to have a gamma distribution with parameters θ (> 0), κ (> 0), denoted by X ∼ Ga(κ, θ), if its pdf is given by Note that for θ = 1, the above distribution reduces to a standard gamma distribution and in this case we write X ∼ Ga(κ).

Definition A.2
A random variable X is said to have a beta type 1 distribution with parameters (a, b), a > 0, b > 0, denoted as X ∼ B1(a, b), if its pdf is given by where B(a, b) is the beta function.

Definition A.3
A random variable X is said to have a beta type 2 (inverted beta) distribution with parameters (a, b), denoted as X ∼ B2(a, b), a > 0, b > 0, if its pdf is given by