0
$\begingroup$

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

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$: scatter plot of X

And of Gamma[0.4,1]:

scatter plot of Gamma[a,1]

2 Answers 2

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

  • 0
    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.2012-02-29
0

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] ListPlot[RandomReal[GammaDistribution[a, 1], 1000]]