updates to BetaBinom, yaml, tests, documentation

.github/workflows/check.yaml

@@ -79,7 +79,7 @@ jobs:
           install.packages('devtools')
           install.packages('pracma')
           install.packages('numDeriv')
-          devtools::install_github("nimble-dev/nimble", ref = "ADoak", subdir = "packages/nimble") # Remove this line once AD is released in NIMBLE
+                    devtools::install_github("nimble-dev/nimble", ref = "ADoak", subdir = "packages/nimble") # Remove this line once AD is released in NIMBLE
         shell: Rscript {0}
 
       - name: Session info

NAMESPACE

@@ -1,7 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
-export(dBetaBinom)
-export(dBetaBinom_One)
+export(dBetaBinom_s)
+export(dBetaBinom_v)
 export(dCJS_ss)
 export(dCJS_sv)
 export(dCJS_vs)
@@ -52,8 +52,8 @@ export(dOcc_s)
 export(dOcc_v)
 export(nimBetaFun)
 export(nimNmixPois_logFac)
-export(rBetaBinom)
-export(rBetaBinom_One)
+export(rBetaBinom_s)
+export(rBetaBinom_v)
 export(rCJS_ss)
 export(rCJS_sv)
 export(rCJS_vs)

R/dBetaBinom.R

@@ -1,60 +1,59 @@
 # dBetaBinom
 #' A beta binomial distribution and beta function for use in \code{nimble} models
 #'
-#' \code{dBetaBinom} and \code{dBetaBinom_One} provide a beta binomial
+#' \code{dBetaBinom_v} and \code{dBetaBinom_s} provide a beta binomial
 #' distribution that can be used directly from R or in \code{nimble}
 #' models. These are also used by beta binomial variations of dNmixture distributions.
 #' \code{nimBetaFun} is the beta function.
 #'
 #' @name dBetaBinom
-#' @aliases nimBetaFun dBetaBinom dBetaBinom_One rBetaBinom rBetaBinom_One
+#' @aliases nimBetaFun dBetaBinom_v dBetaBinom_s rBetaBinom_v rBetaBinom_s
 #'
 #' @author Ben Goldstein and Perry de Valpine
 #'
 #' @param x vector of integer counts.
 #' @param N number of trials, sometimes called "size".
-#' @param shape1 shape1 parameter of the beta-binomial distribution.
-#' @param shape2 shape2 parameter of the beta-binomial distribution.
+#' @param shape1 shape1 parameter of the beta distribution.
+#' @param shape2 shape2 parameter of the beta distribution.
 #' @param log TRUE or 1 to return log probability. FALSE or 0 to return
 #'   probability.
 #' @param n number of random draws, each returning a vector of length
 #'   \code{len}. Currently only \code{n = 1} is supported, but the argument
 #'   exists for standardization of "\code{r}" functions.
-#' @param a shape1 argument of the beta function nimBetaFun.
-#' @param b shape2 argument of the beta function nimBetaFun.
+#' @param len length of \code{x}.
+#' @param a shape1 argument of the beta function.
+#' @param b shape2 argument of the beta function.
 #'
-#' @details These nimbleFunctions provide distributions that can be
-#'     used directly in R or in \code{nimble} hierarchical models (via
-#'     \code{\link[nimble]{nimbleCode}} and
-#'     \code{\link[nimble]{nimbleModel}}). They were originally written for
-#'     the beta binomial N-mixture extensions.
+#' @details These nimbleFunctions provide distributions that can be used
+#'   directly in R or in \code{nimble} hierarchical models (via
+#'   \code{\link[nimble]{nimbleCode}} and \code{\link[nimble]{nimbleModel}}).
+#'   They are used by the beta-binomial variants of the N-mixture distributions
+#'   (\code{\link{dNmixture}}).
 #'
