[R] inconsistency between gamma and choose functions
Thomas Lumley
tlumley at u.washington.edu
Fri Dec 14 00:05:54 CET 2001
On Fri, 14 Dec 2001, Cowpertwait, Paul wrote:
> Please can someone explain why I seem to get these contradictory results?
>
> choose(5,2)
> [1] 10
> gamma(6)/(gamma(3)*gamma(4))
> [1] 10
> gamma(6)/(gamma(3)*gamma(4)) == choose(5,2)
> [1] TRUE
> # all's well so far.
>
> # now look what happens:
> gamma(21)/(gamma(6)*gamma(16)) == choose(20,5)
> [1] FALSE
>
> # check individual terms:
> gamma(21)/(gamma(6)*gamma(16))
> [1] 15504
> choose(20,5)
> [1] 15504
> # so they are the same, BUT we get FALSE when comparing - a contradiction!
> gamma(21)/(gamma(6)*gamma(16)) == choose(20,5)
> [1] FALSE
If you look at
> choose(20,5)-gamma(21)/(gamma(6)*gamma(16))
[1] -2.182787e-11
you see what is happening.
choose(20,5) is exactly 15504 -- eg check that
> choose(20,5)-15504==0
[1] TRUE
but the ratio of gamma functions involves some slight rounding error
> gamma(21)/(gamma(6)*gamma(16))-15504
[1] 2.182787e-11
Now if you look at the precision of double precision arithmetic
> Machine()$double.eps*15504
[1] 3.442580e-12
so the error is about 7 times machine resolution. This is pretty good, as
the numerator and denominator are both quite large.
> gamma(21)
[1] 2.432902e+18
Even given full 16-digit accuracy for gamma(x) the division would cause
some error. Furthermore in computing choose() we know that the result is
an integer, and so we round it (see src/nmath/choose.c). We can't do
this for the gamma function as it can be used for non-integer
arguments.
In fact there's another source of difference because choose() actually
uses the logarithm of the gamma function and exponentiates the result.
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list