R help - Dec 2001 - inconsistency between gamma and choose functions

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

# the problem seems to have root in the gamma function, because:
choose(20,5)
[1] 15504
choose(20,5) == 15504
[1] TRUE
# but,
gamma(21)/(gamma(6)*gamma(16)) == 15504
[1] FALSE
# and yet ..
gamma(21)/(gamma(6)*gamma(16))
[1] 15504


# a function to 'compare' shows FALSE starts to appear when n >= 10.
Why the
inconsistency?

compare <- function (n, k) choose(n,k) =gamma(n+1)/(gamma(k+1)*gamma(n-k+1))
compare(5,2)
[1] TRUE
compare(10,2)
[1] FALSE
compare(9,2)
[1] TRUE
compare(9,3)
[1] TRUE
#etc. 

_____________________________________

Paul S.P. Cowpertwait
IIMS, Massey University, Albany, 
Private Bag 102 904 
North Shore Mail Centre
Auckland, NZ
Tel (+64) (9) 443 9799 ext 9488

http://www.massey.ac.nz/~pscowper
_____________________________________
 
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

Hi,

the problem is much more general than R and gamma().  There are two ways to
store numbers in computers -- as integers and as reals.  The integer
calculations are exact (so long you do not receive overflow), calculations
with real numbers generally isn't (gamma() works with real numbers).  The
precision is always finite and many of the nice mathematical formulas fail
to work.  You should never expect that sin(x)^2 + cos(x)^2 == 1 or
sqrt(x)^2 == x although it sometimes happens.  Rather than writing

if( real1 == real2) {

one should always consider

if( abs( real1 - real2) < very.small.number) {

This is a fundamental problem with digital computing but at least I do not
see any other efficent ways how the computers could work (I do not know much
about quantum computing, however).


Regards,

Ott Toomet




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
>
> # the problem seems to have root in the gamma function, because:
> choose(20,5)
> [1] 15504
> choose(20,5) == 15504
> [1] TRUE
> # but,
> gamma(21)/(gamma(6)*gamma(16)) == 15504
> [1] FALSE
> # and yet ..
> gamma(21)/(gamma(6)*gamma(16))
> [1] 15504
>
>
> # a function to 'compare' shows FALSE starts to appear when n >=
10. Why the
> inconsistency?
>
> compare <- function (n, k) choose(n,k) =>
gamma(n+1)/(gamma(k+1)*gamma(n-k+1))
> compare(5,2)
> [1] TRUE
> compare(10,2)
> [1] FALSE
> compare(9,2)
> [1] TRUE
> compare(9,3)
> [1] TRUE
> #etc.
>
> _____________________________________
>
> Paul S.P. Cowpertwait
> IIMS, Massey University, Albany,
> Private Bag 102 904
> North Shore Mail Centre
> Auckland, NZ
> Tel (+64) (9) 443 9799 ext 9488
>
> http://www.massey.ac.nz/~pscowper
> _____________________________________
>
>
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
> 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
>
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
>
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._