-#' The beta binomial distribution is equivalent to a binomial distribution whose
-#' probability is itself a beta distributed random variable.
+#' The beta binomial is the marginal distribution of a  binomial distribution whose
+#' probability follows a  beta distribution.
 #'
 #' The probability mass function of the beta binomial is
 #' \code{choose(N, x) * B(x + shape1, N - x + shape2) /
 #' B(shape1, shape2)}, where \code{B(shape1, shape2)} is the beta function.
 #'
-#' The beta binomial distribution is provided in two forms. \code{dBetaBinom} and
-#' \code{rBetaBinom} are used when \code{x} is a vector (i.e. \code{length(x) > 1}),
-#' in which case the parameters \code{alpha} and \code{beta} must also be vectors.
-#' When \code{x} is scalar, \code{dBetaBinom_One} and \code{rBetaBinom_One} are
-#' used.
+#' \code{nimBetaFun(shape1, shape2)} calculates \code{B(shape1, shape2)}.
+#'
+#' The beta binomial distribution is provided in two forms. \code{dBetaBinom_v} and
+#' is when \code{shape1} and \code{shape2} are vectors.
+#' \code{dBetaBinom_s} is used when \code{shape1} and \code{shape2} are scalars.
+#' In both cases, \code{x} is a vector.
 #'
 #' @seealso For beta binomial N-mixture models, see \code{\link{dNmixture}}.
 #'   For documentation on the beta function, use \code{?stats::dbeta}
 #'
 #' @examples
-#' # Calculate a beta binomial probability
-#' dBetaBinom(x = c(4, 0, 0, 3), N = 10,
+#' # Calculate a beta binomial probability with different shape1 and shape2 for each x[i]
+#' dBetaBinom_v(x = c(4, 0, 0, 3), N = 10,
 #'   shape1 = c(0.5, 0.5, 0.3, 0.5), shape2 = c(0.2, 0.4, 1, 1.2))
-#' # Same for case with one observation
-#' dBetaBinom_One(x = 3, N = 10, shape1 = 0.5, shape2 = 0.5, log = TRUE)
-#' @export
+#' # or with constant shape1 and shape2
+#' dBetaBinom_s(x = c(4, 0, 0, 3), N = 10, shape1 = 0.5, shape2 = 0.5, log = TRUE)
 
-NULL
-# nimbleOptions(checkNimbleFunction = FALSE)
 
 ##### Beta binomial support functions #####
 #' @rdname dBetaBinom
