From Exponential Distributions to Weibull Distribution (CDF) Problem: 
Let $X_1$ and $X_2$ be two independent exponential random variables with the PDFs $f_{X_i}(x_i)={1\over \lambda_i} \exp(-\frac{x_i}{\lambda_i})$ (where $i=1,2$). 
Also, let $Y=\frac{(X_1)^2 X_2}{a}$.
I want to find $(Y\leq x)$ i.e. $F_Y(x)=\frac{(X_1)^2 X_2}{a} \leq x$.
My attempted sol (1): 
$$\eqalign{&=(X_1)^2  \leq \frac{a  x} {X_2}\\
&=\int_0^\infty     X_1  \leq \sqrt{\frac{a  x} {z_2}}    \quad f_{X_2}(z_2) dz_2\\
&= {1\over \lambda_2} \int_0^\infty     \left(1-\exp\big(-{\sqrt\frac{a  x} {z_2 \lambda_1^2}}\big)\right) \exp(-\frac{z_2}{\lambda_2})   \quad dz_2\\ 
&=1-{1\over \lambda_2} \int_0^\infty \exp\left(-{\sqrt\frac{c} {z_2}}-\frac{z_2}{\lambda_2}\right) dz_2\tag{1}}$$
I know that $\int_0^\infty \exp\left(-{\frac{\beta} {4z_2}}-{z_2 \gamma}\right) dz_2 = \sqrt{β\over\gamma}K_1(\sqrt{\beta\gamma})$ from Table of Integrals, Series and Products, 7th edition - equation §3.324.1]. However, the final form of above equation contains $\sqrt{}$ and therefore cannot be solved by using §3.324.1.
So if you guys can comment or provide any kind of help that would be very helpful.
My attempted sol (2): 
$$\eqalign{&=(X_1)^2  \leq \frac{a  x} {X_2}\\
&=\int_0^\infty     X_2  \leq {\frac{a  x} {z_1^2}}    \quad f_{X_1}(z_1) dz_1\\
&= {1\over \lambda_1} \int_0^\infty     \left(1-\exp\big(-{\frac{a  x} {z_1^2 \lambda_2}}\big)\right) \exp(-\frac{z_1}{\lambda_1})   \quad dz_1\\ 
&=1-{1\over \lambda_1} \int_0^\infty \exp\left(-{\frac{c} {z_1^2}}-\frac{z_1}{\lambda_1}\right) dz_1\tag{1}}$$
Once again to the best my knowledge this above equation doesn't submit to any closed form solution. So I am stuck here....
Since, X is exponential r.v with mean $\lambda$, then $X^{1\over\gamma}$ is a Weibull (γ, β) random variable. Can we solve it this way? by using the CDF or pdf of weibull during conditioning?
Any kind of help will be very much appreciated.
 A: We can forget about $\lambda_1$ and $\lambda_2$ because $X_1/\lambda_1$ and $X_2/\lambda_2$ are always exponentially distributed with unit mean and the variable $Y=X_1^2X_2$ is related to $(X_1/\lambda_1)^2(X_2/\lambda_2)$ by a constant factor $\lambda_1^2\lambda_2$. So W. L. O. G, let's consider $\lambda_1=\lambda_2=1$. The probability density of $Y$ is given by
$$p(y)=\int_0^\infty dx_1\int_0^\infty dx_2\,e^{-x_1-x_2}\delta(x_1^2x_2-y)
=\int_0^\infty \frac{dx_1}{x_1^2}\,\exp\left(-x_1-\frac{y}{x_1^2}\right),$$
where $y>0$. Therefore, the complementary CDF can be found by
$$P(Y>y)=\int_y^\infty p(y')dy'=\int_0^\infty dx_1\exp\left(-x_1-\frac{y}{x_1^2}\right).$$
Mathematica says the integral gives the Meijer G-function. See http://reference.wolfram.com/language/ref/MeijerG.html. So no further simplification possible, I think. Saddle point approximation exists for large $y$.
A: As another answer has suggested the probability density function you are working with doesn't seem to have a particularly nice closed form, but in case it is of some use to you we can also calculate the characteristic function of $Y = X_1^2 X_2$, where $\lambda_1 = \lambda_2 = 1$ by
\begin{align*}
\varphi_Y(s) &= \int_0^{\infty} e^{isy} f(y) dy \\
&= \int_0^{\infty} e^{isy} \int_0^{\infty}\frac{1}{2w^{3/2}}e^{-\sqrt{w}}e^{-y/w}dwdy \\
&=\int_0^{\infty}\frac{1}{2w^{3/2}}e^{-\sqrt{w}}\int_0^{\infty} e^{(is - 1/w)y}dydw \\
&= \int_0^{\infty}\frac{1}{2w^{3/2}}e^{-\sqrt{w}}\frac{w}{1-isw}dw\\
&= \int_0^{\infty}\frac{1}{1 - isu^2}e^{-u}du\\
&=\int_0^{\infty} \frac{e^{-u}}{s^2 u^4 + 1}du +i\int_0^{\infty}\frac{su^2e^{-u}}{s^2u^4 + 1}du.
\end{align*}
Which summing over the roots of $s^2u^4 + 1 = 0$ is given by

\begin{align*} \varphi_Y(s) &= \frac{1}{4}\left(\sum_{r \,: \,s^2r^4 +
1 = 0} \frac{1}{s^2}\frac{e^{-r}
\mbox{Ei}(u)}{r^3}\bigg|_{u=0}^{u=\infty} +
\frac{i}{s}\frac{e^{-r}\mbox{Ei}(u)}{r}\bigg|_{u=0}^{u=\infty}\right).
\end{align*}

So the general case can then be recovered from $\varphi_Y(\lambda_1^2 \lambda_2/a)$, and you could then approximate the CDF with Edgeworth/Gram-Charlier series or Saddlepoint approximations.
