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$ .


  • $\begingroup$ 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. $\endgroup$
    – whuber
    Jun 23 '14 at 13:13
  • $\begingroup$ Could you please provide me with some references?(I need the PDF or the MGF of this product) $\endgroup$
    – tam
    Jun 23 '14 at 13:33
  • $\begingroup$ Brute-force application of a formula for the ratio PDF will do it. $\endgroup$
    – whuber
    Jun 23 '14 at 13:47
  • $\begingroup$ @whuber if you please, I need the detailed calculation (it's complicated..) $\endgroup$
    – tam
    Jun 24 '14 at 9:25
  • $\begingroup$ 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. $\endgroup$
    – whuber
    Jun 24 '14 at 13:56


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

enter image description here

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

enter image description here

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

enter image description here

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}$

enter image description here

All done.

  • $\begingroup$ Thank you @wolfies for your help, I really appreciate it! $\endgroup$
    – tam
    Jun 23 '14 at 14:33
  • $\begingroup$ +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). $\endgroup$
    – whuber
    Jun 23 '14 at 14:48
  • $\begingroup$ @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]. $\endgroup$
    – wolfies
    Jun 23 '14 at 14:53
  • $\begingroup$ 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! $\endgroup$
    – whuber
    Jun 23 '14 at 15:01
  • $\begingroup$ 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. $\endgroup$
    – 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")

enter image description here

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).

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")

enter image description here

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). $$

hyperg_2F1(alpha, K, beta+alpha, 0.2/theta)
## [1] 1.121871
## [1] 1.121846

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.