[Rd] qbinom (PR#13711)
maechler at stat.math.ethz.ch
maechler at stat.math.ethz.ch
Mon May 25 14:15:23 CEST 2009
>>>>> "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
More information about the R-devel
mailing list