[Rd] qbinom (PR#13711)

[maechler at stat.math.ethz.ch]

>>>>> "PS" == Petr Savicky <savicky at cs.cas.cz>
>>>>>     on Sat, 23 May 2009 18:22:26 +0200 writes:

    PS> On Wed, May 20, 2009 at 11:10:11PM +0200, wolfgang.resch at gmail.com wrote:
    PS> ...
    >> Strange behavior of qbinom:
    >> 
    >> > qbinom(0.01, 5016279, 1e-07)
    >> [1] 0
    >> > qbinom(0.01, 5016279, 2e-07)
    >> [1] 16
    >> > qbinom(0.01, 5016279, 3e-07)
    >> [1] 16
    >> > qbinom(0.01, 5016279, 4e-07)
    >> [1] 16
    >> > qbinom(0.01, 5016279, 5e-07)
    >> [1] 0
    >> 

    PS> There is a bug in function do_search() in file src/nmath/qbinom.c. This
    PS> function contains a cycle
    PS> for(;;) {
    PS>   if(y == 0 ||
    PS>   (*z = pbinom(y - incr, n, pr, /*l._t.*/TRUE, /*log_p*/FALSE)) < p)
    PS>   return y;
    PS>   y = fmax2(0, y - incr);
    PS> }
    PS> When this cycle stops, *z contains pbinom(y - incr, ...), but is used
    PS> as if it is pbinom(y, ...) for the resulting y.

    PS> In the example qbinom(p=0.01, size=5016279, prob=4e-07), we get at 
    PS> some step y=15 as a result of a search left with incr=50, so we have
    PS> *z=pbinom(-35, ...)=0. Hence, y=15 is treated as too low and is increased
    PS> to 16. Since 16 is detected to be sufficient, the search stops with y=16,
    PS> which is wrong.

    [............]

Thanks to Wolfgang and Petr,

Petr's analysis of the problem seems right on spot to me,
and I'm currently testing the patch (and will also add some
regression tests along Petr's example)
which will make it into R-patched (to be R 2.9.1 in a while) and
R-devel.

Gratefully,
Martin Maechler, ETH Zurich