[R] inconsistency between gamma and choose functions

Ott Toomet siim at obs.ee
Fri Dec 14 02:37:16 CET 2001


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



More information about the R-help mailing list