2
$\begingroup$

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?

$\endgroup$
  • 5
    $\begingroup$ 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). $\endgroup$ – user10525 Jul 14 '12 at 17:16
  • $\begingroup$ Indeed, it works, I feel stupid at the moment. Thanks a lot. $\endgroup$ – ehmicky Jul 14 '12 at 17:17
  • $\begingroup$ Thanks, I didn't know if answering my own question was considered a bad thing on this website. $\endgroup$ – ehmicky Jul 15 '12 at 19:19
  • 1
    $\begingroup$ @ehmicky: Not at all; in fact, it's encouraged. $\endgroup$ – cardinal Jul 15 '12 at 19:48
4
$\begingroup$

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
$\endgroup$
  • 3
    $\begingroup$ $X$ and $Y$ also should be independent, which is worth mentioning. $\endgroup$ – cardinal Jul 15 '12 at 19:47

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.