version 1.0.1

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: nbconv
 Title: Evaluate Arbitrary Negative Binomial Convolutions
-Version: 1.0.0
+Version: 1.0.1
 Authors@R: 
     person(
       "Gregory", "Bedwell",
@@ -15,9 +15,9 @@ Encoding: UTF-8
 RoxygenNote: 7.2.0
 License: GPL (>= 3)
 NeedsCompilation: no
-Packaged: 2023-04-26 12:49:11 UTC; gbedwell
+Packaged: 2023-07-06 18:01:40 UTC; gbedwell
 Author: Gregory Bedwell [aut, cre, cph]
     (<https://orcid.org/0000-0003-0456-0032>)
 Maintainer: Gregory Bedwell <gregoryjbedwell@gmail.com>
 Repository: CRAN
-Date/Publication: 2023-04-27 12:50:02 UTC
+Date/Publication: 2023-07-06 18:20:02 UTC

MD5

@@ -1,20 +1,20 @@
-3075f708f3a4c1fc5edc2fc10deb16db *DESCRIPTION
+700019b87a6cd286204b4715fd26c1f7 *DESCRIPTION
 0ab9635778542a5905b28d6cd7a339b0 *NAMESPACE
-458f726e6fe0e86e5f1608bc9f46ef3f *NEWS.md
-99bb41a8d47f9e5c9686f6c54cac4188 *R/dnbconv.R
-5994f87d65766a9d48f50c8558ccaed5 *R/nb_sum_exact.R
-00b27d04a6e966f7858af01fc0600e0f *R/nb_sum_moments.R
-44707a8116d7eb1fc833a718a81105d9 *R/nb_sum_saddlepoint.R
-1a264257abe9facfbcae593ccec786a6 *R/nbconv_params.R
-7afb7decba2428d80c348d90eff1620c *R/pnbconv.R
-d9ab1b6104179ee145bb0f61b9bc365b *R/qnbconv.R
-b6d18007ed2e927fab8741ec53751f80 *R/rnbconv.R
+519bf552041da368a349de3d65a3b399 *NEWS.md
+57c85610915e0b31b321dc8994112125 *R/dnbconv.R
+0006fdea71e5282265cb4b07485c6758 *R/nb_sum_exact.R
+aa85374549301c9c4328b52d726de485 *R/nb_sum_moments.R
+74b04cc73c5e146ebbdc4beddb5a5f80 *R/nb_sum_saddlepoint.R
+f1b00372d54c05d088bc8d5a69dd22a2 *R/nbconv_params.R
+fec6a7c01bddb53bf82b7f7347e8b034 *R/pnbconv.R
+33f08820f839ada95a2191e80fe37f95 *R/qnbconv.R
+5fad3df49125bee2cd05fc5b64e69c92 *R/rnbconv.R
 3fa298bd3de98df62a9d1dd55c6e732d *README.md
-baf003c448485c4e5358634418af1913 *man/dnbconv.Rd
-118ecda1aa82761fb4d4dd4680edb4b5 *man/nb_sum_exact.Rd
-b2179398035882a63d87cb97f17a62e4 *man/nb_sum_moments.Rd
-50be4396a808bf89f96351b294a0b6a2 *man/nb_sum_saddlepoint.Rd
-3cb368dc0901eae43070ca9ef63be6a2 *man/nbconv_params.Rd
-45cf1613f62b02b98495af264cc459c1 *man/pnbconv.Rd
-2e987ebbeca2edf473e1b00532dc372a *man/qnbconv.Rd
-fc820024b18e01a9f1e19b075846d749 *man/rnbconv.Rd
+487d9c688940afb10f1c12f75d2c58e2 *man/dnbconv.Rd
+7cab25be638f556ee552eb070db4104c *man/nb_sum_exact.Rd
+38e36d322325a4af4e3443b0e59eb5a8 *man/nb_sum_moments.Rd
+476e5acc33372943a377fae200020ef3 *man/nb_sum_saddlepoint.Rd
+84c6f17a0be03f165250ad014a64a56f *man/nbconv_params.Rd
+627e0172d738fa8701935c4426494598 *man/pnbconv.Rd
+532fa2a7cec0a5b43683968e397b6d52 *man/qnbconv.Rd
+4b26b1481cec0f478d46697dfbaec76c *man/rnbconv.Rd

NEWS.md

@@ -1,3 +1,14 @@
 # nbconv 0.0.0.9000
 
 * Added a `NEWS.md` file to track changes to the package.
+
+# nbconv 1.0.0
+
+* nbconv is on CRAN!
+
+# nbconv 1.0.1
+
+* Fixed bug in nb_sum_saddlepoint() that could lead to NaN errors when solving the saddlepoint equation with stats::uniroot().
+* Added parallelization in rnbconv() to speed up large-scale sampling.
+* Updated notation in nb_sum_saddlepoint() for clarity.
+* Added additional checks to ensure: 0 < ps <= 1, mus >= 0, phis > 0.

R/dnbconv.R

@@ -1,3 +1,5 @@
+#' Probability mass function
+#'
 #' Calculates the PMF for the convolution of arbitrary negative binomial random variables.
 #'
 #'@param counts The counts over which the convolution is evaluated. Should be a vector.
@@ -22,26 +24,42 @@
 dnbconv <- function(counts, mus, ps, phis, method = c("exact", "moments", "saddlepoint"),
                     n.terms = 1000, n.cores = 1, tolerance = 1e-3, normalize = TRUE){
 
-  method <- match.arg(method, c("exact", "moments", "saddlepoint"))
+  method <- match.arg( method )
 
   if( method != "exact" & method != "moments" & method != "saddlepoint"){
     stop("method must be one of 'exact', 'moments', or 'saddlepoint'.", call. = FALSE)
   }
 
   if (!missing(ps) & !missing(mus)){
-    stop("'mus' and 'ps' both specified", call. = FALSE)
+    stop("mus and ps both specified", call. = FALSE)
   }
 
   if (missing(ps) & missing(mus)){
-    stop("One of 'mus' and 'ps' must be specified", call. = FALSE)
+    stop("One of mus and ps must be specified", call. = FALSE)
+  }
+
+  if ( !missing( ps ) ){
+    if ( any( ps <= 0 ) | any( ps > 1 ) ){
+      stop("ps must be 0 < ps <= 1", call. = FALSE)
+    }
+  }
+
+  if ( any( phis <= 0 ) ){
+    stop("phis must be > 0.", call. = FALSE)
   }
 
-  if( method == "exact" ){
+  if ( !missing( mus ) ){
+    if ( any( mus < 0 ) ){
+      stop("mus must be > 0.", call. = FALSE)
+    }
+  }
+
+  if ( method == "exact" ){
     if (missing(ps) & !missing(mus)){
       ps <- phis / ( phis + mus )
     }
     if (length(ps) != length(phis)){
-      stop("'ps' and 'phis' must have the same length.", call. = FALSE)
+      stop("ps and phis must have the same length.", call. = FALSE)
     }
   }
 
@@ -50,7 +68,7 @@ dnbconv <- function(counts, mus, ps, phis, method = c("exact", "moments", "saddl
       mus <- phis*(1 - ps)/ps
     }
     if (length(mus) != length(phis)){
-      stop("'mus' and 'phis' must have the same length.", call. = FALSE)
+      stop("mus and phis must have the same length.", call. = FALSE)
     }
   }
 
@@ -61,3 +79,4 @@ dnbconv <- function(counts, mus, ps, phis, method = c("exact", "moments", "saddl
 
   return( pmf )
 }
+

R/nb_sum_exact.R

@@ -1,3 +1,5 @@
+#' Furman's PMF
+#'
 #' Implements Furman's exact PMF for the evaluation of the sum of arbitrary NB random variables. Called by other functions. Not intended to be run alone.
 #'
 #'@param phis Vector of individual dispersion parameters. Equivalent to 'size' in dnbinom.
@@ -18,6 +20,7 @@
 #'
 nb_sum_exact <- function(phis, ps, n.terms = 1000, counts, n.cores = 1, tolerance = 1e-3){
   # Implements the PMF described in https://ssrn.com/abstract=1650365
+  # Adapted from https://github.com/slundberg/NBConvolution.jl/blob/master/src/furman.jl
 
   qs <- 1 - ps
   pmax <- max(ps)

R/nb_sum_moments.R

@@ -1,3 +1,5 @@
+#' Method of moments
+#'
 #' Implements the method of moments approximation for the sum of arbitrary NB random variables. Called by other functions. Not intended to be run alone.
 #'
 #'@param mus Vector of individual mean values.

R/nb_sum_saddlepoint.R

@@ -1,3 +1,5 @@
+#' Saddlepoint approximation
+#'
 #' Implements the saddlepoint approximation for the sum of arbitrary NB random variables. Called by other functions. Not intended to be run alone.
 #'
 #' Inspired by https://www.martinmodrak.cz/2019/06/20/approximate-densities-for-sums-of-variables-negative-binomials-and-saddlepoint/
@@ -22,13 +24,13 @@
 nb_sum_saddlepoint <- function(mus, phis, counts, normalize = TRUE, n.cores = 1){
 
   # Inspired by https://www.martinmodrak.cz/2019/06/20/approximate-densities-for-sums-of-variables-negative-binomials-and-saddlepoint/
-  
+
   saddlepoint_calc <- function(mus, phis, counts){
 
     K  <-  function(t) { sum( phis * ( log( phis ) - log( phis + mus * ( 1 - exp( t ) ) ) ) ) }
-    dK <-  function(t) { logSumExp( log( phis ) + log( mus ) + t - log( phis + mus - mus * exp( t ) ) ) }
-    ddK <- function(t) { logSumExp( log( phis ) + log( mus ) + log( phis + mus ) + t - 2 * log( phis + mus - mus * exp( t ) ) ) }
-    pmf_eq <- function(t, x) { -0.5 * ( log( 2 * pi ) + ddK( t ) ) + K( t ) - t * x }
+    ldK <-  function(t) { logSumExp( log( phis ) + log( mus ) + t - log( phis + mus - mus * exp( t ) ) ) }
+    lddK <- function(t) { logSumExp( log( phis ) + log( mus ) + log( phis + mus ) + t - 2 * log( phis + mus - mus * exp( t ) ) ) }
+    pmf_eq <- function(t, x) { -0.5 * ( log( 2 * pi ) + lddK( t ) ) + K( t ) - t * x }
 
     if ( min( counts ) == 0 ){
       pmf0 <- prod( pnbinom( 0, size = phis, mu = mus ) )
@@ -38,9 +40,10 @@ nb_sum_saddlepoint <- function(mus, phis, counts, normalize = TRUE, n.cores = 1)
 
     pmf <- sapply(X = counts,
                   FUN = function(x) {
-                    t <- uniroot(function(t) { dK(t) - log(x) },
+                    t <- uniroot(function(t) { ldK(t) - log(x) },
                                  lower = -1e2,
                                  upper = min( log( phis / mus + 1 ) ),
+                                 f.upper = min( log( phis / mus + 1 ) ),
                                  extendInt = "yes",
                                  tol = sqrt( .Machine$double.eps ) )$root
                     pmf <- pmf_eq(t, x)
@@ -69,7 +72,7 @@ nb_sum_saddlepoint <- function(mus, phis, counts, normalize = TRUE, n.cores = 1)
                          FUN = function(y) {
                            pmf <- saddlepoint_calc(mus = mus,
                                                    phis = phis,
-                                                   counts = y)
+                                                   counts = y )
                            return(pmf) },
                          mc.cores = n.cores)
 

R/nbconv_params.R

@@ -1,3 +1,5 @@
+#' Summary statistics
+#'
 #' Calculates distribution parameters for the convolution of arbitrary negative binomial random variables.
 #'
 #'@param mus Vector of individual mean values
@@ -16,11 +18,27 @@
 nbconv_params <- function(mus, phis, ps){
 
   if (!missing(ps) & !missing(mus)){
-    stop("'mus' and 'ps' both specified", call. = FALSE)
+    stop("mus and ps both specified", call. = FALSE)
   }
 
   if (missing(ps) & missing(mus)){
-    stop("One of 'mus' and 'ps' must be specified", call. = FALSE)
+    stop("One of mus and ps must be specified", call. = FALSE)
+  }
+
+  if ( !missing( ps ) ){
+    if ( any( ps <= 0 ) | any( ps > 1 ) ){
+      stop("ps must be 0 < ps <= 1", call. = FALSE)
+    }
+  }
+
+  if ( any( phis <= 0 ) ){
+    stop("phis must be > 0.", call. = FALSE)
+  }
+
+  if ( !missing( mus ) ){
+    if ( any( mus < 0 ) ){
+      stop("mus must be > 0.", call. = FALSE)
+    }
   }
 
   if (missing(mus) & !missing(ps)){
@@ -45,6 +63,7 @@ nbconv_params <- function(mus, phis, ps){
 
   mean <- k1
   sigma2 <- k2
+  sigma <- sqrt ( sigma2 )
   skewness <- k3 / k2^(3/2)
   ekurtosis <- k4 / k2^2
 
@@ -53,7 +72,7 @@ nbconv_params <- function(mus, phis, ps){
 
   K.mean <- ( mean * pmax / qmax ) - sum( phis )
 
-  params <- c( mean = mean, sigma2 = sigma2, skewness = skewness, ekurtosis = ekurtosis, K.mean = K.mean)
+  params <- c( mean = mean, sigma2 = sigma2, skewness = skewness, ex.kurtosis = ekurtosis, K.mean = K.mean )
 
   return( params )
 }

R/pnbconv.R

@@ -1,3 +1,5 @@
+#' Cumulative distribution function
+#'
 #' Calculates the CDF for the convolution of arbitrary negative binomial random variables.
 #'
 #'@param quants Vector of quantiles.
@@ -24,18 +26,34 @@ pnbconv <- function(quants, mus, ps, phis, method = c("exact", "moments", "saddl
 
   counts <- 0:max( quants )
 
-  method <- match.arg(method, c("exact", "moments", "saddlepoint"))
+  method <- match.arg( method )
 
   if( method != "exact" & method != "moments" & method != "saddlepoint"){
     stop("method must be one of 'exact', 'moments', or 'saddlepoint'.", call. = FALSE)
   }
 
   if (!missing(ps) & !missing(mus)){
-    stop("'mus' and 'ps' both specified", call. = FALSE)
+    stop("mus and ps both specified", call. = FALSE)
   }
 
   if (missing(ps) & missing(mus)){
-    stop("One of 'mus' and 'ps' must be specified", call. = FALSE)
+    stop("One of mus and ps must be specified", call. = FALSE)
+  }
+
+  if ( !missing( ps ) ){
+    if ( any( ps <= 0 ) | any( ps > 1 ) ){
+      stop("ps must be 0 < ps <= 1", call. = FALSE)
+    }
+  }
+
+  if ( any( phis <= 0 ) ){
+    stop("phis must be > 0.", call. = FALSE)
+  }
+
+  if ( !missing( mus ) ){
+    if ( any( mus < 0 ) ){
+      stop("mus must be > 0.", call. = FALSE)
+    }
   }
 
   if( method == "exact" ){

R/qnbconv.R

@@ -1,3 +1,5 @@
+#' Quantile function
+#'
 #' Calculates the quantile function for the convolution of arbitrary negative binomial random variables.
 #'
 #'@param probs Vector of target (cumulative) probabilities.
@@ -27,20 +29,36 @@ qnbconv <- function(probs, counts, mus, ps, phis, method = c("exact", "moments",
 
     counts <- 0:max( counts )
 
-    method <- match.arg(method, c("exact", "moments", "saddlepoint"))
+    method <- match.arg( method )
 
     if( method != "exact" & method != "moments" & method != "saddlepoint"){
       stop("method must be one of 'exact', 'moments', or 'saddlepoint'.", call. = FALSE)
     }
 
     if (!missing(ps) & !missing(mus)){
-      stop("'mus' and 'ps' both specified", call. = FALSE)
+      stop("mus and ps both specified", call. = FALSE)
     }
 
     if (missing(ps) & missing(mus)){
-      stop("One of 'mus' and 'ps' must be specified", call. = FALSE)
+      stop("One of mus and ps must be specified", call. = FALSE)
+    }
+
+    if ( !missing( ps ) ){
+      if ( any( ps <= 0 ) | any( ps > 1 ) ){
+        stop("ps must be 0 < ps <= 1", call. = FALSE)
+      }
     }
 
+    if ( any( phis <= 0 ) ){
+      stop("phis must be > 0.", call. = FALSE)
+    }
+
+    if ( !missing( mus ) ){
+      if ( any( mus < 0 ) ){
+        stop("mus must be > 0.", call. = FALSE)
+        }
+      }
+
     if( method == "exact" ){
       if (missing(ps) & !missing(mus)){
         ps <- phis/(phis + mus)