From af0b33745334b6b2623f5ae1ba5fdeae8764535b Mon Sep 17 00:00:00 2001 From: hiddevandebeek Date: Fri, 24 Jul 2026 13:10:03 +0200 Subject: [PATCH 1/6] Allow a continuous size (dispersion) in nbinom and nbinomMu In the negative binomial's mean/dispersion parameterisation `size` is a real dispersion parameter, not a count, and stats::dnbinom() has always accepted a continuous value. `size` was stored in an Eigen::VectorXi and assigned with (int)(size), so the C API entry points -- which model solves and focei reach directly, bypassing the R-level assertions -- silently returned the log-likelihood at trunc(size) for size > 1 and aborted the process for 0 < size < 1. The truncation also made the log-likelihood a step function of size, so the dispersion could not be estimated even when its true value was an integer. Store `size` as a double, replace the `size > INT_MAX` bound with `size <= 0` (NB2 requires a strictly positive dispersion, and the integer bound is meaningless for a double), and relax the R assertions to assertNumeric(). llikBinom() is unchanged: there `size` is a number of trials and the integer restriction is correct. Two tests asserting NA for `size > INT_MAX` encoded the old restriction and now assert the stats::dnbinom() value, which the new code matches exactly; new tests cover continuous size and its derivatives. --- NEWS.md | 14 ++++++++++++ R/llik.R | 12 +++++----- man/llikNbinom.Rd | 2 +- man/llikNbinomMu.Rd | 2 +- src/llikNbinom.cpp | 16 +++++++------- src/llikNbinom2.cpp | 12 +++++----- tests/testthat/test-llik.R | 45 +++++++++++++++++++++++++++++++++----- 7 files changed, 75 insertions(+), 28 deletions(-) diff --git a/NEWS.md b/NEWS.md index a9d374f..7e09dc7 100644 --- a/NEWS.md +++ b/NEWS.md @@ -3,6 +3,20 @@ * Fix ABI issues in rxode2 and the nlmixr2 ecosystem by requiring RcppParallel 6.0.0 +* `llikNbinom()` and `llikNbinomMu()` now accept a continuous (non-integer) + `size`. In the negative binomial's mean/dispersion parameterisation `size` + is a real dispersion parameter, not a count, and `stats::dnbinom()` has + always allowed it. Previously `size` was stored in an integer vector and + silently truncated, so a `size` above 1 returned the log-likelihood at + `trunc(size)` and a `size` between 0 and 1 truncated to 0 and aborted; the + R-level `assertIntegerish(size)` rejected such values outright. Truncation + also made the log-likelihood a step function of `size`, so the dispersion + could not be estimated even when its true value was an integer. This + blocked `nbinomMu()` models with continuous overdispersion in `nlmixr2`. + `size` must now be strictly positive; non-positive values return `NA`. The + `INT_MAX` bound added in 2.0.15 still applies to `x` everywhere, and to + `size` for `rxLlikBinom` where `size` is a number of trials. + # rxode2ll 2.0.15 * Fix signed integer overflow: loop indices in `llikXxxInternal()` Rcpp diff --git a/R/llik.R b/R/llik.R index 6b33805..29d3f38 100644 --- a/R/llik.R +++ b/R/llik.R @@ -88,9 +88,9 @@ llikBinom <- function(x, size, prob, full=FALSE) { } #' Calculate the log likelihood of the negative binomial function (and its derivatives) -#' +#' #' @param x Number of successes -#' @param size Size of trial +#' @param size Dispersion parameter; may be non-integer #' @param prob probability of success #' #' @inheritParams llikNorm @@ -108,7 +108,7 @@ llikBinom <- function(x, size, prob, full=FALSE) { #' llikNbinom <- function(x, size, prob, full=FALSE) { checkmate::assertIntegerish(x, min.len=0, lower=0, any.missing=FALSE) - checkmate::assertIntegerish(size, min.len=0, lower=0, any.missing=FALSE) + checkmate::assertNumeric(size, min.len=0, lower=0, any.missing=FALSE, finite=TRUE) checkmate::assertNumeric(prob, min.len=0, lower=0, upper=1, any.missing=FALSE, finite=TRUE) .df <- try(data.frame(x=x, size=size, prob=prob), silent=TRUE) if (inherits(.df, "try-error")) { @@ -123,8 +123,8 @@ llikNbinom <- function(x, size, prob, full=FALSE) { #' #' @param x Number of successes #' -#' @param size Size of trial -#' +#' @param size Dispersion parameter; may be non-integer +#' #' @param mu mu parameter for negative binomial #' #' @inheritParams llikNorm @@ -144,7 +144,7 @@ llikNbinom <- function(x, size, prob, full=FALSE) { #' llikNbinomMu <- function(x, size, mu, full=FALSE) { checkmate::assertIntegerish(x, min.len=0, lower=0, any.missing=FALSE) - checkmate::assertIntegerish(size, min.len=0, lower=0, any.missing=FALSE) + checkmate::assertNumeric(size, min.len=0, lower=0, any.missing=FALSE, finite=TRUE) checkmate::assertNumeric(mu, min.len=0, lower=0, any.missing=FALSE, finite=TRUE) .df <- try(data.frame(x=x, size=size, mu=mu), silent=TRUE) if (inherits(.df, "try-error")) { diff --git a/man/llikNbinom.Rd b/man/llikNbinom.Rd index 5c07827..b4b56b4 100644 --- a/man/llikNbinom.Rd +++ b/man/llikNbinom.Rd @@ -9,7 +9,7 @@ llikNbinom(x, size, prob, full = FALSE) \arguments{ \item{x}{Number of successes} -\item{size}{Size of trial} +\item{size}{Dispersion parameter; may be non-integer} \item{prob}{probability of success} diff --git a/man/llikNbinomMu.Rd b/man/llikNbinomMu.Rd index 94c2fec..19350c0 100644 --- a/man/llikNbinomMu.Rd +++ b/man/llikNbinomMu.Rd @@ -9,7 +9,7 @@ llikNbinomMu(x, size, mu, full = FALSE) \arguments{ \item{x}{Number of successes} -\item{size}{Size of trial} +\item{size}{Dispersion parameter; may be non-integer} \item{mu}{mu parameter for negative binomial} diff --git a/src/llikNbinom.cpp b/src/llikNbinom.cpp index db9075c..7504bc3 100644 --- a/src/llikNbinom.cpp +++ b/src/llikNbinom.cpp @@ -9,9 +9,9 @@ // R , sigma=scale struct nbinom_llik { const Eigen::VectorXi y_; - const Eigen::VectorXi N_; - - nbinom_llik(const Eigen::VectorXi& y, Eigen::VectorXi& N) : y_(y), N_(N) { } + const Eigen::VectorXd N_; // dispersion; real-valued, not a count + + nbinom_llik(const Eigen::VectorXi& y, Eigen::VectorXd& N) : y_(y), N_(N) { } template Eigen::Matrix operator()(const Eigen::Matrix& theta) const { @@ -19,7 +19,7 @@ struct nbinom_llik { Eigen::Matrix lp(y_.size()); for (Eigen::Index i = 0; i < y_.size(); ++i) { - T mu = (double)(N_[i])*(1.0-p)/p; + T mu = N_[i]*(1.0-p)/p; lp[i] = stan::math::neg_binomial_2_lpmf(y_[i], mu, N_[i]); } return lp; @@ -27,7 +27,7 @@ struct nbinom_llik { }; -stanLl llik_nbinom(Eigen::VectorXi& y, Eigen::VectorXi& N, Eigen::VectorXd& params) { +stanLl llik_nbinom(Eigen::VectorXi& y, Eigen::VectorXd& N, Eigen::VectorXd& params) { rx_stan_math_thread_init_rev_autodiff(); nbinom_llik f(y, N); Eigen::VectorXd fx; @@ -70,7 +70,7 @@ static inline void llikNbinomFull(double* ret, double x, double size, double pro return; } if (x < 0.0 || x > static_cast(INT_MAX) || - size < 0.0 || size > static_cast(INT_MAX)) { + size <= 0.0) { ret[0] = isNbinom; ret[1] = x; ret[2] = size; @@ -80,10 +80,10 @@ static inline void llikNbinomFull(double* ret, double x, double size, double pro return; } Eigen::VectorXi y(1); - Eigen::VectorXi N(1); + Eigen::VectorXd N(1); Eigen::VectorXd params(1); y(0) = (int)(x); - N(0) = (int)(size); + N(0) = size; params(0) = prob; stanLl ll = llik_nbinom(y, N, params); ret[0] = isNbinom; diff --git a/src/llikNbinom2.cpp b/src/llikNbinom2.cpp index 766f981..1150624 100644 --- a/src/llikNbinom2.cpp +++ b/src/llikNbinom2.cpp @@ -10,8 +10,8 @@ struct nbinomMu_llik { const Eigen::VectorXi y_; - const Eigen::VectorXi N_; - nbinomMu_llik(const Eigen::VectorXi& y, Eigen::VectorXi& N) : y_(y), N_(N) { } + const Eigen::VectorXd N_; // dispersion; real-valued, not a count + nbinomMu_llik(const Eigen::VectorXi& y, Eigen::VectorXd& N) : y_(y), N_(N) { } template Eigen::Matrix operator()(const Eigen::Matrix& theta) const { @@ -24,7 +24,7 @@ struct nbinomMu_llik { } }; -stanLl llik_nbinomMu(Eigen::VectorXi& y, Eigen::VectorXi& N, Eigen::VectorXd& params) { +stanLl llik_nbinomMu(Eigen::VectorXi& y, Eigen::VectorXd& N, Eigen::VectorXd& params) { rx_stan_math_thread_init_rev_autodiff(); nbinomMu_llik f(y, N); Eigen::VectorXd fx; @@ -68,7 +68,7 @@ static inline void llikNbinomMuFull(double* ret, double x, double size, double m return; } if (x < 0.0 || x > static_cast(INT_MAX) || - size < 0.0 || size > static_cast(INT_MAX)) { + size <= 0.0) { ret[0] = isNbinomMu; ret[1] = x; ret[2] = size; @@ -78,10 +78,10 @@ static inline void llikNbinomMuFull(double* ret, double x, double size, double m return; } Eigen::VectorXi y(1); - Eigen::VectorXi N(1); + Eigen::VectorXd N(1); Eigen::VectorXd params(1); y(0) = (int)(x); - N(0) = (int)(size); + N(0) = size; params(0) = mu; stanLl ll = llik_nbinomMu(y, N, params); ret[0] = isNbinomMu; diff --git a/tests/testthat/test-llik.R b/tests/testthat/test-llik.R index 8cb53e6..d1af9b6 100644 --- a/tests/testthat/test-llik.R +++ b/tests/testthat/test-llik.R @@ -39,6 +39,31 @@ test_that("log-liklihood tests for NbinomMu (including derivatives)", { expect_equal(fromR$fx, dnbinom(et$time, size=100, mu=40, log=TRUE)) }) +test_that("nbinom size may be continuous (non-integer dispersion)", { + # size used to be truncated to an int, silently returning the value at + # trunc(size) for size > 1 and aborting the process for 0 < size < 1. + size <- c(0.05, 0.3, 0.7, 0.99, 1.5, 2.5, 3.9) + x <- rep(2L, length(size)) + expect_equal(llikNbinomMu(x, size, rep(5, length(size)))$fx, + dnbinom(2L, size=size, mu=5, log=TRUE)) + expect_equal(llikNbinom(x, size, rep(0.4, length(size)))$fx, + dnbinom(2L, size=size, prob=0.4, log=TRUE)) +}) + +test_that("nbinom derivatives are correct for continuous size", { + h <- 1e-5 + size <- c(0.3, 0.7, 1.5, 4.25) + x <- rep(2L, length(size)) + dMuNum <- (dnbinom(2L, size=size, mu=5 + h, log=TRUE) - + dnbinom(2L, size=size, mu=5 - h, log=TRUE)) / (2 * h) + expect_equal(llikNbinomMu(x, size, rep(5, length(size)))$dMu, dMuNum, + tolerance=1e-5) + dProbNum <- (dnbinom(2L, size=size, prob=0.4 + h, log=TRUE) - + dnbinom(2L, size=size, prob=0.4 - h, log=TRUE)) / (2 * h) + expect_equal(llikNbinom(x, size, rep(0.4, length(size)))$dProb, dProbNum, + tolerance=1e-5) +}) + test_that("log-liklihood tests for beta (including derivatives)", { et <- data.frame(time=seq(1e-4, 1-1e-4, length.out=21)) et$shape1 <- 0.5 @@ -293,11 +318,15 @@ test_that("llikNbinomInternal returns NA for x > INT_MAX", { expect_true(is.na(res$dProb)) }) -test_that("llikNbinomInternal returns NA for size > INT_MAX", { +test_that("llikNbinomInternal accepts size > INT_MAX (size is a real dispersion)", { big <- 2^31 res <- llikNbinomInternal(5, big, 0.5) - expect_true(is.na(res$fx)) - expect_true(is.na(res$dProb)) + expect_equal(res$fx, stats::dnbinom(5, size = big, prob = 0.5, log = TRUE)) +}) + +test_that("llikNbinomInternal returns NA for size <= 0", { + expect_true(is.na(llikNbinomInternal(5, 0, 0.5)$fx)) + expect_true(is.na(llikNbinomInternal(5, -1, 0.5)$fx)) }) test_that("llikNbinomMuInternal returns NA for x > INT_MAX", { @@ -307,11 +336,15 @@ test_that("llikNbinomMuInternal returns NA for x > INT_MAX", { expect_true(is.na(res$dMu)) }) -test_that("llikNbinomMuInternal returns NA for size > INT_MAX", { +test_that("llikNbinomMuInternal accepts size > INT_MAX (size is a real dispersion)", { big <- 2^31 res <- llikNbinomMuInternal(5, big, 40) - expect_true(is.na(res$fx)) - expect_true(is.na(res$dMu)) + expect_equal(res$fx, stats::dnbinom(5, size = big, mu = 40, log = TRUE)) +}) + +test_that("llikNbinomMuInternal returns NA for size <= 0", { + expect_true(is.na(llikNbinomMuInternal(5, 0, 40)$fx)) + expect_true(is.na(llikNbinomMuInternal(5, -1, 40)$fx)) }) ## Large-vector test for R_xlen_t fix From b00650d2ce9bd3ffdd2fe0653574f745ceff0a7a Mon Sep 17 00:00:00 2001 From: mattfidler Date: Sat, 25 Jul 2026 15:45:58 -0500 Subject: [PATCH 2/6] Guard the nbinom mean so out-of-domain inputs return NA, not an abort `neg_binomial_2_lpmf()` requires a positive finite mean. Both nbinom entry points could reach it with an out-of-domain mean, and because `rxLlikNbinom()` / `rxLlikNbinomMu()` are `extern "C"` -- called directly by rxode2 solves and focei, with no handler in between -- the resulting C++ exception escaped and aborted the R process. Dropping the `size > INT_MAX` bound in the size/prob form opened a new case: `mu = size*(1-prob)/prob` can now overflow to `Inf` for a large finite `size` with a small `prob`, which the old bound had rejected. Guard the derived mean itself rather than re-bounding `size`, which also covers the pre-existing `prob == 0` (`mu = Inf`), `prob == 1` (`mu = 0`) and, via the C API, `prob` outside `[0, 1]` (`mu < 0`). The mean/dispersion form has the same hole for a non-positive `mu`; guard that too. Tests cover each case, and NEWS.md records the new `NA` returns. Co-Authored-By: Claude Opus 5 (1M context) --- NEWS.md | 9 +++++++++ src/llikNbinom.cpp | 9 ++++++++- src/llikNbinom2.cpp | 5 ++++- tests/testthat/test-llik.R | 23 +++++++++++++++++++++++ 4 files changed, 44 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index 7e09dc7..ab9172f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -17,6 +17,15 @@ `INT_MAX` bound added in 2.0.15 still applies to `x` everywhere, and to `size` for `rxLlikBinom` where `size` is a number of trials. +* `llikNbinom()` and `llikNbinomMu()` now return `NA` instead of aborting the R + process when the mean handed to the underlying likelihood is out of domain. + This covers `prob` of 0 or 1 (and, through the C API, `prob` outside + `[0, 1]`), a non-positive `mu`, and -- because `size` is no longer bounded by + `INT_MAX` -- a large `size` with a small `prob` whose product overflows. + Previously these threw a C++ exception that escaped the `rxLlikNbinom()` / + `rxLlikNbinomMu()` entry points used by `rxode2` solves and `focei`, killing + the session. + # rxode2ll 2.0.15 * Fix signed integer overflow: loop indices in `llikXxxInternal()` Rcpp diff --git a/src/llikNbinom.cpp b/src/llikNbinom.cpp index 7504bc3..abe61bd 100644 --- a/src/llikNbinom.cpp +++ b/src/llikNbinom.cpp @@ -69,8 +69,15 @@ static inline void llikNbinomFull(double* ret, double x, double size, double pro ret[5] = NA_REAL; return; } + // neg_binomial_2_lpmf() needs a positive finite mean, which this + // parameterisation reaches through mu = size*(1-prob)/prob: that is Inf at + // prob == 0, 0 at prob == 1, negative outside [0, 1], and -- now that size is + // no longer bounded by INT_MAX -- can overflow for a large size with a small + // prob. Stan throws on all of those, and the exception would escape the + // extern "C" entry points below and abort the R process, so return NA here. + double meanPar = size*(1.0-prob)/prob; if (x < 0.0 || x > static_cast(INT_MAX) || - size <= 0.0) { + size <= 0.0 || !R_finite(meanPar) || meanPar <= 0.0) { ret[0] = isNbinom; ret[1] = x; ret[2] = size; diff --git a/src/llikNbinom2.cpp b/src/llikNbinom2.cpp index 1150624..e7f58da 100644 --- a/src/llikNbinom2.cpp +++ b/src/llikNbinom2.cpp @@ -67,8 +67,11 @@ static inline void llikNbinomMuFull(double* ret, double x, double size, double m ret[5] = NA_REAL; return; } + // neg_binomial_2_lpmf() needs a positive finite mean; a non-positive mu makes + // it throw, and that exception would escape the extern "C" entry points below + // and abort the R process, so return NA here instead. if (x < 0.0 || x > static_cast(INT_MAX) || - size <= 0.0) { + size <= 0.0 || mu <= 0.0) { ret[0] = isNbinomMu; ret[1] = x; ret[2] = size; diff --git a/tests/testthat/test-llik.R b/tests/testthat/test-llik.R index d1af9b6..9907785 100644 --- a/tests/testthat/test-llik.R +++ b/tests/testthat/test-llik.R @@ -329,6 +329,22 @@ test_that("llikNbinomInternal returns NA for size <= 0", { expect_true(is.na(llikNbinomInternal(5, -1, 0.5)$fx)) }) +test_that("llikNbinomInternal returns NA for an out-of-domain prob", { + # mu = size*(1-prob)/prob is Inf at prob == 0, 0 at prob == 1 and negative + # outside [0, 1]; each makes neg_binomial_2_lpmf() throw, which aborts the + # process when it escapes rxLlikNbinom() + expect_true(is.na(llikNbinomInternal(5, 10, 0)$fx)) + expect_true(is.na(llikNbinomInternal(5, 10, 1)$fx)) + expect_true(is.na(llikNbinomInternal(5, 10, -0.5)$fx)) + expect_true(is.na(llikNbinomInternal(5, 10, 1.5)$fx)) +}) + +test_that("llikNbinomInternal returns NA when size/prob overflows the mean", { + # size is no longer bounded by INT_MAX, so size*(1-prob)/prob can overflow + # even though both size and prob are finite and in range + expect_true(is.na(llikNbinomInternal(5, 1e300, 1e-10)$fx)) +}) + test_that("llikNbinomMuInternal returns NA for x > INT_MAX", { big <- 2^31 res <- llikNbinomMuInternal(big, 10, 40) @@ -347,6 +363,13 @@ test_that("llikNbinomMuInternal returns NA for size <= 0", { expect_true(is.na(llikNbinomMuInternal(5, -1, 40)$fx)) }) +test_that("llikNbinomMuInternal returns NA for mu <= 0", { + # neg_binomial_2_lpmf() throws on a non-positive mean, which aborts the + # process when it escapes rxLlikNbinomMu() + expect_true(is.na(llikNbinomMuInternal(5, 10, 0)$fx)) + expect_true(is.na(llikNbinomMuInternal(5, 10, -1)$fx)) +}) + ## Large-vector test for R_xlen_t fix ## Skipped in normal test runs: requires ~17 GB per numeric vector ## (~103 GB total for 6 vectors), which exceeds typical system RAM. From 5796d8ba3e5f9410715d84f904e9e38d9708ac26 Mon Sep 17 00:00:00 2001 From: mattfidler Date: Sat, 25 Jul 2026 16:36:42 -0500 Subject: [PATCH 3/6] drop RcppParallel requirement --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 8522a8b..a0fd84f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -14,7 +14,7 @@ Suggests: Imports: Rcpp (>= 1.0.8), checkmate, - RcppParallel (>= 6.0.0) + RcppParallel Description: Provides the log-likelihoods with gradients from 'stan' (Carpenter et al (2015), ) needed for generalized log-likelihood estimation in 'nlmixr2' From 6cb51c256906aebb1a2d778f4d20cfef39867c8e Mon Sep 17 00:00:00 2001 From: mattfidler Date: Sat, 25 Jul 2026 16:38:27 -0500 Subject: [PATCH 4/6] Increment version number to 2.0.16.9000 --- DESCRIPTION | 2 +- NEWS.md | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index a0fd84f..0b1ae9f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,5 +1,5 @@ Package: rxode2ll -Version: 2.0.16 +Version: 2.0.16.9000 Title: Log-Likelihood Functions for 'rxode2' Authors@R: c( person("Matthew L.","Fidler", diff --git a/NEWS.md b/NEWS.md index ab9172f..30f19d9 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,5 @@ +# rxode2ll (development version) + # rxode2ll 2.0.16 * Fix ABI issues in rxode2 and the nlmixr2 ecosystem by requiring From 05691e83714e290c0613bd84d9738ceb65bcee98 Mon Sep 17 00:00:00 2001 From: mattfidler Date: Sat, 25 Jul 2026 16:39:22 -0500 Subject: [PATCH 5/6] fix news --- NEWS.md | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/NEWS.md b/NEWS.md index 30f19d9..dbbcc45 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,10 +1,5 @@ # rxode2ll (development version) -# rxode2ll 2.0.16 - -* Fix ABI issues in rxode2 and the nlmixr2 ecosystem by requiring - RcppParallel 6.0.0 - * `llikNbinom()` and `llikNbinomMu()` now accept a continuous (non-integer) `size`. In the negative binomial's mean/dispersion parameterisation `size` is a real dispersion parameter, not a count, and `stats::dnbinom()` has @@ -28,6 +23,12 @@ `rxLlikNbinomMu()` entry points used by `rxode2` solves and `focei`, killing the session. + +# rxode2ll 2.0.16 + +* Fix ABI issues in rxode2 and the nlmixr2 ecosystem by requiring + RcppParallel 6.0.0 + # rxode2ll 2.0.15 * Fix signed integer overflow: loop indices in `llikXxxInternal()` Rcpp From 7f5389538582e4adef6e03506d9907930b3ca129 Mon Sep 17 00:00:00 2001 From: mattfidler Date: Sat, 25 Jul 2026 18:29:02 -0500 Subject: [PATCH 6/6] Return R's likelihood at the degenerate nbinom points At prob == 1 (size/prob form) and mu == 0 (mean/dispersion form) the negative binomial collapses to a point mass at zero, and stats::dnbinom() gives a likelihood of 1 at x == 0 and 0 elsewhere. neg_binomial_2_lpmf() needs a strictly positive mean, so Stan cannot evaluate these points at all; the previous commit made them return NA to stop the C++ exception escaping rxLlikNbinom() / rxLlikNbinomMu() and aborting the R process. Fill in the value R would instead -- 0 at x == 0, -Inf otherwise -- and report the derivative, which does not exist at a point mass, as NA. An unusable gradient is better than an unusable log-likelihood: rxode2ll now matches R at every point where R defines one. Key the size/prob branch on prob == 1 rather than on the derived mean == 0. Those are not the same condition: size*(1-prob)/prob also underflows to exactly 0 for a denormal size with a prob well under 1, where the true log-likelihood is finite (about -745) rather than -Inf. That case stays NA, as before. Tests cover both degenerate points -- against stats::dnbinom(), for a continuous size, and through the R-level interface -- plus the denormal underflow that must not take the degenerate path. Co-Authored-By: Claude Opus 5 (1M context) --- NEWS.md | 17 +++++++++--- src/llikNbinom.cpp | 38 ++++++++++++++++++++++----- src/llikNbinom2.cpp | 24 ++++++++++++++--- tests/testthat/test-llik.R | 54 ++++++++++++++++++++++++++++++++------ 4 files changed, 111 insertions(+), 22 deletions(-) diff --git a/NEWS.md b/NEWS.md index dbbcc45..b9e5548 100644 --- a/NEWS.md +++ b/NEWS.md @@ -16,13 +16,22 @@ * `llikNbinom()` and `llikNbinomMu()` now return `NA` instead of aborting the R process when the mean handed to the underlying likelihood is out of domain. - This covers `prob` of 0 or 1 (and, through the C API, `prob` outside - `[0, 1]`), a non-positive `mu`, and -- because `size` is no longer bounded by - `INT_MAX` -- a large `size` with a small `prob` whose product overflows. - Previously these threw a C++ exception that escaped the `rxLlikNbinom()` / + This covers a `prob` of 0 (and, through the C API, a `prob` outside `[0, 1]` + or a negative `mu`), and -- because `size` is no longer bounded by `INT_MAX` + -- a large `size` with a small `prob` whose product overflows. Previously + these threw a C++ exception that escaped the `rxLlikNbinom()` / `rxLlikNbinomMu()` entry points used by `rxode2` solves and `focei`, killing the session. +* `llikNbinom()` at `prob == 1` and `llikNbinomMu()` at `mu == 0` now match + `stats::dnbinom()`. The distribution is degenerate at zero there, so `fx` is + `0` (a likelihood of 1) at `x == 0` and `-Inf` otherwise. These are the only + points where `stats::dnbinom()` defines a value that the underlying Stan + likelihood cannot compute -- it requires a strictly positive mean. No + derivative exists at a point mass, so `dProb` / `dMu` are `NA`; an unusable + gradient beats a `NA` log-likelihood, and beats the abort these used to + cause. `prob == 0` remains `NA`, matching the `NaN` R returns there. + # rxode2ll 2.0.16 diff --git a/src/llikNbinom.cpp b/src/llikNbinom.cpp index abe61bd..1f4d174 100644 --- a/src/llikNbinom.cpp +++ b/src/llikNbinom.cpp @@ -69,15 +69,41 @@ static inline void llikNbinomFull(double* ret, double x, double size, double pro ret[5] = NA_REAL; return; } + if (x < 0.0 || x > static_cast(INT_MAX) || size <= 0.0) { + ret[0] = isNbinom; + ret[1] = x; + ret[2] = size; + ret[3] = prob; + ret[4] = NA_REAL; + ret[5] = NA_REAL; + return; + } + if (prob == 1.0) { + // prob == 1 makes the distribution degenerate at zero, so the likelihood is + // 1 at x == 0 and 0 elsewhere, matching stats::dnbinom(). Stan cannot + // evaluate this point at all (it needs a strictly positive mean) and there + // is no usable derivative here, so report the value R would and leave the + // gradient NA rather than returning NA for both. Test prob rather than the + // derived mean below: that mean can also underflow to 0 for a denormal size + // with a prob well under 1, where the true log-likelihood is finite and + // large-negative rather than -Inf. + ret[0] = isNbinom; + ret[1] = x; + ret[2] = size; + ret[3] = prob; + ret[4] = (x == 0.0) ? 0.0 : R_NegInf; + ret[5] = NA_REAL; + return; + } // neg_binomial_2_lpmf() needs a positive finite mean, which this // parameterisation reaches through mu = size*(1-prob)/prob: that is Inf at - // prob == 0, 0 at prob == 1, negative outside [0, 1], and -- now that size is - // no longer bounded by INT_MAX -- can overflow for a large size with a small - // prob. Stan throws on all of those, and the exception would escape the - // extern "C" entry points below and abort the R process, so return NA here. + // prob == 0, negative outside [0, 1], and -- now that size is no longer + // bounded by INT_MAX -- can overflow for a large size with a small prob, or + // underflow to 0 for a denormal size. Stan throws on all of those, and the + // exception would escape the extern "C" entry points below and abort the R + // process, so return NA here. double meanPar = size*(1.0-prob)/prob; - if (x < 0.0 || x > static_cast(INT_MAX) || - size <= 0.0 || !R_finite(meanPar) || meanPar <= 0.0) { + if (!R_finite(meanPar) || meanPar <= 0.0) { ret[0] = isNbinom; ret[1] = x; ret[2] = size; diff --git a/src/llikNbinom2.cpp b/src/llikNbinom2.cpp index e7f58da..b05c7fd 100644 --- a/src/llikNbinom2.cpp +++ b/src/llikNbinom2.cpp @@ -67,11 +67,13 @@ static inline void llikNbinomMuFull(double* ret, double x, double size, double m ret[5] = NA_REAL; return; } - // neg_binomial_2_lpmf() needs a positive finite mean; a non-positive mu makes - // it throw, and that exception would escape the extern "C" entry points below - // and abort the R process, so return NA here instead. + // neg_binomial_2_lpmf() needs a positive finite mean; a negative mu (reachable + // through the C API, which bypasses the R-level assertions) makes it throw, + // and that exception would escape the extern "C" entry points below and abort + // the R process, so return NA here instead. mu == 0 is degenerate rather than + // out of domain and is handled just below. if (x < 0.0 || x > static_cast(INT_MAX) || - size <= 0.0 || mu <= 0.0) { + size <= 0.0 || mu < 0.0) { ret[0] = isNbinomMu; ret[1] = x; ret[2] = size; @@ -80,6 +82,20 @@ static inline void llikNbinomMuFull(double* ret, double x, double size, double m ret[5] = NA_REAL; return; } + if (mu == 0.0) { + // mu == 0 makes the distribution degenerate at zero, so the likelihood is 1 + // at x == 0 and 0 elsewhere, matching stats::dnbinom(). Stan cannot + // evaluate this point at all (it needs a strictly positive mean) and there + // is no usable derivative here, so report the value R would and leave the + // gradient NA rather than returning NA for both. + ret[0] = isNbinomMu; + ret[1] = x; + ret[2] = size; + ret[3] = mu; + ret[4] = (x == 0.0) ? 0.0 : R_NegInf; + ret[5] = NA_REAL; + return; + } Eigen::VectorXi y(1); Eigen::VectorXd N(1); Eigen::VectorXd params(1); diff --git a/tests/testthat/test-llik.R b/tests/testthat/test-llik.R index 9907785..a48d970 100644 --- a/tests/testthat/test-llik.R +++ b/tests/testthat/test-llik.R @@ -330,21 +330,45 @@ test_that("llikNbinomInternal returns NA for size <= 0", { }) test_that("llikNbinomInternal returns NA for an out-of-domain prob", { - # mu = size*(1-prob)/prob is Inf at prob == 0, 0 at prob == 1 and negative - # outside [0, 1]; each makes neg_binomial_2_lpmf() throw, which aborts the - # process when it escapes rxLlikNbinom() + # mu = size*(1-prob)/prob is Inf at prob == 0 and negative outside [0, 1]; + # each makes neg_binomial_2_lpmf() throw, which aborts the process when it + # escapes rxLlikNbinom() expect_true(is.na(llikNbinomInternal(5, 10, 0)$fx)) - expect_true(is.na(llikNbinomInternal(5, 10, 1)$fx)) expect_true(is.na(llikNbinomInternal(5, 10, -0.5)$fx)) expect_true(is.na(llikNbinomInternal(5, 10, 1.5)$fx)) }) +test_that("llikNbinom is degenerate at zero for prob == 1", { + # prob == 1 puts all the mass at zero, as stats::dnbinom() does; Stan cannot + # evaluate the point (mu == 0), so the value is filled in and the gradient, + # which does not exist here, is NA + res <- llikNbinomInternal(c(0, 1, 5), rep(10, 3), rep(1, 3)) + expect_equal(res$fx, stats::dnbinom(c(0, 1, 5), size = 10, prob = 1, log = TRUE)) + expect_equal(res$fx, c(0, -Inf, -Inf)) + expect_true(all(is.na(res$dProb))) + # a continuous size is degenerate at zero too + expect_equal(llikNbinomInternal(c(0, 1), c(0.5, 0.5), c(1, 1))$fx, + stats::dnbinom(c(0, 1), size = 0.5, prob = 1, log = TRUE)) + # and it is reachable from the R-level interface + expect_equal(llikNbinom(0L, 10, 1)$fx, 0) +}) + test_that("llikNbinomInternal returns NA when size/prob overflows the mean", { # size is no longer bounded by INT_MAX, so size*(1-prob)/prob can overflow # even though both size and prob are finite and in range expect_true(is.na(llikNbinomInternal(5, 1e300, 1e-10)$fx)) }) +test_that("llikNbinomInternal returns NA when a denormal size underflows the mean", { + # size*(1-prob)/prob underflows to exactly 0 for a denormal size even with a + # prob well under 1. That is not the degenerate prob == 1 point -- the true + # log-likelihood is finite (about -745 here) -- so it must return NA rather + # than the -Inf used at the point mass + expect_true(is.finite(stats::dnbinom(1, size = 5e-324, prob = 0.5, log = TRUE))) + expect_true(is.na(llikNbinomInternal(1, 5e-324, 0.5)$fx)) + expect_true(is.na(llikNbinomInternal(0, 5e-324, 0.5)$fx)) +}) + test_that("llikNbinomMuInternal returns NA for x > INT_MAX", { big <- 2^31 res <- llikNbinomMuInternal(big, 10, 40) @@ -363,13 +387,27 @@ test_that("llikNbinomMuInternal returns NA for size <= 0", { expect_true(is.na(llikNbinomMuInternal(5, -1, 40)$fx)) }) -test_that("llikNbinomMuInternal returns NA for mu <= 0", { - # neg_binomial_2_lpmf() throws on a non-positive mean, which aborts the - # process when it escapes rxLlikNbinomMu() - expect_true(is.na(llikNbinomMuInternal(5, 10, 0)$fx)) +test_that("llikNbinomMuInternal returns NA for mu < 0", { + # neg_binomial_2_lpmf() throws on a negative mean, which aborts the process + # when it escapes rxLlikNbinomMu() expect_true(is.na(llikNbinomMuInternal(5, 10, -1)$fx)) }) +test_that("llikNbinomMu is degenerate at zero for mu == 0", { + # mu == 0 puts all the mass at zero, as stats::dnbinom() does; Stan cannot + # evaluate the point, so the value is filled in and the gradient, which does + # not exist here, is NA + res <- llikNbinomMuInternal(c(0, 1, 5), rep(10, 3), rep(0, 3)) + expect_equal(res$fx, stats::dnbinom(c(0, 1, 5), size = 10, mu = 0, log = TRUE)) + expect_equal(res$fx, c(0, -Inf, -Inf)) + expect_true(all(is.na(res$dMu))) + # a continuous size is degenerate at zero too + expect_equal(llikNbinomMuInternal(c(0, 1), c(0.5, 0.5), c(0, 0))$fx, + stats::dnbinom(c(0, 1), size = 0.5, mu = 0, log = TRUE)) + # and it is reachable from the R-level interface + expect_equal(llikNbinomMu(0L, 10, 0)$fx, 0) +}) + ## Large-vector test for R_xlen_t fix ## Skipped in normal test runs: requires ~17 GB per numeric vector ## (~103 GB total for 6 vectors), which exceeds typical system RAM.