# [Rd] qbinom (PR#13711)

Petr Savicky savicky at cs.cas.cz
Sat May 23 18:22:26 CEST 2009

```On Wed, May 20, 2009 at 11:10:11PM +0200, wolfgang.resch at gmail.com wrote:
...
> 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
>

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

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

A possible correction is to replace the code above by
double newz;
for(;;) {
if(y == 0 ||
(newz = pbinom(y - incr, n, pr, /*l._t.*/TRUE, /*log_p*/FALSE)) < p)
return y;
y = fmax2(0, y - incr);
*z = newz;
}

Patch against R-devel_2009-05-22 is in an attachment.

Verification may be done using the following examples.
b1 <- qbinom(p=0.01, size=5016279, prob=(1:20)*1e-07)
b2 <- qbinom(p=0.01, size=5006279, prob=(1:20)*1e-07)
b3 <- qbinom(p=0.01, size=5000279, prob=(1:20)*1e-07)
p1 <- qpois(p=0.01, lambda=5016279*(1:20)*1e-07)
p2 <- qpois(p=0.01, lambda=5006279*(1:20)*1e-07)
p3 <- qpois(p=0.01, lambda=5000279*(1:20)*1e-07)
cbind(b1, b2, b3, p1, p2, p3)

Wrong
b1 b2 b3 p1 p2 p3
[1,]  0  0  0  0  0  0
[2,] 16  6 50  0  0  0
[3,] 16  6 50  0  0  0
[4,] 16  6 50  0  0  0
[5,]  0  0  0  0  0  0
[6,]  0  0  0  0  0  0
[7,]  0  0  0  0  0  0
[8,]  0  0  0  0  0  0
[9,]  0  0  0  0  0  0
[10,]  1  1  1  1  1  1
[11,]  1  1  1  1  1  1
[12,]  1  1  1  1  1  1
[13,]  1  1  1  1  1  1
[14,]  2  2  2  2  2  2
[15,]  2  2  2  2  2  2
[16,]  2  2  2  2  2  2
[17,] 19  9 53  3  3  3
[18,]  3  3  3  3  3  3
[19,]  3  3  3  3  3  3
[20,]  3  3  3  3  3  3

Correct
b1 b2 b3 p1 p2 p3
[1,]  0  0  0  0  0  0
[2,]  0  0  0  0  0  0
[3,]  0  0  0  0  0  0
[4,]  0  0  0  0  0  0
[5,]  0  0  0  0  0  0
[6,]  0  0  0  0  0  0
[7,]  0  0  0  0  0  0
[8,]  0  0  0  0  0  0
[9,]  0  0  0  0  0  0
[10,]  1  1  1  1  1  1
[11,]  1  1  1  1  1  1
[12,]  1  1  1  1  1  1
[13,]  1  1  1  1  1  1
[14,]  2  2  2  2  2  2
[15,]  2  2  2  2  2  2
[16,]  2  2  2  2  2  2
[17,]  3  3  3  3  3  3
[18,]  3  3  3  3  3  3
[19,]  3  3  3  3  3  3
[20,]  3  3  3  3  3  3

Petr.

-------------- next part --------------
--- R-devel/src/nmath/qbinom.c	2007-07-25 17:54:27.000000000 +0200
+++ R-devel-qbinom/src/nmath/qbinom.c	2009-05-23 17:22:43.538566976 +0200
@@ -36,6 +36,7 @@
static double
do_search(double y, double *z, double p, double n, double pr, double incr)
{
+    double newz;
if(*z >= p) {
/* search to the left */
#ifdef DEBUG_qbinom
@@ -43,9 +44,10 @@
#endif
for(;;) {
if(y == 0 ||
-	       (*z = pbinom(y - incr, n, pr, /*l._t.*/TRUE, /*log_p*/FALSE)) < p)
+	       (newz = pbinom(y - incr, n, pr, /*l._t.*/TRUE, /*log_p*/FALSE)) < p)
return y;
y = fmax2(0, y - incr);
+	    *z = newz;
}
}
else {		/* search to the right */
```