# What is the distribution of a mixture of exponential distributions whose rate parameters follow a gamma distribution?

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

• The question appears clear enough to me: given $r$ and $\beta$, let $\theta$ follow a gamma distribution with parameters $r$ and $\beta$ (whether these are scale or rate params doesn't matter for now), denoted $\theta\sim\Gamma(r, \beta)$. Next, $y$ is exponentially distributed with parameter $\theta$, denoted $y\sim\text{Exp}(\theta)$. That is, we have a compound distribution, AKA a mixture. Q: What is the unconditional distribution of $y$ with parameters $r$ and $\beta$? May 1 '20 at 5:29
• At least that is what the original question seems to have asked, and I have tried to clarify this specific question with my edit and answer it. I will vote to reopen. Can someone clarify to me where this question is unclear? May 1 '20 at 5:30
• @StephanKolassa: Well, $\theta$ wasn't defined; it could be the mean rather than the rate. The variable $n$ & its role were also unexplained. May 2 '20 at 8:57
• @Scortchi-ReinstateMonica: true. But converting between a mean and a rate formulation for the exponential is not overly hard. I assumed $n$ referred to a sample size and removed it in my edit. May 2 '20 at 9:11
• @StephanKolassa: (1) The compound distribution would be something else if it were the means distributed according to a gamma distribution. (2) Quite possibly, but why mention it then? I think your edits do a fine job of clarification, but I can see why the q. was closed in its original version. May 2 '20 at 9:19

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 .

The reference  is to Johnson, N. L.; Kotz, S.; Balakrishnan, N. (1994). "20 Pareto distributions". Continuous univariate distributions. 1 (2nd ed.). New York: Wiley. p. 573

• Yuh, but is that really the question?
– Carl
May 1 '20 at 4:10