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

Question

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?

Answer 1

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