Nested Gamma Distributions

Question

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.

Answer 1

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

Answer 2

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