# If $X \sim \mathrm{Gamma}(\alpha,1)$ and $Y \sim \mathrm{Gamma}(\beta,1)$, is $X/(X+Y) \sim \mathrm{Beta}(\alpha,\beta)$?

According to Wikipedia :

$\mathrm{Beta}(\alpha,\beta) = \mathrm{Gamma}(\alpha,\theta) / (\mathrm{Gamma}(\alpha,\theta) + \mathrm{Gamma}(\beta,\theta))$

However, when I try to simulate this in R :

> m <- 2^11
> a <- rgamma(m,1,1) / (rgamma(m,1,1) + rgamma(m,1,1))
> quantile(a,1:9/10)
0.05576023 0.12110211 0.19886341 0.29088744 0.41851181
0.56698135 0.80151216 1.20748646 2.13406961
> qbeta(1:9/10,1,1)
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9


Those two distributions don't look the same on the lower and upper percentiles to me, even when I increase the simulation loop index m. Simulation works for higher $\alpha$ and $\beta$ though.

I understand than since $\mathrm{Gamma}(1,1)$ is identical to $\mathrm{Exponential}(1)$, dividing one exponential number by the sum of two exponential numbers can produce a result greater than 1, whereas the Beta distribution supports only real numbers on $(0,1)$, so I don't understand how the formula above could be correct for any $\alpha$ and any $\beta$ (even small ones) ?

Did I make a confusion between the Gamma function and the Gamma distribution, or is it something else?

• The problem is that you are generating different numbers in the numerator and denonimator for $Gamma(\alpha,\theta)$, try using the same values and it will work. For example, a<-rgamma(m,1,1), b<-rgamma(m,1,1), c<-a/(a+b) and then calculate quantile(c,1:9/10).
– user10525
Jul 14 '12 at 17:16
• Indeed, it works, I feel stupid at the moment. Thanks a lot. Jul 14 '12 at 17:17
• Thanks, I didn't know if answering my own question was considered a bad thing on this website. Jul 15 '12 at 19:19
• @ehmicky: Not at all; in fact, it's encouraged. Jul 15 '12 at 19:48

The correct formula is: If $X \sim \mathrm{Gamma}(\alpha,\theta)$ and $Y \sim \mathrm{Gamma}(\beta,\theta)$, and $X$ and $Y$ are independent, then:

$\mathrm{Beta}(\alpha,\beta) = X/(X+Y)$

Which requires the variable $X$ being the same in the numerator and denominator.
A correct simulation in R should be:

> m <- 2^16
> x <- rgamma(m,1,1)
> y <- rgamma(m,1,1)
> a <- x/(x+y)
> quantile(a,1:9/10)
0.09974079 0.20052892 0.30086288 0.40084320 0.50010657
0.60012512 0.70155472 0.80055784 0.90050201
> qbeta(1:9/10,1,1)
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

• $X$ and $Y$ also should be independent, which is worth mentioning. Jul 15 '12 at 19:47