Add zero-truncated gamma-Poisson functions

NAMESPACE

@@ -38,6 +38,7 @@ export(dev_binom)
 export(dev_gamma)
 export(dev_gamma_pois)
 export(dev_gamma_pois_zi)
+export(dev_gamma_pois_zt)
 export(dev_lnorm)
 export(dev_neg_binom)
 export(dev_norm)
@@ -67,6 +68,7 @@ export(log_lik_exp)
 export(log_lik_gamma)
 export(log_lik_gamma_pois)
 export(log_lik_gamma_pois_zi)
+export(log_lik_gamma_pois_zt)
 export(log_lik_lnorm)
 export(log_lik_neg_binom)
 export(log_lik_norm)
@@ -106,6 +108,7 @@ export(ran_binom)
 export(ran_gamma)
 export(ran_gamma_pois)
 export(ran_gamma_pois_zi)
+export(ran_gamma_pois_zt)
 export(ran_lnorm)
 export(ran_neg_binom)
 export(ran_norm)
@@ -120,6 +123,7 @@ export(res_binom)
 export(res_gamma)
 export(res_gamma_pois)
 export(res_gamma_pois_zi)
+export(res_gamma_pois_zt)
 export(res_lnorm)
 export(res_neg_binom)
 export(res_norm)

NEWS.md

@@ -2,6 +2,7 @@
 
 # extras 0.8.0.9001
 
+- Added `log_lik_gamma_pois_zt()`, `dev_gamma_pois_zt()`, `ran_gamma_pois_zt()`, and `res_gamma_pois_zt()` for the zero-truncated gamma-Poisson (negative binomial) distribution.
 - Added `log_lik_exp()`.
 - Added `log_lik_beta()`.
 - Added 'log_lik_unif()'.

R/dev.R

@@ -244,6 +244,62 @@ dev_gamma_pois_zi <- function(x, lambda = 1, theta = 0, prob = 0, res = FALSE) {
   dev_res(x, lambda * (1 - prob), dev)
 }
 
+#' Zero-Truncated Gamma-Poisson Deviances
+#'
+#' Saturated log-likelihood is found numerically per observation since the
+#' truncated likelihood has no closed-form maximum.
+#'
+#' @inheritParams params
+#' @param x A whole numeric vector of values greater than or equal to 1.
+#'
+#' @return An numeric vector of the corresponding deviances or deviance residuals.
+#' @family dev_dist # make live when complete
+#' @export
+#'
+#' @examples
+#' dev_gamma_pois_zt(c(1, 3, 4), 3, 2)
+dev_gamma_pois_zt <- function(x, lambda = 1, theta = 0, res = FALSE) {
+  force(x)
+  args_not_na <- !is.na(x + lambda + theta)
+  if (length(theta) == 1) {
+    theta_vec <- theta
+  } else {
+    theta_vec <- theta[args_not_na]
+  }
+  opt_lambda <- rep(NA_real_, length(args_not_na))
+  if (any(args_not_na)) {
+    upper <- max(x[args_not_na])
+    if (!is.finite(upper) || upper <= 0) upper <- 1
+    opt_lambda[args_not_na] <- parallel_optimize(
+      f = make_opt_gamma_pois_zt(x[args_not_na], theta_vec),
+      interval = c(.Machine$double.eps, upper),
+      N = sum(args_not_na)
+    )
+  }
+  dev1 <- log_lik_gamma_pois_zt(x = x, lambda = opt_lambda, theta = theta)
+  dev2 <- log_lik_gamma_pois_zt(x = x, lambda = lambda, theta = theta)
+  dev <- dev1 - dev2
+  dev[!is.na(dev) & dev < 0 & dev > -1e-7] <- 0
+  dev <- dev * 2
+  if (vld_false(res)) {
+    return(dev)
+  }
+  log_p0 <- dnbinom(0, mu = lambda, size = 1 / theta, log = TRUE)
+  trunc_mean <- lambda / -expm1(log_p0)
+  dev_res(x, trunc_mean, dev)
+}
+
+# Objective function for optimization.
+make_opt_gamma_pois_zt <- function(x, theta) {
+  force(x)
+  force(theta)
+  function(lambda) {
+    out <- -log_lik_gamma_pois_zt(x = x, lambda = lambda, theta = theta)
+    out[is.nan(out)] <- Inf
+    out
+  }
+}
+
 #' Log-Normal Deviances
 #'
 #' @inheritParams params

