I want to know the theoretical distribution of a mixture of exponential distributions whose rate parameters are distributed according to a gamma distribution:
$$ y\sim\text{Exp}(\theta), \quad\text{where}\quad \theta\sim\Gamma(r, \beta).$$
Specifically, I am looking at $r=4$ and $\beta=2$.
For all $y \ge 0$ the value of the survival function of $Y$ is
$$S_Y(y\mid\theta) = \Pr(Y \gt y\mid \theta) = \exp(-y\theta)$$
and, taking $\beta$ to be a rate parameter for the Gamma distribution, the probability density function of $\theta$ is proportional to
$$f_\theta(t) \propto t^{r-1}\exp(-\beta t).$$
Consequently the survival function of the mixture distribution is
$$S(y) \propto \int_0^\infty S_Y(y\mid t) f_\theta(t)\,\mathrm{d}t = \int_0^\infty t^{r-1}\exp(-(\beta+y)t)\,\mathrm{d}t.$$
Substituting $u = (\beta+y) t$ gives, with no calculation,
$$S(y) \propto \int_0^\infty \left(\frac{u}{\beta+y}\right)^{r-1}\exp(-u)\,\mathrm{d}\left(\frac{u}{\beta+y}\right) = (\beta+y)^{-r}\int_0^\infty u^{r-1}\exp(-u)\,\mathrm{d}u$$
which is proportional to $(\beta+y)^{-r}.$
The axiom of total probability asserts $S(0)=1$ from which we obtain the implicit constant, giving
$$S(y) = \beta^r\,(\beta+y)^{-r} = \left(1 + \frac{y}{\beta}\right)^{-r},$$
thereby exhibiting $\beta$ as a scale parameter for this mixture variable.
This is a particular kind of Beta prime distribution, sometimes termed a Lomax distribution. (up to scale).
If $\beta$ is intended to be a scale parameter for the Gamma distribution rather than a rate parameter, replace $\beta$ everywhere by $1/\beta$ and, at the end, interpret it as a rate parameter for the mixture variable.
This distribution tends to be positively skewed, so it's better to view the distribution of $\log Y.$ Here are the empirical (black) and theoretical (red) distributions for a simulation of $10^4$ independent realizations of $Y$ (where $r=4$ and $\beta=2$ as in the question):
The agreement is perfect.
The R code that generated this figure illustrates how to code the survival function (S), the density function (its derivative f), and how to simulate from this compound distribution by first producing realizations of $\theta$ and then, for each of them, generating a realization of $Y$ conditional on $\theta.$
S <- function(y, r, beta) 1 / (1 + y/beta)^r # Survival
f <- function(y, r, beta) r * beta^r / (beta + y)^(r+1) # Density
#
# Specify the parameters and simulation size.
#
beta <- 2
r <- 4
n.sim <- 1e4
#
# The simulation.
#
theta <- rgamma(n.sim, r, rate=beta)
y <- rgamma(n.sim, 1, theta)
#
# The plots.
#
par(mfrow=c(1,2))
plot(ecdf(log(y)), xlab=expression(log(y)), ylab="Probability",
main=expression(1-S[Y](y)))
curve(1 - S(exp(y), r, beta), xname="y", add=TRUE, col="Red", lwd=2)
hist(log(y), freq=FALSE, breaks=30, col="#f0f0f0", xlab=expression(log(y)))
curve(f(exp(y), r, beta) * exp(y), xname="y", add=TRUE, col="Red", lwd=2)
par(mfrow=c(1,1))
Wikipedia, under the "Examples" section, informs us that:
Compounding an exponential distribution with its rate parameter distributed according to a gamma distribution yields a Lomax distribution [9].
The reference [9] is to Johnson, N. L.; Kotz, S.; Balakrishnan, N. (1994). "20 Pareto distributions". Continuous univariate distributions. 1 (2nd ed.). New York: Wiley. p. 573