@@ -71,18 +70,19 @@ nimBetaFun <- nimbleFunction(
 
 #' @rdname dBetaBinom
 #' @export
-dBetaBinom <- nimbleFunction(
+dBetaBinom_v <- nimbleFunction(
   run = function(x = double(1),
                  N = double(0),
                  shape1 = double(1),
                  shape2 = double(1),
+                 len = double(),
                  log = integer(0, default = 0)) {
     logprob <- 0
     lgNp1 <- lgamma(N+1)
     for (i in 1:length(x)) {
       logprob <- logprob +
         nimBetaFun(a = x[i] + shape1[i], b = N - x[i] + shape2[i], log = TRUE) -
-        nimBetaFun(a = shape1[i], b = shape2[ i], log = TRUE) +
+        nimBetaFun(a = shape1[i], b = shape2[i], log = TRUE) +
         lgNp1 - (lgamma(x[i] + 1) + lgamma(N - x[i] + 1))
     }
 
@@ -95,34 +95,38 @@ dBetaBinom <- nimbleFunction(
 
 #' @rdname dBetaBinom
 #' @export
-dBetaBinom_One <- nimbleFunction(
-  run = function(x = double(0),
+dBetaBinom_s <- nimbleFunction(
+  run = function(x = double(1),
                  N = double(0),
                  shape1 = double(0),
                  shape2 = double(0),
+                 len = double(),
                  log = integer(0, default = 0)) {
     logprob <- 0
-    logprob <- logprob +
-      nimBetaFun(a = x + shape1, b = N - x + shape2, log = TRUE) -
-      nimBetaFun(a = shape1, b = shape2, log = TRUE) +
-      lgamma(N+1) - (lgamma(x+1) + lgamma(N - x + 1))
-
+    lgNp1 <- lgamma(N+1)
+    lbs1s2 <- nimBetaFun(a = shape1, b = shape2, log = TRUE)
+    for (i in 1:length(x)) {
+      logprob <- logprob +
+        nimBetaFun(a = x[i] + shape1, b = N - x[i] + shape2, log = TRUE) -
+        lbs1s2 +
+        lgNp1 - (lgamma(x[i] + 1) + lgamma(N - x[i] + 1))
+    }
     if (log) return(logprob)
     return(exp(logprob))
     returnType(double(0))
   },
   buildDerivs = list(run=list())
 )
 
-
 #' @rdname dBetaBinom
 #' @export
 #' @importFrom stats rbeta
-rBetaBinom <- nimbleFunction(
+rBetaBinom_v <- nimbleFunction(
   run = function(n = double(0),
                  N = double(0),
                  shape1 = double(1),
-                 shape2 = double(1)) {
+                 shape2 = double(1),
+                 len = double()) {
     p <- numeric(length(shape1))
     for (i in 1:length(shape1)) {
       p[i] <- rbeta(1, shape1[i], shape2[i])
@@ -135,16 +139,17 @@ rBetaBinom <- nimbleFunction(
 #' @rdname dBetaBinom
 #' @export
 #' @importFrom stats rbeta
-rBetaBinom_One <- nimbleFunction(
+rBetaBinom_s <- nimbleFunction(
   run = function(n = double(0),
                  N = double(0),
                  shape1 = double(0),
-                 shape2 = double(0)) {
-    p <- rbeta(1, shape1, shape2)
-    x <- rbinom(1, N, p)
+                 shape2 = double(0),
+                 len = double()) {
+    p <- numeric(length=len)
+    for (i in 1:len) {
+      p[i] <- rbeta(1, shape1, shape2)
+    }
+    x <- rbinom(len, N, p)
     return(x)
-    returnType(double())
+    returnType(double(1))
   })
-
-# nimbleOptions(checkNimbleFunction = TRUE)
-

R/dCJS.R

@@ -379,7 +379,7 @@ rCJS_sv <- nimbleFunction(
                  len = integer(0, default = 0)) {
     if (n != 1) stop("rCJS only works for n = 1")
     if (len < 2)
-      stop("len must be non-negative.")
+      stop("len must be greater than 1.")
     if (length(probCapture) != len)
       stop("Length of probCapture is not the same as len.")
     ans <- numeric(length = len, init = FALSE)

R/dNmixture.R

@@ -450,7 +450,7 @@ dNmixture_BBP_v <- nimbleFunction(
     }
     Nmin <- max( max(x), Nmin ) ## set Nmin to at least the largest x
     logProb <- dNmixture_BBP_steps(x, beta-x, lambda, s, Nmin, Nmax,
-                                   dBetaBinom(x, Nmin, alpha, beta, log = TRUE))
+                                   dBetaBinom_v(x, Nmin, alpha, beta, len, log = TRUE))
     if (log) return(logProb)
     else return(exp(logProb))
     returnType(double())
@@ -484,7 +484,7 @@ dNmixture_BBP_s <- nimbleFunction(
     }
     Nmin <- max( max(x), Nmin ) ## set Nmin to at least the largest x
     logProb <- dNmixture_BBP_steps(x, beta-x, lambda, s, Nmin, Nmax,
-                                   dBetaBinom(x, Nmin, rep(alpha, len), rep(beta, len), log = TRUE))
+                                   dBetaBinom_s(x, Nmin, alpha, beta, len, log = TRUE))
     if (log) return(logProb)
     else return(exp(logProb))
     returnType(double())
@@ -537,7 +537,7 @@ dNmixture_BBNB_v <- nimbleFunction(
     }
     Nmin <- max( max(x), Nmin ) ## set Nmin to at least the largest x
     logProb <- dNmixture_BBNB_steps(x, beta-x,lambda,theta,s,Nmin,Nmax,
-                                    dBetaBinom(x, Nmin, alpha, beta, log = TRUE))
+                                    dBetaBinom_v(x, Nmin, alpha, beta, len, log = TRUE))
     if (log) return(logProb)
     else return(exp(logProb))
     returnType(double())
@@ -575,7 +575,7 @@ dNmixture_BBNB_s <- nimbleFunction(
     }
     Nmin <- max( max(x), Nmin ) ## set Nmin to at least the largest x
     logProb <- dNmixture_BBNB_steps(x, beta-x,lambda,theta,s,Nmin,Nmax,
-                                    dBetaBinom(x, Nmin, rep(alpha, len), rep(beta, len), log = TRUE))
+                                    dBetaBinom_s(x, Nmin, alpha, beta, len, log = TRUE))
     if (log) return(logProb)
     else return(exp(logProb))
     returnType(double())
@@ -703,7 +703,7 @@ rNmixture_BBP_v <- nimbleFunction(
     beta <- s - prob * s
 
     trueN <- rpois(1, lambda = lambda)
-    ans <- rBetaBinom(n = 1, N = trueN, shape1 = alpha, shape2 = beta)
+    ans <- rBetaBinom_v(n = 1, N = trueN, shape1 = alpha, shape2 = beta, len = len)
 
     return(ans)
     returnType(double(1))
@@ -727,8 +727,8 @@ rNmixture_BBP_s <- nimbleFunction(
     beta <- s - prob * s
 
     trueN <- rpois(1, lambda = lambda)
-    ans <- rBetaBinom(n = 1, N = trueN,
-                      shape1 = rep(alpha, len), shape2 = rep(beta, len))
+    ans <- rBetaBinom_s(n = 1, N = trueN,
+                      shape1 = alpha, shape2 = beta, len = len)
 
     return(ans)
     returnType(double(1))
@@ -751,9 +751,9 @@ rNmixture_BBP_oneObs <- nimbleFunction(
     beta <- s - prob * s
 
     trueN <- rpois(1, lambda = lambda)
-    ans <- rBetaBinom_One(n = 1, N = trueN, shape1 = alpha, shape2 = beta)
+    ans <- rBetaBinom_s(n = 1, N = trueN, shape1 = alpha, shape2 = beta, len = 1)
 
-    return(ans)
+    return(ans[1])
     returnType(double())
   })
 
@@ -778,7 +778,7 @@ rNmixture_BBNB_v <- nimbleFunction(
     p <- 1 / (1 + theta * lambda)
 
     trueN <- rnbinom(1, size = r, prob = p)
-    ans <- rBetaBinom(n = 1, N = trueN, shape1 = alpha, shape2 = beta)
+    ans <- rBetaBinom_v(n = 1, N = trueN, shape1 = alpha, shape2 = beta, len = len)
 
     return(ans)
     returnType(double(1))
@@ -805,8 +805,8 @@ rNmixture_BBNB_s <- nimbleFunction(
     p <- 1 / (1 + theta * lambda)
 
     trueN <- rnbinom(1, size = r, prob = p)
-    ans <- rBetaBinom(n = 1, N = trueN,
-                      shape1 = rep(alpha, len), shape2 = rep(beta, len))
+    ans <- rBetaBinom_s(n = 1, N = trueN,
+                      shape1 = alpha, shape2 = beta, len = len)
 
     return(ans)
     returnType(double(1))
@@ -832,7 +832,7 @@ rNmixture_BBNB_oneObs <- nimbleFunction(
     p <- 1 / (1 + theta * lambda)
 
     trueN <- rnbinom(1, size = r, prob = p)
-    ans <- rBetaBinom_One(n = 1, N = trueN, shape1 = alpha, shape2 = beta)
-    return(ans)
+    ans <- rBetaBinom_s(n = 1, N = trueN, shape1 = alpha, shape2 = beta, len = 1)
+    return(ans[1])
     returnType(double())
   })

R/dNmixtureAD.R

@@ -50,12 +50,11 @@
 #' In the AD system some kinds of values are "baked in" (cannot be changed) to
 #' the AD calculations from the first call, unless and until the AD calculations
 #' are reset. For all variants of the \code{dNmixtureAD} distributions, the
-#' sizes of the inputs and the data values themselves are baked in. These can be
-#' different for different iterations through a for loop (or nimble model
-#' declarations with different indices, for example), but the sizes and data
-#' values for each specific iteration will be "baked in" after the first call.
-#' \bold{In other words, it is assumed that \code{x} are data and are not going
-#' to change.}
+#' sizes of the inputs as well as \code{Nmin} and \code{Nmax} are baked in.
+#' These can be different for different iterations through a for loop (or nimble
+#' model declarations with different indices, for example), but the sizes and
+#' \code{Nmin} and \code{Nmax} values for each specific iteration will be
+#' "baked in" after the first call.
 #'
 #' @return The probability (or likelihood) or log probability of an observation
 #'   vector \code{x}.
@@ -258,7 +257,7 @@ dNmixtureAD_BBP_v <- nimbleFunction(
     }
     Nmin <- ADbreak(max( max(x), Nmin )) ## set Nmin to at least the largest x
     logProb <- dNmixture_BBP_steps(x, beta-x, lambda, s, Nmin, Nmax,
-                                   dBetaBinom(x, Nmin, alpha, beta, log = TRUE))
+                                   dBetaBinom_v(x, Nmin, alpha, beta, len = len, log = TRUE))
     if (log) return(logProb)
     else return(exp(logProb))
     returnType(double())
@@ -286,11 +285,11 @@ dNmixtureAD_BBP_s <- nimbleFunction(
     if (Nmin == -1 | Nmax == -1) {
       stop("Dynamic choice of Nmin/Nmax is not supported for beta binomial N-mixtures.")
     }
-    Clen <- 0L
-    Clen <- ADbreak(len)
+    #Clen <- 0L
+    #Clen <- ADbreak(len)
     Nmin <- ADbreak(max( max(x), Nmin )) ## set Nmin to at least the largest x
     logProb <- dNmixture_BBP_steps(x, beta-x, lambda, s, Nmin, Nmax,
-                                   dBetaBinom(x, Nmin, rep(alpha, Clen), rep(beta, Clen), log = TRUE))
+                                   dBetaBinom_v(x, Nmin, alpha, beta, len = len, log = TRUE))
     if (log) return(logProb)
     else return(exp(logProb))
     returnType(double())
@@ -342,7 +341,7 @@ dNmixtureAD_BBNB_v <- nimbleFunction(
     }
     Nmin <- ADbreak(max( max(x), Nmin )) ## set Nmin to at least the largest x
     logProb <- dNmixture_BBNB_steps(x, beta-x,lambda,theta,s,Nmin,Nmax,
-                                    dBetaBinom(x, Nmin, alpha, beta, log = TRUE))
+                                    dBetaBinom_v(x, Nmin, alpha, beta, len = len, log = TRUE))
     if (log) return(logProb)
     else return(exp(logProb))
     returnType(double())
@@ -376,11 +375,11 @@ dNmixtureAD_BBNB_s <- nimbleFunction(
     if (Nmin == -1 | Nmax == -1) {
       stop("Dynamic choice of Nmin/Nmax is not supported for beta binomial N-mixtures.")
     }
-    Clen <- 0L
-    Clen <- ADbreak(len)
+#    Clen <- 0L
+#    Clen <- ADbreak(len)
     Nmin <- ADbreak(max( max(x), Nmin )) ## set Nmin to at least the largest x
     logProb <- dNmixture_BBNB_steps(x, beta-x,lambda,theta,s,Nmin,Nmax,
-                                    dBetaBinom(x, Nmin, rep(alpha, Clen), rep(beta, Clen), log = TRUE))
+                                    dBetaBinom_s(x, Nmin, alpha, beta, len = len, log = TRUE))
     if (log) return(logProb)
     else return(exp(logProb))
     returnType(double())