How do I show that $Y=2\sqrt{X_1X_2}\sim$Gamma$(2p,1)$? Suppose that $X_1\sim \operatorname{Gamma}(p,1)$ and independently, $X_2\sim \operatorname{Gamma}(p+1/2,1)$. Show that $Y=2\sqrt{X_1X_2}\sim\operatorname{Gamma}(2p,1)$.
This problem followed a section on bivariate transformations, so I made attempts to solve this problem using that method. However, I was getting an integral that does not properly integrate, as was my professor.
I also attempted to use CDFs and mgfs, but I was not able to make any progress. How would you solve this problem?
 A: If $X$ is a Gamma $\mathcal G(p,1)$ variate, its density is
$$f_X(x)= \frac{x^{p-1}}{(p-1)!}e^{-x}\mathbb I_{\mathbb R^*_+}(x)$$
Therefore, the probability density of $Y_1=\sqrt{X_1}$ is
$$f_{Y_1}(y)=\frac{y^{2(p-1)}}{(p-1)!}e^{-y^2}\mathbb I_{\mathbb R^*_+}(y)\times\overbrace{\left|\frac{\text{d}x}{\text{d}y}\right|}^{\text{Jacobian}}
=2y\times\frac{y^{2(p-1)}}{(p-1)!}e^{-y^2}\mathbb I_{\mathbb R^*_+}(y)
=\frac{2y^{2p-1}}{(p-1)!}e^{-y^2}\mathbb I_{\mathbb R^*_+}(y)$$
And the probability density of $Y_2=\sqrt{X_2}$ is
$$f_{Y_2}(y)=\frac{2y^{\overbrace{2p+1-1}^{2p}}}{\Gamma(p+1/2)}e^{-y^2}\mathbb I_{\mathbb R^*_+}(y)$$
Hence the joint density of $(Y_1,Y_2)$ is
$$g(y_1,y_2)=\frac{4y_1^{2p-1}y_2^{2p}}{(p-1)!\Gamma(p+1/2)}e^{-y_1^2-y_2^2}\mathbb I_{\mathbb R^*_+}(y_1)\mathbb I_{\mathbb R^*_+}(y_1)$$
Considering the change of variables from $(y_1,y_2)$ to $(z=y_1y_2,y_2)$, the joint density of $(Z,Y_2)$ is
$$h(z,y_2)=g(z/y_2,y_2)\overbrace{\left|\frac{\text{d}(y_1,y_2)}{\text{d}(z,y_2)}\right|}^{\text{Jacobian}}=g(z/y_2,y_2)\left|\frac{\text{d}y_1}{\text{d}z}\right|=g(z/y_2,y_2)y_2^{-1}$$
and the density of $Z$ is the marginal
\begin{align*}
f_Z(z) &= \int_0^\infty g(z/y_2,y_2)y_2^{-1}\,\text{d}y_2\\
       &=\int_0^\infty \frac{4z^{2p-1}y_2^{2p-(2p-1)}}{(p-1)!\Gamma(p+1/2)}e^{-z^2y_2^{-2}-y_2^2}y_2^{-1}\,\text{d}y_2\\
       &= \frac{4z^{2p-1}}{(p-1)!\Gamma(p+1/2)}\int_0^\infty e^{-z^2y_2^{-2}-y_2^2}\,\text{d}y_2\\
       &= \frac{4z^{2p-1}}{(p-1)!\Gamma(p+1/2)} \frac{\sqrt{\pi}}{2}e^{-2z}
\end{align*}
[the last integral is formula 3.325 in Gradshteyn & Ryzhik, 2007]
Hence the density of $S=2Z$ is
$$f_S(s)=\frac{\sqrt{\pi}2^{1-2p}s^{2p-1}}{(p-1)!\Gamma(p+1/2)}e^{-s}=\frac{s^{2p-1}}{\Gamma(2p)}e^{-s}$$
[where the constant simplifies by the Legendre duplication formula]
A: I think the moment generating function approach works fine, but easier if we consider MGF of $\ln Y$.
Assuming of course $\mathsf{Gamma}(p,1)$ refers to shape $p$ parameterization, i.e. with density $$f(x)=\frac{e^{-x}x^{p-1}}{\Gamma(p)}\mathbf1_{x>0}$$ with $p>0$ as in @Xi'an's answer.
We have
\begin{align}
E\left[e^{t\ln Y}\right]&=E\left[Y^t\right]
\\&=2^tE\left[X_1^{t/2}\right]E\left[X_2^{t/2}\right]
\end{align}
For $t>-2p$ where $p>0$, clearly
$$E\left[X_1^{t/2}\right]=\frac{\Gamma\left(p+\frac t2\right)}{\Gamma(p)}$$
And $$E\left[X_2^{t/2}\right]=\frac{\Gamma\left(p+\frac t2+\frac12\right)}{\Gamma\left(p+\frac12\right)}$$
Using Legendre's duplication formula,
\begin{align}
E\left[e^{t\ln Y}\right]&=2^t\cdot \frac{\Gamma(2p+t)\sqrt\pi/2^{2p+t-1}}{\Gamma(2p)\sqrt\pi/2^{2p-1}}
\\&=\frac{\Gamma(2p+t)}{\Gamma(2p)}
\end{align}
This is the MGF of the logarithm of a $\mathsf{Gamma}(2p,1)$ distribution evaluated at $t$ (or simply the $t$th order raw moment about $0$) where $t\in (-2p,\infty)$. As the MGF exists in an open interval containing $0$, we can conclude that $Y\sim \mathsf{Gamma}(2p,1)$ by uniqueness theorem of MGF.  In effect we see that this Gamma distribution is uniquely determined by its moments.
A: ... Comment continued: R code for the simple simulation is as shown
below. Unfortunately, the simulation gives no clue how to work the
problem. (See @Xi'an's Answer.)
set.seed(1023)
p = 2;  m = 10^6
x1 = rgamma(m,p,1);  x2 = rgamma(m,p+.5,1)
y = 2*sqrt(x1*x2)
hist(y, br=60, prob=T, col="skyblue2")
 curve(dgamma(x,2*p,1), add=T, col="red", n=10001)


