[Rd] [R] choose(n, k) as n approaches k

peter dalgaard pd@|gd @end|ng |rom gm@||@com
Tue Jan 14 16:07:44 CET 2020


Yep, that looks wrong (probably want to continue discussion over on R-devel)

I think the culprit is here (in src/nmath/choose.c)
 
   if (k < k_small_max) {
        int j;
        if(n-k < k && n >= 0 && R_IS_INT(n)) k = n-k; /* <- Symmetry */
        if (k <  0) return 0.;
        if (k == 0) return 1.;
        /* else: k >= 1 */

if n is a near-integer, then k can become non-integer and negative. In your case, 

n == 4 - 1e-7
k == 4
n - k == -1e-7 < 4
n >= 0 
R_IS_INT(n) = TRUE (relative diff < 1e-7 is allowed)

so k gets set to

n - k == -1e-7

which is less than 0, so we return 0. However, as you point out, 1 would be more reasonable and in accordance with the limit as n -> 4, e.g.

> factorial(4 - 1e-10)/factorial(1e-10)/factorial(4) -1
[1] -9.289025e-11

I guess that the fix could be as simple as replacing n by R_forceint(n) in the k = n - k step.

-pd



> On 14 Jan 2020, at 00:33 , Wright, Erik Scott <ESWRIGHT using pitt.edu> wrote:
> 
> This struck me as incorrect:
> 
>> choose(3.999999, 4)
> [1] 0.9999979
>> choose(3.9999999, 4)
> [1] 0
>> choose(4, 4)
> [1] 1
>> choose(4.0000001, 4)
> [1] 4
>> choose(4.000001, 4)
> [1] 1.000002
> 
> Should base::choose(n, k) check whether n is within machine precision of k and return 1?
> 
> Thanks,
> Erik
> 
> ***
> sessionInfo()
> R version 3.6.0 beta (2019-04-15 r76395)
> Platform: x86_64-apple-darwin15.6.0 (64-bit)
> Running under: macOS High Sierra 10.13.6
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: pd.mes using cbs.dk  Priv: PDalgd using gmail.com



More information about the R-devel mailing list