From 967289ddc051c05c9568307eb747842a1ab6dbdf Mon Sep 17 00:00:00 2001 From: maechler Date: Wed, 1 May 2024 08:43:29 +0000 Subject: [PATCH] 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 | 6 +++++- src/nmath/qDiscrete_search.h | 28 +++++++++++++++++----------- tests/d-p-q-r-tst-2.R | 16 ++++++++++++++++ 3 files changed, 38 insertions(+), 12 deletions(-) diff --git a/doc/NEWS.Rd b/doc/NEWS.Rd index 7ebcc743a61..2dbbd42d427 100644 --- a/doc/NEWS.Rd +++ b/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{[["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}. } } } diff --git a/src/nmath/qDiscrete_search.h b/src/nmath/qDiscrete_search.h index 0bdbcb49eea..df7bdf1ad6b 100644 --- a/src/nmath/qDiscrete_search.h +++ b/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) := (x, ), find p-quantile y(p) :<==> P(y) < p <= P(y) + For P(x) := (x, ), find p-quantile y = y(p) :<==> P(y-1) < p <= P(y) @param y current guess @param *z := (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 */ \ diff --git a/tests/d-p-q-r-tst-2.R b/tests/d-p-q-r-tst-2.R index 97df2e02e16..2c7fe4c0623 100644 --- a/tests/d-p-q-r-tst-2.R +++ b/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")