On Wed, May 20, 2009 at 11:10:11PM +0200, wolfgang.re...@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.
--- 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 */
______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel