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?
a<-rgamma(m,1,1)
,b<-rgamma(m,1,1)
,c<-a/(a+b)
and then calculatequantile(c,1:9/10)
. $\endgroup$