R/log-lik.R

@@ -205,6 +205,23 @@ log_lik_gamma_pois_zi <- function(x, lambda = 1, theta = 0, prob = 0) {
   log(lpois)
 }
 
+#' Zero-Truncated Gamma-Poisson Log-Likelihood
+#'
+#' @inheritParams params
+#' @param x A whole numeric vector of values greater than or equal to 1.
+#'
+#' @return An numeric vector of the corresponding log-likelihoods.
+#' @family log_lik_dist
+#' @export
+#'
+#' @examples
+#' log_lik_gamma_pois_zt(c(1, 3, 4), 3, 1)
+log_lik_gamma_pois_zt <- function(x, lambda = 1, theta = 0) {
+  log_lik <- dnbinom(x, mu = lambda, size = 1 / theta, log = TRUE)
+  log_p0 <- dnbinom(0, mu = lambda, size = 1 / theta, log = TRUE)
+  log_lik - log1p(-exp(log_p0))
+}
+
 #' Log-Normal Log-Likelihood
 #'
 #' @inheritParams params

R/ran.R

@@ -97,6 +97,40 @@ ran_gamma_pois_zi <- function(n = 1, lambda = 1, theta = 0, prob = 0) {
   ran_neg_binom(n = n, lambda = lambda, theta = theta) * ran_bern(n, prob = 1 - prob)
 }
 
+#' Zero-Truncated Gamma-Poisson Random Samples
+#'
+#' Drawn via rejection sampling so all returned values are at least 1.
+#'
+#' @inheritParams params
+#' @return A numeric vector of the random samples.
+#' @family ran_dist
+#' @export
+#'
+#' @examples
+#' ran_gamma_pois_zt(10, lambda = 3, theta = 1)
+ran_gamma_pois_zt <- function(n = 1, lambda = 1, theta = 0) {
+  chk_whole_number(n)
+  chk_gte(n)
+  if (n == 0) {
+    return(integer(0))
+  }
+  lambda <- rep_len(lambda, n)
+  theta <- rep_len(theta, n)
+  out <- integer(n)
+  remaining <- seq_len(n)
+  while (length(remaining) > 0) {
+    draws <- as.integer(stats::rnbinom(
+      length(remaining),
+      mu = lambda[remaining],
+      size = 1 / theta[remaining]
+    ))
+    accepted <- draws > 0
+    out[remaining[accepted]] <- draws[accepted]
+    remaining <- remaining[!accepted]
+  }
+  out
+}
+
 #' Log-Normal Random Samples
 #'
 #' @inheritParams params

R/res.R

@@ -173,6 +173,37 @@ res_gamma_pois_zi <- function(x, lambda = 1, theta = 0, prob = 0, type = "dev",
   )
 }
 
+#' Zero-Truncated Gamma-Poisson Residuals
+#'
+#' Centred on the truncated mean and standardised by the truncated variance.
+#'
+#' @inheritParams params
+#' @param x A whole numeric vector of values greater than or equal to 1.
+#'
+#' @return An numeric vector of the corresponding residuals.
+#' @family res_dist
+#' @export
+#'
+#' @examples
+#' res_gamma_pois_zt(c(1, 2, 5), 2, 1)
+res_gamma_pois_zt <- function(x, lambda = 1, theta = 0, type = "dev", simulate = FALSE) {
+  chk_string(type)
+  if (!vld_false(simulate)) {
+    x <- ran_gamma_pois_zt(length(x), lambda = lambda, theta = theta)
+  }
+  log_p0 <- dnbinom(0, mu = lambda, size = 1 / theta, log = TRUE)
+  one_minus_p0 <- -expm1(log_p0)
+  trunc_mean <- lambda / one_minus_p0
+  trunc_var <- (lambda + lambda^2 * (1 + theta)) / one_minus_p0 - trunc_mean^2
+  switch(type,
+    data = x,
+    raw = x - trunc_mean,
+    standardized = (x - trunc_mean) / sqrt(trunc_var),
+    dev = dev_gamma_pois_zt(x, lambda = lambda, theta = theta, res = TRUE),
+    chk_subset(x, c("data", "raw", "dev", "standardized"))
+  )
+}
+
 #' Log-Normal Residuals
 #'
 #' @inheritParams params