[Rd] qbinom (PR#13711)

[Petr Savicky]

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 */