I would greatly appreciate if somebody could confirm or negate my result to the following problem. I am especially not sure about "putting it all together" step.
Generate $U_1,U_2,U_3 \sim U[0,1]$, i.i.d
If $U_1 > \frac{e}{\alpha+e}$, go to 4
Let $Y=U_2^{\frac{1}{\alpha}}.$ If $U_3 > e^{-Y}$ go to 1, else go to 5.
Let $Y=1-\log{U_2}$. If $U_3>Y^{\alpha-1}$ go to 1.
Return X=Y
Find the distribution of $X$.
And this is my solution:
Route 1.
$\mathbb{P}(U_1>\frac{e}{\alpha+e})=\frac{a}{\alpha+e}$
$Y=1-\log{U_2} \text{ (change of variables)} \implies f_Y(y)=e^{1-y} \mathbf{1}_{[1,\infty)}(y)$
$\mathbb{P}({U_3\leq Y^{\alpha-1}})=\mathbb{P}({U_3\leq (1-\log{U_2})^{\alpha-1}})=\mathbb{E}[(1-\log{U_2})^{\alpha-1}]=\mathbb{E}[Y^{\alpha-1}]=e\Gamma(\alpha,1)$
Route 2.
$\mathbb{P}(U_1\leq\frac{e}{\alpha+e})=\frac{e}{\alpha+e}$
$Y=U_2^{\frac{1}{\alpha}} \text{ (change of variables)} \implies f_Y(y)=\alpha y^{\alpha-1}\mathbf{1}_{[0,1]}(y)$
$\mathbb{P}(U_3\leq e^{{U_2}^{\frac{1}{\alpha}}}) = \mathbb{E}[e^{{U_2}^{\frac{1}{\alpha}}}]$=$\displaystyle \alpha \int_0^1 e^{-y}y^{\alpha-1}\;dy=\alpha \gamma(\alpha,1)$
Where $\Gamma(\alpha,1):=\displaystyle \int_1^{\infty}y^{\alpha-1}e^{-y}\;dy$; $\gamma(\alpha,1):=\displaystyle \int_0^1 y^{\alpha-1}e^{-y}\;dy$
Putting it all together:
$f_X(y)=\frac{\alpha}{\alpha+e}\Big[ e\Gamma(\alpha,1)e^{1-y} \mathbf{1}_{[1,\infty)}(y) +(1-e\Gamma(\alpha,1))f_X(y)\Big]+$ $+\frac{e}{\alpha+e}\Big[ \alpha \gamma(\alpha,1) \alpha y^{\alpha-1}\mathbf{1}_{[0,1]}(y) +(1-\alpha \gamma(\alpha,1)f_X(y)\Big]$
And the final answer:
$f_X(y)=\frac{\Gamma(\alpha,1)e^{1-y} \mathbf{1}_{[1,\infty)}(y)+ \gamma(\alpha,1) \alpha y^{\alpha-1}\mathbf{1}_{[0,1]}(y)}{\Gamma(\alpha)}$
Some analysis of the answer: Good thing - it integrates to 1. Bad thing - I was expecting Gamma distribution.
Here's a scatter plot of X (from the algorithm), $\alpha=0.4$:
And of Gamma[0.4,1]: