fix discrete inversion p |--> q; for qbinom() and more: PR#18711

git-svn-id: https://svn.r-project.org/R/trunk@86504 00db46b3-68df-0310-9c12-caf00c1e9a41

doc/NEWS.Rd

@@ -49,7 +49,7 @@
   \subsection{C-LEVEL FACILITIES}{
     \itemize{
       \item The non-API and hidden entry points \code{Rf_setIVector},
-      \code{Rf_setRVector} and \code{Rf_setSVector} have been removed.     
+      \code{Rf_setRVector} and \code{Rf_setSVector} have been removed.
     }
   }
 
@@ -89,6 +89,10 @@
       \item The error message from \code{<POSIXlt>[["hour"]]} and similar
       now mentions \code{*[[, "hour"]]}, as wished for in \PR{17409} and
       proposed by \I{Michael Chirico}.
+
+      \item \code{qbinom()} and potentially \code{qpois()}, \code{qnbinom()},
+      no longer sometimes fail accurate inversion (of \code{pbinom()},
+      etc), thanks to Christopher Chang's report and patch in \PR{18711}.
     }
   }
 }

src/nmath/qDiscrete_search.h

@@ -1,6 +1,6 @@
 /*
  *  Mathlib : A C Library of Special Functions
- *  Copyright (C) 2000-2021 The R Core Team
+ *  Copyright (C) 2000-2024 The R Core Team
  *  Copyright (C) 2005-2021 The R Foundation
  *
  *  This program is free software; you can redistribute it and/or modify
@@ -18,7 +18,7 @@
  *  https://www.R-project.org/Licenses/
  */
 
-/* This is #included  from ./qnbinom.c , .........
+/* This is #included  from ./qpois.c ./qbinom.c, ./qnbinom{,_mu}.c
 */
 
 #define PST_0(a, b) a ## b
@@ -31,7 +31,7 @@
 #define _qDIST_  PASTE(q, _thisDIST_)
 /**
    For a discrete distribution on  the integers,
-   For P(x) := <pDist>(x, <distPars>),  find p-quantile  y(p)   :<==>   P(y) < p <= P(y)
+   For P(x) := <pDist>(x, <distPars>),  find p-quantile  y = y(p)   :<==>   P(y-1) < p <= P(y)
 
    @param y    current guess
    @param *z   := <pDist>(y, ..)
@@ -90,38 +90,44 @@ static double DO_SEARCH_FUN(_dist_PARS_DECL_)
     }
     else { // (lower_tail, *z < p)  or  (upper tail, *z >= p): search to the __right__
 	for(int iter = 0; ; iter++) {
+	    double prevy = y;
+	    double newz = -1.; // -Wall
 #ifndef MATHLIB_STANDALONE
 	    if(iter % 10000 == 0) R_CheckUserInterrupt();
 #endif
 	    y += incr;
 #ifdef _dist_MAX_y
 	    if(y < _dist_MAX_y)
-		*z = P_DIST(y, _dist_PARS_);
+		newz = P_DIST(y, _dist_PARS_);
 	    else if(y > _dist_MAX_y)
 		y = _dist_MAX_y;
 #else
-	    *z = P_DIST(y, _dist_PARS_);
+	    newz = P_DIST(y, _dist_PARS_);
 #endif
 
 	    if(
 #ifdef _dist_MAX_y
 		y == _dist_MAX_y ||
 #endif
-		ISNAN(*z) || (lower_tail ? (*z >= p) : (*z < p)))
+		ISNAN(newz) || (lower_tail ? (newz >= p) : (newz < p)))
 	    {
 		R_DBG_printf("  new y=%.15g, z=%g = " AS_CHAR(_pDIST_) "(y,*) %s;"
-			     " ==> search() returns after %d iter.\n", y, *z,
-			     ISNAN(*z) ? "is NaN" : (lower_tail ? ">= p" : "< p"), iter);
-		return y;
+			     " ==> search() returns after %d iter.\n", y, newz,
+			     ISNAN(newz) ? "is NaN" : (lower_tail ? ">= p" : "< p"), iter);
+		if (incr <= 1) {
+		    *z = newz;
+		    return y;
+		}
+		return prevy;
 	    }
+	    *z = newz;
 	}
     }
 } // do_search()
 
 
 /*
-* Note : "same" code in qbinom.c, qnbinom.c __FIXME__ NOT YET for qpois() ??
-* FIXME: This is far from optimal [cancellation for p ~= 1, etc]:
+* Note : called in qbinom.c, qnbinom.c but not (yet) qpois.c -- NB: only DBG_print()ing; *no other* effect
 */
 #define q_DISCRETE_01_CHECKS() do {					\
     double p_n; /* p  in the "normal" (lower_tail=TRUE, log.p=FALSE) scale */ \

tests/d-p-q-r-tst-2.R

@@ -856,6 +856,7 @@ stopifnot(exprs = {
     pbeta(1,  0, 0) == 1   # gave 0.5
 })
 
+
 ## PR#18640 -- stirlerr(x) concerns (for *non* half-integer x) -- visible in dgamma()
 sh <- 465/32 # = 14.53125
 x0 <- 1/4 + 8:20
@@ -870,5 +871,20 @@ relE * 2^53 # 2  2 -2  0  0  2 -1  0  0  2  0  0  2 //  was in {-95 : -91}  in R
 stopifnot(abs(relE) < 9e-16) # max{x86_64}: 2.22e-16
 
 
+## PR#18711 -- qbinom() - inversion of pbinom()
+##             but probably also fixing qpois(), qnbinom() cases
+sz <- 6040:6045
+prb <- 0.995
+(qb6 <- qbinom(p = 0.05, size = sz, prob = prb))
+(pqb6   <- pbinom(qb6,   size = sz, prob = prb))
+(pqb6_1 <- pbinom(qb6-1, size = sz, prob = prb))
+stopifnot(exprs = {
+    qb6 == c(6001:6004,6004:6005) # not in R 4.4.0, nor 4.1.0
+    1 > pqb6 & pqb6 >= 0.05       #  "
+X    0.05 > pqb6_1 & pqb6_1 >= 0.035# "
+})
+## was wrong for R versions in [4.1.1, 4.4.0]
+
+
 
 cat("Time elapsed: ", proc.time() - .ptime,"\n")