[Rd] qnbinom with small size is slow

[Constantin Ahlmann-Eltze]

Hi Martin,

thanks for verifying. I agree that the Cornish-Fisher seems to struggle
with the small size parameters, but I also don't have a good idea how to
replace it.

But I think fixing do_search() is possible:

I think the problem is that when searching to the left y is decremented
only if `pnbinom(y - incr, n, pr, /*l._t.*/TRUE, /*log_p*/FALSE)) < p` is
FALSE. I think the solution is to move the update of y before the if.
However, I need to make this slightly awkward check if incr == 1, so that
the return in line 123 and the do-while block at the end of qnbinom() do
not need to be modified.

diff --git a/src/nmath/qnbinom.c b/src/nmath/qnbinom.c
index b313ce56b2..16845d9373 100644
--- a/src/nmath/qnbinom.c
+++ b/src/nmath/qnbinom.c
@@ -49,10 +49,18 @@ do_search(double y, double *z, double p, double n,
double pr, double incr)
 {
     if(*z >= p) {	/* search to the left */
 	for(;;) {
+        y = fmax2(0, y - incr);
 	    if(y == 0 ||
-	       (*z = pnbinom(y - incr, n, pr, /*l._t.*/TRUE, /*log_p*/FALSE)) < p)
-		return y;
-	    y = fmax2(0, y - incr);
+	       (*z = pnbinom(y, n, pr, /*l._t.*/TRUE, /*log_p*/FALSE)) < p){
+            if(incr == 1){
+                // we know that the search is stopped if incr == 1
+                // and we know that the correct result is just right
+                // of the current y
+                return y + 1;
+            }else{
+                return y;
+            }
+        }
 	}
     }
     else {		/* search to the right */


With this patch, we get the expected result

> qnbinom(0.9999, mu = 3, size = 1e-4)
[1] 7942

I have updated the pull request at https://github.com/r-devel/r-svn/pull/11
and it is currently checking if the change breaks anything.

Best,
Constantin


Am 20.08.20 um 22:27 schrieb Martin Maechler:

Constantin Ahlmann-Eltze via R-devel
    on Mon, 10 Aug 2020 10:05:36 +0200 writes:

    > Thanks Ben for verifying the issue. It is always reassuring to hear
    > when others can reproduce the problem.

    > I wrote a small patch that fixes the issue
    > (https://github.com/r-devel/r-svn/pull/11):

    > diff --git a/src/nmath/qnbinom.c b/src/nmath/qnbinom.c
    > index b313ce56b2..d2e8d98759 100644
    > --- a/src/nmath/qnbinom.c
    > +++ b/src/nmath/qnbinom.c
    > @@ -104,6 +104,7 @@ double qnbinom(double p, double size, double prob,
    > int lower_tail, int log_p)
    > /* y := approx.value (Cornish-Fisher expansion) :  */
    > z = qnorm(p, 0., 1., /*lower_tail*/TRUE, /*log_p*/FALSE);
    > y = R_forceint(mu + sigma * (z + gamma * (z*z - 1) / 6));
    > +    y = fmax2(0.0, y);

    > z = pnbinom(y, size, prob, /*lower_tail*/TRUE, /*log_p*/FALSE);

    > I used the https://github.com/r-devel/r-svn repo and its continuous
    > integration tools to check that it doesn't break any existing tests:
    > https://github.com/r-devel/r-svn/actions/runs/201327042

    > I have also requested a Bugzilla-account, but haven't heard
anything back yet.

    > Best,
    > Constantin

Thank you for the report, and Ben for his experiment.

And, indeed in this case, this returns  0  much more  quickly.

Note that this could be even much more quickly: The
Cornish-Fisher expansion is not really of much use here, ...
and quick check would just see that pnbinom(0, size, prob) >

Note however, that in other cases, results for  small 'size'
are *still* not good  (and *not* influenced by your patch !!),
e.g.,

## Other examples, not giving 0, are fast already but  *in*accurate:
qnbinom(.9999, mu=3, size=1e-4)
## [1] 8044

## but
str(ur1 <- uniroot(function(q) pnbinom(q, mu=3, size=1e-4) - 0.9999,
c(7000,8000)))
## List of 5
##  $ root      : num 7942
##  $ f.root    : num 1.52e-09
##  $ iter      : int 18
##  $ init.it   : int NA
##  $ estim.prec: num 6.49e-05

## and this of course does not change when asking for more precision :

str(ur2 <- uniroot(function(q) pnbinom(q, mu=3, size=1e-4) - 0.9999,
c(7000,8000), tol=1e-12))
## List of 5
##  $ root      : num 7942  <<< correct is 7942, not 8044 !!!
##  $ f.root    : num 1.52e-09
##  $ iter      : int 47
##  $ init.it   : int NA
##  $ estim.prec: num 7.28e-12

----------

so, in principle the C-internal  search() function really should be
improved for such  ( somewhat extreme!! )  cases.
or ... ?? ... a different approximation should be used for such
extreme small 'size' (and  prob := size/(size+mu)  ) ...

Martin Maechler
ETH Zurich   and  R Core team

	[[alternative HTML version deleted]]