4
$\begingroup$

I am required to show that the following accept/reject algoritmh produces observations from a beta distribution with parameters $\alpha$ and $\beta$.


  1. Generate $U_1$ and $U_2$ iid uniform(0,1) random variables. Set $V_1=U^{1/{\alpha}}$ and $V_2=U_2 ^{1/{\beta}}$
  2. Set $W=V_1+V_2$. If $W\leq 1 $, set $X=V_1 /W $; else go to step 1.
  3. Deliver $X$

Thus, I need to find the conditional distribution of step 2. Now we have:

$$P \left[X \leq x \right]= P \left[ V_1 /W \leq x | W \leq 1 \right]= \frac{P \left[V_1 /W \leq x \cap W \leq 1 \right]}{P \left[ W \leq 1 \right]} $$

A simple transformation shows that $V_1 \sim B \left( \alpha ,1 \right) $ and $V_2 \sim B \left(1, \beta \right)$ but other than that I do not know how to proceed as I cannot simplify the fraction. Had I known the distribution of $W=V_1+V_2$ things would have been easier but this is not the case here.

Can you please help me compute these probabilities? Thank you.

$\endgroup$
4
  • $\begingroup$ I think you have made a mistake in the algorithm. Consider if $\alpha = \beta = 1$. Then $V_1 \sim Uni(0,1)$ and $V_2 \sim Uni(0,1)$. Then we should get $Beta(1,1)$ (which is equal to $Uni(0,1)$) at the end, but it's easy to see using samples that $\frac{V_1}{V_1+V_2}$ is not Uniform. I believe it should be: $V_1 = (\ln{U_1})^{\frac{1}{\alpha}}$, and similar for $V_2$. $\endgroup$ – Cam.Davidson.Pilon Dec 31 '13 at 23:42
  • $\begingroup$ @Cam.Davidson.Pilon You can find the algorithm in page 272 of Introduction to Mathematical Statistics, 7th edition by Hogg, Craig McKean. Unless they are mistaken, I am not either. $\endgroup$ – JohnK Dec 31 '13 at 23:45
  • $\begingroup$ @Cam.Davidson.Pilon you have neglected to make use of the fact that they are rejecting many of the $V_1 / (V_1 + V_2)$. The following R-code shows that this is correct in your case: V <- rbeta(10000, 1, 1); Z <- rbeta(10000, 1, 1); Y <- ifelse(V + Z <= 1, V / (V + Z), NA); plot.density(W, na.rm = TRUE); $\endgroup$ – guy Dec 31 '13 at 23:47
  • $\begingroup$ Ah, yes, my mistake -- I did forget to perform the rejection. $\endgroup$ – Cam.Davidson.Pilon Dec 31 '13 at 23:56
5
$\begingroup$

First, it should be a Beta($\alpha, 1$) and Beta($\beta, 1$). Look at the numerator.

$$ P(V_1 / W \le x \cap W \le 1) = \int_0^1 \int_0 ^ {wx} f_{V_1}(v) f_{V_2}(w - v) \ dv \ dw, $$ using a substitution to get $f_{V_1, W}$ from $f_{V_1, V_2}$. Next, interchange order of integration

$$ \int_0 ^ x \int_{v/x} ^ 1 \alpha v^{\alpha - 1} \beta (w - v)^{\beta - 1} \ dw \ dv = \int_0 ^ x \alpha v^{\alpha - 1}\left[(1 - v)^\beta - v^\beta \left(\frac{1 - x}{x}\right)^\beta\right]. $$ If $g(x)$ is the density of our proposed rejection sampler, then $g(x)$ is equal to the derivative of this expression, up-to a factor of $P(W \le 1)$,

$$ g(x) \propto \frac{d}{dx} \int_0 ^ x \alpha v^{\alpha - 1}\left[(1 - v)^\beta - v^\beta \left(\frac{1 - x}{x}\right)^\beta\right]. $$

By Leibniz rule,

$$ g(x) \propto \alpha \underbrace{x^{\alpha - 1} \left[(1 - x)^\beta - x^\beta \left(\frac{1 - x}{x}\right)^\beta\right]}_{= 0} + \int_0 ^ x \frac{d}{dx} \alpha v^{\alpha - 1}\left[(1 - v)^\beta - v^\beta \left(\frac{1 - x}{x}\right)^\beta\right] \ dv \\ \propto \int_0 ^ x v^{\alpha + \beta - 1} \frac{d}{dx} \left(\frac{1 - x}{x}\right)^\beta \ dv. $$ This final integral can be evaluated to get

$$ g(x) \propto x^{\alpha - 1}(1 - x)^{\beta - 1}. $$

This argument is valid for $x \in [0, 1]$ while the density is trivially $0$ otherwise. This is the kernel of a Beta($\alpha, \beta$) distribution, so the conclusion follows.

$\endgroup$
3
  • $\begingroup$ What about the denominator though? Can we merely overlook it? $\endgroup$ – JohnK Dec 31 '13 at 23:51
  • 1
    $\begingroup$ @JohnK you don't need it. We have shown the density is equal $x^{\alpha - 1}(1-x)^{\beta-1}$ times some constant. Because we know that $g(x)$ integrates to 1 we can recover the constant $\Gamma(\alpha+\beta)/[\Gamma(\alpha)\Gamma(\beta)]$. You can actually redo the whole thing keeping track of all the constants if you want, it might be a good exercise - you can get the denominator in closed form easily enough, i.e. That isn't the hard part. $\endgroup$ – guy Dec 31 '13 at 23:59
  • 1
    $\begingroup$ I threw away three terms in calculating $g(x) \propto x^{\alpha-1}(1 - x)^{\beta-1}$; $\frac{\alpha\beta}{\alpha + \beta}$ if you follow along the argument and understand the steps. These terms are associated with the numerator. We can calculate in closed form $P(W \le 1) = \frac{\alpha \Gamma(\alpha)\Gamma(\beta + 1)}{\Gamma(\alpha + \beta + 1)}$ directly from $P(W \le 1) = \int_0 ^ 1 \int_v ^ 1 \alpha\beta v^{\alpha - 1}(w - v)^{\beta - 1}$. $P(W \le 1)$ cancels exactly with $\alpha \beta / (\alpha + \beta)$ to give the normalizing constant. $\endgroup$ – guy Jan 1 '14 at 0:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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