Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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]

share|improve this question

2 Answers 2

up vote 1 down vote accepted

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

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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.