Question:
Let $X|Y = y\sim\text{Poisson}(y)$ and $Y\sim\text{Gamma}(\alpha, \lambda)$. Find join density $f_{X,Y}(x,y)$ and find the probability density function of $X$ (simplify until there are no integrals).
Answer:
So, I believe that the joint density function for this will be:
$$f_{X}(x|y) \cdot f_{Y}(y) = f_{X,Y}(x,y) = \frac{y^x \cdot e^{-y}}{x!} \cdot \frac{\lambda e^{-\lambda y}(\lambda y)^{\alpha-1}}{\Gamma(\alpha)}.$$
Next, in order to find the pdf of $X$, I need to $$\int_{-\infty}^{\infty} f_{X,Y}(x,y)\,dy.$$ I tried to use integration by parts, but I couldn't find any answer. Does anyone know the approach to integrating $f_{X,Y}(x,y)\,dy$?
Regarding the joint density of $X$ and $Y$, you must be clear about the support of the distribution: your computation is correct, but you should also write $$X \in \{0, 1, 2, \ldots\} \cap Y \in (0,\infty).$$
As for the marginal distribution of $X$, you would indeed perform the integration with respect to $Y$ but only on the support of $Y$. If we accept that the integral of a gamma density with rate $\lambda$ and shape $\alpha$ is $1$ over its support, then $$\int_{y=0}^\infty e^{-y} \frac{y^x}{x!} \frac{\lambda^\alpha y^{\alpha-1} e^{-\lambda y}}{\Gamma(\alpha)} \, dy = \frac{\lambda^\alpha}{x! \, \Gamma(\alpha)} \int_{y=0}^\infty y^{x+\alpha-1} e^{-(\lambda+1)y} \, dy.$$ Now if we recognize this integrand as being proportional to a gamma distribution with shape $\alpha^* = x + \alpha$, and rate $\lambda^* = \lambda+1$, we note we would need to have a factor of $(\lambda^*)^{\alpha^*}$ in the numerator and $\Gamma(\alpha^*)$ in the denominator; thus $$\Pr[X = x] = \frac{\lambda^\alpha \Gamma(x+\alpha)}{(\lambda+1)^{x+\alpha} x! \, \Gamma(\alpha)} \int_{y=0}^\infty \frac{(\lambda^*)^{\alpha^*} y^{\alpha^* - 1} e^{-\lambda^* y}}{\Gamma(\alpha^*)} \, dy,$$ and now because this integral is $1$, and $x$ is a nonnegative integer, we find $$\Pr[X = x] = \binom{x+\alpha-1}{x} \frac{\lambda^\alpha}{(\lambda+1)^{x+\alpha}}, \quad x = 0, 1, 2, \ldots.$$ What kind of distribution is this? It looks like it could be negative binomial. Indeed, if we let $$p = \frac{1}{\lambda+1}, \quad 1 - p = \frac{\lambda}{\lambda+1},$$ we find that $$(1-p)^{\alpha} p^x = \frac{\lambda^\alpha}{(\lambda+1)^{x+\alpha}},$$ as desired; thus the marginal distribution of $X$ is negative binomial with rate $r = \alpha$, and probability of "success" $p = (1+\lambda)^{-1}$; where $X$ counts the random number of "successes" that are observed in a series of IID Bernoulli trials until the $r^{\rm th}$ failure (if $\alpha$ is a positive integer, which it need not be).