Generate random numbers with bivariate gamma distribution in R

Question

How to generate random numbers with bivariate gamma distribution. The density is:

F_{(X, Y)}(x, y) = α^p+qx^p-1(y-x)^q-1e^-αy / [Γ(p) Γ(q)], 𝕀_{0≤ x≤ y}

With y>x>0, α>0, p>0 and q>0.

I did not find any package on R that does this and nothing in literature.

Answer 1

This is straightforward:

1. Generate X~ Gamma(p,alpha) (alpha being the rate parameter in your formulation)

2. Generate W~ Gamma(q,alpha), independent of X

3. Calculate Y=X+W

4. (X,Y) have the required bivariate distribution.

in R (assuming p,q,alpha and n are already defined):

x <- rgamma(n,p,alpha)
y <- x + rgamma(n,q,alpha)

generates n values from the bivariate distribution with parameters p,q,alpha