From maechler at stat.math.ethz.ch  Mon May 25 14:15:23 2009
From: maechler at stat.math.ethz.ch (maechler at stat.math.ethz.ch)
Date: Mon, 25 May 2009 14:15:23 +0200 (CEST)
Subject: [Rd] qbinom (PR#13711)
Message-ID: <20090525121523.9EB0A282C765@mail.pubhealth.ku.dk>

>>>>> "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



[Rd] qbinom (PR#13711)
maechler at stat.math.ethz.ch