From: Martin Maechler <maechler_at_stat.math.ethz.ch>

Date: Thu, 17 Dec 2009 15:10:49 +0100

* >>> lgamma(x) and lfactorial(x) are defined to return
*

* >>>
*

>>> ln|Gamma(x)| {= log(abs(gamma(x)))} or ln|Gamma(x+1)| respectively.

* >>>
*

* >>> Unfortunately, we haven't chosen the analogous definition for
*

* >>> lchoose().
*

* >>>
*

* >>> So, currently
*

* >>>
*

* >>> > lchoose(1/2, 1:10)
*

* >>> [1] -0.6931472 -2.0794415 NaN -3.2425924 NaN -3.8869494
*

* >>> [7] NaN -4.3357508 NaN -4.6805913
*

* >>> Warning message:
*

* >>> In lchoose(n, k) : NaNs produced
*

* >>> >
*

* >>>
*

* >>> which (the NaN's) is not particularly useful.
*

* >>> (I have use case where I really have to workaround those NaNs.)
*

* >>>
*

* >>> I herebey propose to *amend* the definition of lchoose() such
*

* >>> that it behaves analogously to lgamma() and lfactorial(),
*

* >>> i.e., to return
*

* >>>
*

* >>> log(abs(choose(.,.))
*

* >>>
*

* >>> Your comments are welcome.
*

* >>> Martin Maechler, ETH Zurich
*

R-devel_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel Received on Thu 17 Dec 2009 - 14:14:01 GMT

Date: Thu, 17 Dec 2009 15:10:49 +0100

>>>>> "MM" == Martin Maechler <maechler_at_stat.math.ethz.ch> >>>>> on Tue, 15 Dec 2009 14:54:18 +0100 writes:

>>>>> "PS" == Petr Savicky <savicky_at_cs.cas.cz> >>>>> on Tue, 15 Dec 2009 10:52:43 +0100 writes:

PS> On Tue, Dec 15, 2009 at 09:49:28AM +0100, Martin Maechler wrote:

>>> ln|Gamma(x)| {= log(abs(gamma(x)))} or ln|Gamma(x+1)| respectively.

PS> I do not have a strong opinion, whether lchoose() should be log(choose(.,.)) PS> or log(abs(choose(.,.))), although i would slightly prefere the latter (with abs()). PS> One of the reasons is that this would simplify the code of the function PS> lchoose() in src/nmath/choose.c. It produces the NA by explicit testing PS> the sign, so this test could be simply removed.

PS> Let me also point out that the current implementation seems to contain PS> a bug, which causes that it is neither log(choose(.,.)) nor log(abs(choose(.,.))).

PS> x <- choose(1/2, 1:10) PS> y <- lchoose(1/2, 1:10) PS> cbind(choose=x, lchoose=y, log.choose=log(abs(x))) PS> choose lchoose log.choose PS> # [1,] 0.500000000 -0.6931472 -0.6931472 PS> # [2,] -0.125000000 -2.0794415 -2.0794415 PS> # [3,] 0.062500000 NaN -2.7725887 PS> # [4,] -0.039062500 -3.2425924 -3.2425924 PS> # [5,] 0.027343750 NaN -3.5992673 PS> # [6,] -0.020507812 -3.8869494 -3.8869494 PS> # [7,] 0.016113281 NaN -4.1281114 PS> # [8,] -0.013092041 -4.3357508 -4.3357508 PS> # [9,] 0.010910034 NaN -4.5180723 PS> # [10,] -0.009273529 -4.6805913 -4.6805913

PS> So, besides x[1], we get NA for those cases, when choose() is positive. PS> The reason seems to be the test at line 89 of src/nmath/choose.c, which is

PS> if (fmod(floor(n-k+1), 2.) == 0) /* choose() < 0 */ PS> return ML_NAN;

PS> but it should be negated. I tested it using

PS> if (fmod(floor(n-k+1), 2.) != 0) /* choose() < 0 */ PS> return ML_NAN;

MM> Indeed !!! MM> Pretty "interesting" that this thinko has not been detected for MM> quite a few years.... MM> This, of course, is an even more compelling reason to implement MM> the change of return log(abs(choose(.,.)),MM> and at the moment, I'd even plan to "backport" that to R "2.10.1 MM> patched", as the current behavior is clearly bugous.

This has now happened;

svn revisions 50775 and 50776 {R-devel and R-patched}.

Martin

R-devel_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel Received on Thu 17 Dec 2009 - 14:14:01 GMT

*
This archive was generated by hypermail 2.2.0
: Fri 18 Dec 2009 - 02:11:10 GMT
*