Finding a distribution of a random variable generated using a Monte Carlo method

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.

1. Generate $U_1,U_2,U_3 \sim U[0,1]$, i.i.d

2. If $U_1 > \frac{e}{\alpha+e}$, go to 4

3. Let $Y=U_2^{\frac{1}{\alpha}}.$ If $U_3 > e^{-Y}$ go to 1, else go to 5.

4. Let $Y=1-\log{U_2}$. If $U_3>Y^{\alpha-1}$ go to 1.

5. 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]$$

$$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]:

-

Let us compare the probabilities of the possible outputs of one run of the algorithm:

• The path 1.-2.-3.-5. produces $X=U_2^{1/\alpha}$ and has probability $\frac{\mathrm e}{\alpha+\mathrm e}$ (for the step 2.-3.) times $\mathrm e^{-X}$ (for the step 3.-5.) of happening.
• The path 1.-2.-4.-5. produces $X=1-\log U_2$ and has probability $\frac{\alpha}{\alpha+\mathrm e}$ (for the step 2.-4.) times $X^{\alpha-1}$ (for the step 4.-5.) of happening.
• All the other paths go back to 1.

Since the runs are i.i.d., this shows that after a geometrically distributed number of unsuccessful runs (whose number is irrelevant), the successful run will produce a result $X$ whose density is proportional to the function $g$ such that $$g(x)=\frac{\mathrm e}{\alpha+\mathrm e}\cdot\mathrm e^{-x}\cdot f_1(x)+\frac{\alpha}{\alpha+\mathrm e}\cdot x^{\alpha-1}\cdot f_2(x),$$ where $f_1$ is the density of $U_2^{1/\alpha}$ and $f_2$ is the density of $1-\log U_2$. Standard computations yield $$f_1(x)=\alpha x^{\alpha-1}[0\lt x\leqslant1],\qquad f_2(x)=\mathrm e^{1-x}[x\gt1].$$ Hence $g(x)=\frac{\alpha\mathrm e}{\alpha+\mathrm e} x^{\alpha-1}\mathrm e^{-x}[x\gt0]$, which is proportional to the density of a gamma random variable of parameter $\alpha$, that is, $\Gamma(\alpha)^{-1}x^{\alpha-1}\mathrm e^{-x}[x\gt0]$. Finally, the distribution of $X$ is gamma with parameter $\alpha$.

-
Thank you for a complete answer. I can see where I went wrong. I put $\mathbb{P}(U_3\leq e^{{U_2}^{\frac{1}{\alpha}}})$ (=$e\Gamma(\alpha,1)$) instead of $x^{a-1}$. This is the step I am not sure about. – Tom Artiom Fiodorov Feb 29 '12 at 21:24

This Mathematica 7 code was used for generation of the plots. I include this for completion.

b = Range[1000];
a = 0.4;
For[i = 1, i <= 1000, i++,
X = 0;
While[X == 0,
U1 = RandomReal[];
U2 = RandomReal[];
U3 = RandomReal[];
If[U1 > E/(E + a),
Y = 1 - Log[U2];
If[U3 <= Y^{a + 1}, X = Y],
Y = (U2)^(1/a);
If[U3 <= Exp[Y], X = Y]
]
]
i;
b[[i]] = X
]
ListPlot[b]