Bivariate Gamma Distribution Overview

• The bivariate gamma distribution is a continuous model for two interdependent nonnegative variables, defined by a positive shape parameter, a positive-definite covariance matrix, and a non-centrality component.
• It features series expansions for both the joint PDF and CDF, allowing for tractable computation and clearer insights into correlated gamma variates through integral representations.
• The distribution is intrinsically linked to non-central Wishart matrices, highlighting its significance in multivariate statistical analysis and applications to correlated gamma-type data.

The bi-variate gamma distribution, also frequently referred to as the two-variate (non-central) multivariate gamma law, generalizes the univariate gamma distribution to two dependent, non-negative random variables. The construction is naturally linked to the diagonal elements of non-central Wishart matrices and is parameterized by a positive shape parameter, a positive-definite scale/covariance matrix, and a symmetric non-centrality matrix. This distribution is of fundamental importance both as a mathematical object and as a tool in multivariate statistical analysis, distributional theory, and the study of correlated gamma-type random variables (Royen, 2016).

1. Definition and Parameterization

Let $\vec{X} = (X_1, X_2)^{\prime}$ be a random vector characterized by:

• Shape parameter: $\alpha > 0$, sometimes expressed in the Wishart context as degrees of freedom $\nu = 2\alpha$.
• Scale (covariance) matrix: $\Sigma = [\sigma_{ij}]_{2\times2} > 0$. It admits the factorization $\Sigma = W^{-2} R W^{-2}$, where $W = \operatorname{diag}(w_1, w_2)$ with $w_i = \sigma_{ii}^{-1/2}$, and $R = (r_{ij})$ is a correlation matrix such that $R_{ii}=1$, $R_{12}=R_{21}=\rho$, $|\rho|<1$.
• Non-centrality matrix: $\Delta = [\delta_{ij}]_{2\times2} \geq 0$, symmetric positive semi-definite, with frequent consideration of rank-one forms $\Delta = d d^\prime$ for some $d = (d_1, d_2)^\prime$.

The Laplace transform of the joint distribution is given by: $$\mathcal{L}_{X}(t_1, t_2) = E\bigl[e^{-t_1 X_1 - t_2 X_2}\bigr] = |I + R T|^{-\alpha} \exp\left( -\mathrm{tr}\left[ T (I + R T)^{-1} \Delta \right]\right)$$ where $T = \operatorname{diag}(t_1, t_2) \geq 0$. This characterizes the distribution completely, extending the univariate gamma and non-central chi-square families (Royen, 2016).

2. Series Expansions for the Joint Density and Distribution Functions

For the bivariate gamma law, the CDF and PDF admit absolutely convergent double-series expansions. Introducing auxiliary quantities:

• $W = \operatorname{diag}(w_1, w_2)$,
• $B = W R W - I$,
• $D = W^{-1} \Sigma^{-1} \Delta \Sigma^{-1} W^{-1}$,
• $g_\alpha(x) = x^{\alpha-1} e^{-x}/\Gamma(\alpha)$ (central gamma density),
• $P(\alpha, x) = \gamma(\alpha, x)/\Gamma(\alpha)$ (regularized lower incomplete gamma),

the joint CDF is

$$G(x_1, x_2) = \exp\!\left(-\operatorname{tr}\left(\Delta \Sigma^{-1}\right)\right) \sum_{n_1=0}^\infty \sum_{n_2=0}^\infty \frac{d(n_1, n_2)}{n_1! n_2!} \; P(\alpha+n_1; w_1^2 x_1) \, P(\alpha+n_2; w_2^2 x_2)$$

where, for $p=2$,

$$d(n_1, n_2) = \begin{cases} \frac{(B_{12})^{n_1+n_2} (n_1 + n_2)!}{n_1! \, n_2!}, & n_1 + n_2 \text{ even}, \ 0, & n_1 + n_2 \text{ odd}. \end{cases}$$

Differentiating term-by-term yields the joint PDF: $$f(x_1, x_2) = \exp\left(-\operatorname{tr}(\Delta \Sigma^{-1})\right) \sum_{n_1, n_2 \geq 0} \frac{d(n_1, n_2)}{n_1! n_2!} g_{\alpha + n_1}(w_1^2 x_1) \, g_{\alpha + n_2}(w_2^2 x_2)$$ The structure of the series reflects how correlation and non-centrality couple the marginal gamma laws (Royen, 2016).

3. Integral Representations

The cumulative distribution function also admits a single-integral (Fourier-type) representation by virtue of every $2 \times 2$ correlation matrix being "one-factorial." Define:

• $$\rho = \sigma_{12} / (\sigma_{11}\sigma_{22})^{1/2}$$,
• $y = \rho e^{i\varphi}$, $\bar y = \rho e^{-i\varphi}$,
• $$d_{ij} = [W^{-1} \Sigma^{-1} \Delta \Sigma^{-1} W^{-1}]_{ij}$$,
• $\Lambda = \operatorname{tr}(\Delta \Sigma^{-1})$,
• $$\widetilde{P}(\alpha; u, v) = \exp(-v) \sum_{k=0}^\infty \frac{v^k}{k!} P(\alpha+k; u)$$,

Then, the CDF has the form

$$G(x_1, x_2) = (1-\rho^2)^{\alpha} e^{-\Lambda} \frac{1}{2\pi} \int_0^{2\pi} \exp\left(d_{11} y + d_{22} \bar y\right) (1-\rho^2 y)^{-\alpha} \,\widetilde{P}\left(\alpha; w_2^2 x_2 (1-\rho^2 \bar y), d_{22} \bar y\right) \, d\varphi$$

This integral is highly tractable, converges rapidly, and reveals the analytic structure imposed by both the correlation coefficient and the non-centrality matrix. Differentiation under the integral sign provides an analogous formula for the joint PDF (Royen, 2016).

4. Influence of Non-Centrality and Correlation

The non-centrality matrix $\Delta$ appears in the formulae solely through the scalar $\Lambda = \operatorname{tr}(\Sigma^{-1}\Delta)$ and the entries $d_{ij}$ of $D = W^{-1} \Sigma^{-1} \Delta \Sigma^{-1} W^{-1}$. The correlation $\rho$ influences (i) the prefactor $(1-\rho^2)^{\alpha}$, (ii) the power-law and exponential terms in the integrand (via $y = \rho e^{i\varphi}$), and (iii) the scale factors $w_i = \sigma_{ii}^{-1/2}$.

In the case of rank-one non-centrality, $\Delta = dd'$, the exponent in the integrand specializes to $d_{11}y + d_{22}\bar y + 2\,\Re\{d_{12}y\}$, further simplifying computational realization. This form offers insight into the coupling mechanism between marginals as driven jointly by non-centrality and correlation, particularly for practical simulation or evaluation (Royen, 2016).

5. Special Cases and Computational Considerations

Distinct important specializations include:

• Central case ($\Delta = 0$): $\Lambda=0$, $d_{ij}=0$, eliminating the exponential term; the result reduces to the Appell $F_1$ (or Gauss hypergeometric) representation for the central bivariate gamma/chi-square law.
• Equal-scale case ($\sigma_{11} = \sigma_{22}$): $w_1 = w_2$, yielding a symmetric appearance of $y$ in $x_1$ and $x_2$.
• Independent marginals ($\rho=0$): The integral collapses to the product $P(\alpha; w_1^2 x_1) P(\alpha; w_2^2 x_2)$, i.e., independence in the joint law.
• Rank-one non-centrality ($\Delta = dd'$): The CDF admits the simpler single-sum expansion:

$$G(x_1, x_2) = e^{-\|W^{-1}d\|^2} \sum_{k=0}^{\infty} \frac{\left(\|W^{-1}d\|^2\right)^k}{k!} P(\alpha+k; w_1^2 x_1) P(\alpha+k; w_2^2 x_2)$$

In practical calculations, the integral over $\varphi \in [0, 2\pi]$ is truncated or mapped to $[-\pi, \pi]$ and computed using Gaussian or Fourier-type quadrature; convergence is rapid, since all singularities are exterior to the integration path $|y|=|\rho|<1$. This structure affords both (i) a series expansion in non-centrality, positive and real, and (ii) a compact, analytically transparent, single-integral representation (Royen, 2016).

6. Connections to Wishart Distributions and Broader Significance

The bi-variate gamma law emerges as the diagonal marginal distribution of a $2 \times \nu$ non-central Wishart matrix $W_2(\nu, \Sigma, \Delta)$ with $\nu = 2\alpha$. The general framework extends to higher-dimensional (p-variate) settings, with $p-1$-dimensional integral representations for the CDF. For $p=2$, these reduce to the explicit, computationally tractable forms described. This constructs a systematic link between marginal gamma-type variables with arbitrary non-centrality and correlation, bridging chi-square, Wishart, and multivariate gamma-analytic methods. Alternative formulas exist in the case of "one-factorial" correlation matrices, leading to further simplifications for $p=3$ and certain structured correlation models.

The principal results and representations outlined here, as developed in the work of Royen (Royen, 2016), are foundational for computation, theory, and application of correlated gamma variates in multivariate statistics and probabilistic analysis.