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!     
 A: 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.
A: 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

