# Accept/reject for beta distribution

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.

• 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$. – Cam.Davidson.Pilon Dec 31 '13 at 23:42
• @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. – JohnK Dec 31 '13 at 23:45
• @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); – guy Dec 31 '13 at 23:47
• Ah, yes, my mistake -- I did forget to perform the rejection. – Cam.Davidson.Pilon Dec 31 '13 at 23:56

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.

• What about the denominator though? Can we merely overlook it? – JohnK Dec 31 '13 at 23:51
• @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. – guy Dec 31 '13 at 23:59
• 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. – guy Jan 1 '14 at 0:16