Nested Gamma Distributions Let $a, c$ and $d$ be positive real numbers. Let us also assume that $b \sim \Gamma (c, d) $. Next, let us define random variable $X \sim \Gamma (a, b)$. Then we would like to know the probability density function of $X$ given $b$.
I wanted to use somehow the following identity
$$ \mathbb{P} ( X \leq x ) = \int_{0}^{\infty} \mathbb{P} ( X \leq x | b = z) f_{b} (z) \ dz .  $$
However, I don't know whether I am expressing $\mathbb{P} ( X \leq x | b = z)$ correctly and I get some difficult integrals. I tried to differentiate the identity above with respect to $x$ to get the probability density function, but I got something resembling the probability density function of Beta distribution. I would appreciate any help or advice.
 A: A bit too long for a comment, but here's an idea. Denote a probability distribution function of some random variable $x$ conditional on some parameters $a,b$ as $f(x| a,b).$
Your variable $b$ is distributed as $f(b| c,d)$. Conditional on a particular value of $b$, your variable $x$ is distributed as $f(x| a, b).$ The distribution of $x$ unconditional on $b$ is then a weighted average  of $f(x|a,b)$ across all possible values of $b$:
$$ f(x| a,c,d) = \int_0^\infty f(b| c,d) f(x| a,b) db. $$
Provided $f(x|a,b)$ has the form
$$ f(x|a,b) = \frac{b^a}{\Gamma(a)}x^{a-1}e^{-b x},$$
your integral is then
$$f(x|a,c,d) = \int_0^\infty \frac{d^c}{\Gamma(c)}b^{c-1}e^{-d b} \frac{b^a}{\Gamma(a)}x^{a-1}e^{-b x} db \\
 = \frac{d^c x^{a-1}}{\Gamma(a)\Gamma(c)}\int_0^\infty b^{a+c-1}e^{-(d+x)b}db
\\
= \frac{d^c x^{a-1}}{\Gamma(a)\Gamma(c)}\frac{1}{(d+x)^{a+c}} \int_0^\infty u^{a+c-1}e^{-u}du
\\ = \frac{d^c x^{a-1}}{(d+x)^{a+c}}\frac{\Gamma(a+c)}{\Gamma(a)\Gamma(c)}
$$
or something along these lines. This is pretty close to a Beta distribution. Maybe it should be one! But I have to go to bed. Good luck lol
A: You want the marginal (unconditional) density for a hierarchical model for a gamma density whose rate parameter is also gamma:  $$X \mid B \sim \operatorname{Gamma}(a,B), \\ B \sim \operatorname{Gamma}(c,d), \\ f_{X\mid B}(x\mid b) = \frac{b^a x^{a-1} e^{-bx}}{\Gamma(a)}, \\ f_B(b) = \frac{d^c b^{c-1} e^{-db}}{\Gamma(c)}.$$  We note that the marginal density is
$$\begin{align}
f_X(x) &= \int_{b=0}^\infty f_{X \mid B}(x \mid b) f_B(b) \, db \\
&= \frac{d^c x^{a-1}}{\Gamma(a)\Gamma(c)} \int_{b=0}^\infty b^{a+c-1} e^{-(x+d)b} \, db \\
&= \frac{\Gamma(a+c) d^c x^{a-1}}{\Gamma(a)\Gamma(c)(x+d)^{a+c}} \int_{b=0}^\infty \frac{(x+d)^{a+c} b^{a+c-1} e^{-(x+d)b}}{\Gamma(a+c)} \, db \\
&= \frac{\Gamma(a+c) d^c x^{a-1}}{\Gamma(a)\Gamma(c) (x+d)^{a+c}}. \tag{1}
\end{align}$$
The support is on $x \in (0,\infty)$, so this is not beta distributed.  Instead, this is known as a generalized beta prime distribution or compound gamma distribution.  The latter is a special case of the former with $p = 1$.
