>>>>> 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
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
