Marginal Density Function, Gamma and Beta distributions 
If $Y\sim\operatorname{Gamma}(\gamma,\delta)$ and $Z\sim\operatorname{Beta}(\alpha,\beta)$ then their density functions are, respectively,
  $$
f_Y(y)=\frac{\delta^\gamma}{\Gamma(\gamma)}y^{\gamma-1}e^{-\delta y},\quad y>0,\quad\gamma>0,\quad\delta>0
$$
  and
  $$
f_Z(z)=\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}z^{\alpha-1}(1-z)^{\beta-1},\quad 0\leq z\leq 1,\quad\alpha>0,\quad\beta>0.
$$
  Consider $X_1$ and $X_2$ having $\operatorname{Gamma}(a+b,1)$ and $\operatorname{Beta}(a,b)$ distributions, respectively, where $a,b>0$.  Assume that $X_1$ and $X_2$ are independent.

How do i find the marginal density functions of $Y_1 = X_1X_2$ and $Y_2 = X_1(1-X_2)$? 
I know that the marginal density function can be derived from the joint density, but since the joint is not given, how do I create it? 
Also how do I manipulate the gamma function? first time I have come across it.
 A: Plugging in the definition, $X_1$ following $\operatorname{Gamma}(a+b,1)$ means its density is
$$f_{X_1}(x_1) = \frac1{ \Gamma(a+b)}\, x_1^{a+b-1} e^{-x_1} \qquad \text{for}~~ 0 < x_1 < \infty $$
The density of $X_2$ is
$$f_{X_2}(x_2) = \frac{ \Gamma(a+b) }{ \Gamma(a) \Gamma(b) }\, x_2^{a-1} (1 - x_2)^{b-1} \qquad \text{for}~~ 0<x_2<1$$
The fact that $X_1 \perp X_2$ means their joint density is just the direct product
$$f_{X_1X_2}(x_1,\, x_2) = \frac{ x_1^{a+b-1} e^{-x_1} \, x_2^{a-1} (1 - x_2)^{b-1} }{ \Gamma(a) \Gamma(b) } \qquad \text{for}~~ \begin{cases}
0<x_1<\infty \\
0<x_2<1 \end{cases} \tag{1}\label{joint density}$$
The 2-dim transformation is
$$\begin{cases}
Y_1 = X_1 X_2 \\[1.5ex]
Y_2  = X_1 (1 - X_2)
\end{cases}  \Longleftrightarrow \begin{cases}
X_1 = Y_1 + Y_2 \\[2ex]
X_2 = \dfrac{ Y_1 }{ Y_1 + Y_2}
\end{cases} \qquad \text{where}~~ \begin{cases}
0<y_1<\infty \\
0<y_2<\infty \end{cases}$$
with the Jacobian (of the inverse mapping) as
$$J = \left| \begin{matrix} \dfrac{ \partial x_1}{ \partial y_1} & \dfrac{ \partial x_1}{ \partial y_2} \\
\dfrac{ \partial x_2}{ \partial y_1} & \dfrac{ \partial x_2}{ \partial y_2}\end{matrix} \right| = \left| \begin{matrix} 1 & 1 \\
\dfrac{ y_2 }{ (y_1 +y_2)^2 } & \dfrac{ -y_1 }{ (y_1 +y_2)^2 } \end{matrix} \right| = \frac{-1}{ y_1 + y_2 }$$
The transformed joint density for $Y_1$ and $Y_2$ is
\begin{align}
f_{Y_1Y_2}( y_1 ,~y_2 ) &= |J| \cdot f_{X_1X_2}( x_1,\, x_2)\Bigg|_{\substack{x_1 = y_1+y_2 \\ x_2 = \frac{y_1}{y_1 + y_2}}} \qquad \text{, plug in Eq.(\ref{joint density})}\\
&= \frac1{ y_1 + y_2} \cdot \frac{ (y_1 + y_2)^{a+b-1} e^{-(y_1 + y_2)}  }{ \Gamma(a) \Gamma(b) }\, \left(\frac{y_1}{ y_1 + y_2}\right)^{a-1} \left(\frac{y_2}{ y_1 + y_2} \right)^{b-1} \\
&= \frac1{ \Gamma(a) \Gamma(b) }\, y_1^{a-1} y_2^{b-1} e^{-(y_1 + y_2)}  \qquad \text{for}~~0<y_1<\infty ,~0<y_2<\infty
\end{align}
The marginal density of $Y_1$ can be obtained from the joint as
\begin{align}
f_{Y_1}(y_1) &= \int_{y_2 = 0}^{\infty} f_{Y_1Y_2}( y_1 ,~y_2 ) \,\mathrm{d}y_2  \\
&= \frac1{\Gamma(a)} y_1^{a-1} e^{-y_1} \int_{y_2 = 0}^{\infty} \frac1{\Gamma(b)} y_2^{b-1} e^{-y_2}  \,\mathrm{d}y_2  \qquad \small\begin{aligned}[c]
&\text{integral is just the} \\ 
&\text{kernel of Gamma distribution}\end{aligned} \\
&= \frac1{\Gamma(a)} y_1^{a-1} e^{-y_1} 
\end{align}
Thus one identifies the distribution of $Y_1$ as $\operatorname{Gamma}(a,1)$.
Similarly, or noting the symmetry in the joint $f_{Y_1Y_2}( y_1 ,~y_2 )$, we have $Y_2$ follows $\operatorname{Gamma}(b,1)$.
