diff --git a/DESCRIPTION b/DESCRIPTION index 8522a8b..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", @@ -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' diff --git a/NEWS.md b/NEWS.md index a9d374f..b9e5548 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,38 @@ +# rxode2ll (development version) + +* `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. + +* `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 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 * Fix ABI issues in rxode2 and the nlmixr2 ecosystem by requiring 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..1f4d174 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; @@ -69,8 +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 || size > static_cast(INT_MAX)) { + 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, 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 (!R_finite(meanPar) || meanPar <= 0.0) { ret[0] = isNbinom; ret[1] = x; ret[2] = size; @@ -80,10 +113,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..b05c7fd 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; @@ -67,8 +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 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 || size > static_cast(INT_MAX)) { + size <= 0.0 || mu < 0.0) { ret[0] = isNbinomMu; ret[1] = x; ret[2] = size; @@ -77,11 +82,25 @@ 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::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..a48d970 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,55 @@ 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("llikNbinomInternal returns NA for an out-of-domain prob", { + # 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, -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", { @@ -307,11 +376,36 @@ 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)) +}) + +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