# Product of Gamma by Beta rv

If $X$ has a beta distribution $\beta(\alpha,b)$, $Y$ has a gamma distribution $\Gamma (K,\theta)$ and $X$ is independent of $Y$. What is the distribution of the product $P=XY$ .

Thanks!

• The PDF can be obtained as a linear combination of confluent hypergeometric functions. If that form is not useful to you, then please edit your question to include a description of what you need to know about the product distribution. – whuber Jun 23 '14 at 13:13
• Could you please provide me with some references?(I need the PDF or the MGF of this product) – tam Jun 23 '14 at 13:33
• Brute-force application of a formula for the ratio PDF will do it. – whuber Jun 23 '14 at 13:47
• @whuber if you please, I need the detailed calculation (it's complicated..) – tam Jun 24 '14 at 9:25
• Actually, it looks pretty simple: expand the integrand $f_X(x)f_Y(z/x)/x$ as a Maclaurin series in $z$ and integrate term-by-term from $0$ to $1$. The result clearly follows the pattern of development of hypergeometric functions, so it remains only to determine their parameters, which can be done by inspection. – whuber Jun 24 '14 at 13:56

Given:

• $X \sim \text{Beta}(a,b)$ with pdf $f(x)$:

• $Y \sim \text{Gamma}(k,\theta)$ with pdf $g(y)$:

Solution: Then, the pdf of the product $Z = X Y$ can be automatically derived via:

where I am using the TransformProduct function from the mathStatica package for Mathematica, and where Hypergeometric1F1 denotes the Kummer confluent hypergeometric function: http://reference.wolfram.com/mathematica/ref/Hypergeometric1F1.html

This formulation works nicely, except for certain combinations of integer values of the parameters (indeterminate - please see discussion below). [If say $a = 4$ and $k = 3$, just enter $k$ as 3.0000001 and it will side-step the issue.]

Quick Monte Carlo check

It is always a good idea to check symbolic solutions with Monte Carlo methods. Here is a quick comparison of the exact theoretical solution derived above (dashed RED curve) against an empirical Monte Carlo simulation of the pdf of the product (squiggly BLUE), when ${a = 3, b = 6, k = 2.2, \theta = 5}$

All done.

• Thank you @wolfies for your help, I really appreciate it! – tam Jun 23 '14 at 14:33
• +1 The reason you run into trouble for certain values of the parameters is that Gamma has poles where $a-k$ and $a+b-k$ are negative integers and the hypergeometric functions have corresponding poles and zeros to cancel them, leading to finite values. You can sneak up to the correct value by taking limits (which Mathematica will do); they will involve derivatives of the confluent hypergeometric function with respect to its parameters (the first two arguments). – whuber Jun 23 '14 at 14:48
• @whuber Many thanks -> Mathematica spat out the condition $a>k$. However, the problem with integer-valued $k$ occurs even when $a > k$; moreover, the condition $a>k$ appeared superfluous [checked with Monte Carlo]. – wolfies Jun 23 '14 at 14:53
• It's not superfluous; in fact, it's not enough. If (a) any of the Gamma arguments $a+b-k$ or $k-a$ or $a-k$ are negative integers or (b) any of the second arguments of the hypergeometric function $1-k+a$ and $1+k-a$ are non-positive integers, the result is indeterminate. That leaves only the possibility $k=a$ when $k$ and $a$ are integral! – whuber Jun 23 '14 at 15:01
• Makes sense. Thanks. Still, the condition $a>k$ is superfluous or too general, when the required constraint is that: either $a$ or $k$ must not be an integer, or similar. – wolfies Jun 23 '14 at 15:10

This distribution is called the Gamma-Inverse Beta distribution in this paper. It is available in the R package brr.

nsims <- 1e6
alpha <- 3
beta <- 5
K <- 6
theta <- 4
sims <- rgamma(nsims, shape = K, rate = theta) * rbeta(nsims, alpha, beta)

plot(density(sims, to=3))
curve(brr::dGIB(x, K, beta, alpha, theta), # note that alpha and beta are swapped
add = TRUE, col = "red", lwd = 2, lty = "dashed")


Here is its density function:

> brr::dGIB
function (x, a, alpha, beta, rho)
{
exp(lnpoch(beta, alpha) - lngamma(a)) * rho^a * x^(a - 1) *
exp(-rho * x) * hyperg_U(alpha, a - beta + 1, rho * x)
}


$$\frac{{(\beta)}_\alpha}{\Gamma(a)}\rho^a x^{a-1} \exp(-\rho x) U(\alpha, a-\beta+1, \rho x), \qquad x > 0$$ where $$U$$ is the Tricomi hypergeometric function.

There's a generalization of this distribution to random matrices. It is a type II confluent hypergeometric function distribution of kind two (see Gupta & Nagar's book Matrix variate distributions).

library(matrixsampling)
sims2 <- rmatrixCHIIkind2(20000, nu = K, alpha = beta,
beta = K+1-alpha, theta = 1/theta, p = 1)
lines(density(sims2, to=3), col = "orange", lwd = 2, lty = "dashed")


The mgf is given in Gupta & Nagar's book. With your notations, this gives $$\textrm{mgf}(t) = {}_2\!F_1(\alpha, K, \alpha+\beta, t/\theta).$$

library(gsl)
hyperg_2F1(alpha, K, beta+alpha, 0.2/theta)
## [1] 1.121871
mean(exp(0.2*sims))
## [1] 1.121846