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