diff --git a/NAMESPACE b/NAMESPACE index 680601a1..bf0d7dc7 100644 --- a/NAMESPACE +++ b/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) diff --git a/NEWS.md b/NEWS.md index f34816fb..52a4a3aa 100644 --- a/NEWS.md +++ b/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()'. diff --git a/R/dev.R b/R/dev.R index 4c2d17b1..1a5f52e5 100644 --- a/R/dev.R +++ b/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 diff --git a/R/log-lik.R b/R/log-lik.R index 9b2cfe72..c9563dfa 100644 --- a/R/log-lik.R +++ b/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 diff --git a/R/ran.R b/R/ran.R index 8c659585..d2f9f4b7 100644 --- a/R/ran.R +++ b/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 diff --git a/R/res.R b/R/res.R index a1f1496a..f6119522 100644 --- a/R/res.R +++ b/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 diff --git a/man/dev_gamma_pois_zi.Rd b/man/dev_gamma_pois_zi.Rd index 790d700f..c67a9d00 100644 --- a/man/dev_gamma_pois_zi.Rd +++ b/man/dev_gamma_pois_zi.Rd @@ -29,4 +29,8 @@ Zero-Inflated Gamma-Poisson Deviances \examples{ dev_gamma_pois_zi(c(1, 3, 4), 3, 2) } +\seealso{ +Other dev_dist # make live when complete: +\code{\link{dev_gamma_pois_zt}()} +} \concept{dev_dist # make live when complete} diff --git a/man/dev_gamma_pois_zt.Rd b/man/dev_gamma_pois_zt.Rd new file mode 100644 index 00000000..df87402b --- /dev/null +++ b/man/dev_gamma_pois_zt.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dev.R +\name{dev_gamma_pois_zt} +\alias{dev_gamma_pois_zt} +\title{Zero-Truncated Gamma-Poisson Deviances} +\usage{ +dev_gamma_pois_zt(x, lambda = 1, theta = 0, res = FALSE) +} +\arguments{ +\item{x}{A whole numeric vector of values greater than or equal to 1.} + +\item{lambda}{A non-negative numeric vector of means.} + +\item{theta}{A non-negative numeric vector of the dispersion for the mixture +models (student, gamma-Poisson and beta-binomial).} + +\item{res}{A flag specifying whether to return the deviance residual as +opposed to the deviance.} +} +\value{ +An numeric vector of the corresponding deviances or deviance residuals. +} +\description{ +Saturated log-likelihood is found numerically per observation since the +truncated likelihood has no closed-form maximum. +} +\examples{ +dev_gamma_pois_zt(c(1, 3, 4), 3, 2) +} +\seealso{ +Other dev_dist # make live when complete: +\code{\link{dev_gamma_pois_zi}()} +} +\concept{dev_dist # make live when complete} diff --git a/man/log_lik_bern.Rd b/man/log_lik_bern.Rd index cefc8642..440e1440 100644 --- a/man/log_lik_bern.Rd +++ b/man/log_lik_bern.Rd @@ -30,6 +30,7 @@ Other log_lik_dist: \code{\link{log_lik_gamma}()}, \code{\link{log_lik_gamma_pois}()}, \code{\link{log_lik_gamma_pois_zi}()}, +\code{\link{log_lik_gamma_pois_zt}()}, \code{\link{log_lik_lnorm}()}, \code{\link{log_lik_neg_binom}()}, \code{\link{log_lik_norm}()}, diff --git a/man/log_lik_beta.Rd b/man/log_lik_beta.Rd index 565d05dc..7be6aa73 100644 --- a/man/log_lik_beta.Rd +++ b/man/log_lik_beta.Rd @@ -31,6 +31,7 @@ Other log_lik_dist: \code{\link{log_lik_gamma}()}, \code{\link{log_lik_gamma_pois}()}, \code{\link{log_lik_gamma_pois_zi}()}, +\code{\link{log_lik_gamma_pois_zt}()}, \code{\link{log_lik_lnorm}()}, \code{\link{log_lik_neg_binom}()}, \code{\link{log_lik_norm}()}, diff --git a/man/log_lik_beta_binom.Rd b/man/log_lik_beta_binom.Rd index 4fc365e5..3e2e91ae 100644 --- a/man/log_lik_beta_binom.Rd +++ b/man/log_lik_beta_binom.Rd @@ -45,6 +45,7 @@ Other log_lik_dist: \code{\link{log_lik_gamma}()}, \code{\link{log_lik_gamma_pois}()}, \code{\link{log_lik_gamma_pois_zi}()}, +\code{\link{log_lik_gamma_pois_zt}()}, \code{\link{log_lik_lnorm}()}, \code{\link{log_lik_neg_binom}()}, \code{\link{log_lik_norm}()}, diff --git a/man/log_lik_binom.Rd b/man/log_lik_binom.Rd index 53369af5..863dd5b6 100644 --- a/man/log_lik_binom.Rd +++ b/man/log_lik_binom.Rd @@ -32,6 +32,7 @@ Other log_lik_dist: \code{\link{log_lik_gamma}()}, \code{\link{log_lik_gamma_pois}()}, \code{\link{log_lik_gamma_pois_zi}()}, +\code{\link{log_lik_gamma_pois_zt}()}, \code{\link{log_lik_lnorm}()}, \code{\link{log_lik_neg_binom}()}, \code{\link{log_lik_norm}()}, diff --git a/man/log_lik_exp.Rd b/man/log_lik_exp.Rd index 2dbbd29f..95187fa0 100644 --- a/man/log_lik_exp.Rd +++ b/man/log_lik_exp.Rd @@ -29,6 +29,7 @@ Other log_lik_dist: \code{\link{log_lik_gamma}()}, \code{\link{log_lik_gamma_pois}()}, \code{\link{log_lik_gamma_pois_zi}()}, +\code{\link{log_lik_gamma_pois_zt}()}, \code{\link{log_lik_lnorm}()}, \code{\link{log_lik_neg_binom}()}, \code{\link{log_lik_norm}()}, diff --git a/man/log_lik_gamma.Rd b/man/log_lik_gamma.Rd index ae073723..7e1c3b47 100644 --- a/man/log_lik_gamma.Rd +++ b/man/log_lik_gamma.Rd @@ -31,6 +31,7 @@ Other log_lik_dist: \code{\link{log_lik_exp}()}, \code{\link{log_lik_gamma_pois}()}, \code{\link{log_lik_gamma_pois_zi}()}, +\code{\link{log_lik_gamma_pois_zt}()}, \code{\link{log_lik_lnorm}()}, \code{\link{log_lik_neg_binom}()}, \code{\link{log_lik_norm}()}, diff --git a/man/log_lik_gamma_pois.Rd b/man/log_lik_gamma_pois.Rd index 9041c389..17059b56 100644 --- a/man/log_lik_gamma_pois.Rd +++ b/man/log_lik_gamma_pois.Rd @@ -32,6 +32,7 @@ Other log_lik_dist: \code{\link{log_lik_exp}()}, \code{\link{log_lik_gamma}()}, \code{\link{log_lik_gamma_pois_zi}()}, +\code{\link{log_lik_gamma_pois_zt}()}, \code{\link{log_lik_lnorm}()}, \code{\link{log_lik_neg_binom}()}, \code{\link{log_lik_norm}()}, diff --git a/man/log_lik_gamma_pois_zi.Rd b/man/log_lik_gamma_pois_zi.Rd index fed48ce0..145f3772 100644 --- a/man/log_lik_gamma_pois_zi.Rd +++ b/man/log_lik_gamma_pois_zi.Rd @@ -35,6 +35,7 @@ Other log_lik_dist: \code{\link{log_lik_exp}()}, \code{\link{log_lik_gamma}()}, \code{\link{log_lik_gamma_pois}()}, +\code{\link{log_lik_gamma_pois_zt}()}, \code{\link{log_lik_lnorm}()}, \code{\link{log_lik_neg_binom}()}, \code{\link{log_lik_norm}()}, diff --git a/man/log_lik_gamma_pois_zt.Rd b/man/log_lik_gamma_pois_zt.Rd new file mode 100644 index 00000000..cb9b6c2e --- /dev/null +++ b/man/log_lik_gamma_pois_zt.Rd @@ -0,0 +1,45 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/log-lik.R +\name{log_lik_gamma_pois_zt} +\alias{log_lik_gamma_pois_zt} +\title{Zero-Truncated Gamma-Poisson Log-Likelihood} +\usage{ +log_lik_gamma_pois_zt(x, lambda = 1, theta = 0) +} +\arguments{ +\item{x}{A whole numeric vector of values greater than or equal to 1.} + +\item{lambda}{A non-negative numeric vector of means.} + +\item{theta}{A non-negative numeric vector of the dispersion for the mixture +models (student, gamma-Poisson and beta-binomial).} +} +\value{ +An numeric vector of the corresponding log-likelihoods. +} +\description{ +Zero-Truncated Gamma-Poisson Log-Likelihood +} +\examples{ +log_lik_gamma_pois_zt(c(1, 3, 4), 3, 1) +} +\seealso{ +Other log_lik_dist: +\code{\link{log_lik_bern}()}, +\code{\link{log_lik_beta}()}, +\code{\link{log_lik_beta_binom}()}, +\code{\link{log_lik_binom}()}, +\code{\link{log_lik_exp}()}, +\code{\link{log_lik_gamma}()}, +\code{\link{log_lik_gamma_pois}()}, +\code{\link{log_lik_gamma_pois_zi}()}, +\code{\link{log_lik_lnorm}()}, +\code{\link{log_lik_neg_binom}()}, +\code{\link{log_lik_norm}()}, +\code{\link{log_lik_pois}()}, +\code{\link{log_lik_pois_zi}()}, +\code{\link{log_lik_skewnorm}()}, +\code{\link{log_lik_student}()}, +\code{\link{log_lik_unif}()} +} +\concept{log_lik_dist} diff --git a/man/log_lik_lnorm.Rd b/man/log_lik_lnorm.Rd index 7d7d3e02..515f5724 100644 --- a/man/log_lik_lnorm.Rd +++ b/man/log_lik_lnorm.Rd @@ -33,6 +33,7 @@ Other log_lik_dist: \code{\link{log_lik_gamma}()}, \code{\link{log_lik_gamma_pois}()}, \code{\link{log_lik_gamma_pois_zi}()}, +\code{\link{log_lik_gamma_pois_zt}()}, \code{\link{log_lik_neg_binom}()}, \code{\link{log_lik_norm}()}, \code{\link{log_lik_pois}()}, diff --git a/man/log_lik_neg_binom.Rd b/man/log_lik_neg_binom.Rd index 8d47937e..4b7f7adc 100644 --- a/man/log_lik_neg_binom.Rd +++ b/man/log_lik_neg_binom.Rd @@ -33,6 +33,7 @@ Other log_lik_dist: \code{\link{log_lik_gamma}()}, \code{\link{log_lik_gamma_pois}()}, \code{\link{log_lik_gamma_pois_zi}()}, +\code{\link{log_lik_gamma_pois_zt}()}, \code{\link{log_lik_lnorm}()}, \code{\link{log_lik_norm}()}, \code{\link{log_lik_pois}()}, diff --git a/man/log_lik_norm.Rd b/man/log_lik_norm.Rd index 8d1b3c34..c48f3ca2 100644 --- a/man/log_lik_norm.Rd +++ b/man/log_lik_norm.Rd @@ -32,6 +32,7 @@ Other log_lik_dist: \code{\link{log_lik_gamma}()}, \code{\link{log_lik_gamma_pois}()}, \code{\link{log_lik_gamma_pois_zi}()}, +\code{\link{log_lik_gamma_pois_zt}()}, \code{\link{log_lik_lnorm}()}, \code{\link{log_lik_neg_binom}()}, \code{\link{log_lik_pois}()}, diff --git a/man/log_lik_pois.Rd b/man/log_lik_pois.Rd index 26e0de54..d56f178e 100644 --- a/man/log_lik_pois.Rd +++ b/man/log_lik_pois.Rd @@ -30,6 +30,7 @@ Other log_lik_dist: \code{\link{log_lik_gamma}()}, \code{\link{log_lik_gamma_pois}()}, \code{\link{log_lik_gamma_pois_zi}()}, +\code{\link{log_lik_gamma_pois_zt}()}, \code{\link{log_lik_lnorm}()}, \code{\link{log_lik_neg_binom}()}, \code{\link{log_lik_norm}()}, diff --git a/man/log_lik_pois_zi.Rd b/man/log_lik_pois_zi.Rd index 8f378a9f..10e45c06 100644 --- a/man/log_lik_pois_zi.Rd +++ b/man/log_lik_pois_zi.Rd @@ -33,6 +33,7 @@ Other log_lik_dist: \code{\link{log_lik_gamma}()}, \code{\link{log_lik_gamma_pois}()}, \code{\link{log_lik_gamma_pois_zi}()}, +\code{\link{log_lik_gamma_pois_zt}()}, \code{\link{log_lik_lnorm}()}, \code{\link{log_lik_neg_binom}()}, \code{\link{log_lik_norm}()}, diff --git a/man/log_lik_skewnorm.Rd b/man/log_lik_skewnorm.Rd index 93e60f2a..3a1e25e1 100644 --- a/man/log_lik_skewnorm.Rd +++ b/man/log_lik_skewnorm.Rd @@ -38,6 +38,7 @@ Other log_lik_dist: \code{\link{log_lik_gamma}()}, \code{\link{log_lik_gamma_pois}()}, \code{\link{log_lik_gamma_pois_zi}()}, +\code{\link{log_lik_gamma_pois_zt}()}, \code{\link{log_lik_lnorm}()}, \code{\link{log_lik_neg_binom}()}, \code{\link{log_lik_norm}()}, diff --git a/man/log_lik_student.Rd b/man/log_lik_student.Rd index d392447c..1bb261ca 100644 --- a/man/log_lik_student.Rd +++ b/man/log_lik_student.Rd @@ -35,6 +35,7 @@ Other log_lik_dist: \code{\link{log_lik_gamma}()}, \code{\link{log_lik_gamma_pois}()}, \code{\link{log_lik_gamma_pois_zi}()}, +\code{\link{log_lik_gamma_pois_zt}()}, \code{\link{log_lik_lnorm}()}, \code{\link{log_lik_neg_binom}()}, \code{\link{log_lik_norm}()}, diff --git a/man/log_lik_unif.Rd b/man/log_lik_unif.Rd index a2996430..27da9ef6 100644 --- a/man/log_lik_unif.Rd +++ b/man/log_lik_unif.Rd @@ -32,6 +32,7 @@ Other log_lik_dist: \code{\link{log_lik_gamma}()}, \code{\link{log_lik_gamma_pois}()}, \code{\link{log_lik_gamma_pois_zi}()}, +\code{\link{log_lik_gamma_pois_zt}()}, \code{\link{log_lik_lnorm}()}, \code{\link{log_lik_neg_binom}()}, \code{\link{log_lik_norm}()}, diff --git a/man/ran_bern.Rd b/man/ran_bern.Rd index 2ccc2555..3be44ce3 100644 --- a/man/ran_bern.Rd +++ b/man/ran_bern.Rd @@ -29,6 +29,7 @@ Other ran_dist: \code{\link{ran_gamma}()}, \code{\link{ran_gamma_pois}()}, \code{\link{ran_gamma_pois_zi}()}, +\code{\link{ran_gamma_pois_zt}()}, \code{\link{ran_lnorm}()}, \code{\link{ran_neg_binom}()}, \code{\link{ran_norm}()}, diff --git a/man/ran_beta_binom.Rd b/man/ran_beta_binom.Rd index bfca892e..e80beb49 100644 --- a/man/ran_beta_binom.Rd +++ b/man/ran_beta_binom.Rd @@ -42,6 +42,7 @@ Other ran_dist: \code{\link{ran_gamma}()}, \code{\link{ran_gamma_pois}()}, \code{\link{ran_gamma_pois_zi}()}, +\code{\link{ran_gamma_pois_zt}()}, \code{\link{ran_lnorm}()}, \code{\link{ran_neg_binom}()}, \code{\link{ran_norm}()}, diff --git a/man/ran_binom.Rd b/man/ran_binom.Rd index 7084843e..78823fc0 100644 --- a/man/ran_binom.Rd +++ b/man/ran_binom.Rd @@ -31,6 +31,7 @@ Other ran_dist: \code{\link{ran_gamma}()}, \code{\link{ran_gamma_pois}()}, \code{\link{ran_gamma_pois_zi}()}, +\code{\link{ran_gamma_pois_zt}()}, \code{\link{ran_lnorm}()}, \code{\link{ran_neg_binom}()}, \code{\link{ran_norm}()}, diff --git a/man/ran_gamma.Rd b/man/ran_gamma.Rd index 923a88fe..edba7369 100644 --- a/man/ran_gamma.Rd +++ b/man/ran_gamma.Rd @@ -30,6 +30,7 @@ Other ran_dist: \code{\link{ran_binom}()}, \code{\link{ran_gamma_pois}()}, \code{\link{ran_gamma_pois_zi}()}, +\code{\link{ran_gamma_pois_zt}()}, \code{\link{ran_lnorm}()}, \code{\link{ran_neg_binom}()}, \code{\link{ran_norm}()}, diff --git a/man/ran_gamma_pois.Rd b/man/ran_gamma_pois.Rd index 50e0c60b..57d39ef2 100644 --- a/man/ran_gamma_pois.Rd +++ b/man/ran_gamma_pois.Rd @@ -31,6 +31,7 @@ Other ran_dist: \code{\link{ran_binom}()}, \code{\link{ran_gamma}()}, \code{\link{ran_gamma_pois_zi}()}, +\code{\link{ran_gamma_pois_zt}()}, \code{\link{ran_lnorm}()}, \code{\link{ran_neg_binom}()}, \code{\link{ran_norm}()}, diff --git a/man/ran_gamma_pois_zi.Rd b/man/ran_gamma_pois_zi.Rd index 5a70d40c..bb574037 100644 --- a/man/ran_gamma_pois_zi.Rd +++ b/man/ran_gamma_pois_zi.Rd @@ -34,6 +34,7 @@ Other ran_dist: \code{\link{ran_binom}()}, \code{\link{ran_gamma}()}, \code{\link{ran_gamma_pois}()}, +\code{\link{ran_gamma_pois_zt}()}, \code{\link{ran_lnorm}()}, \code{\link{ran_neg_binom}()}, \code{\link{ran_norm}()}, diff --git a/man/ran_gamma_pois_zt.Rd b/man/ran_gamma_pois_zt.Rd new file mode 100644 index 00000000..fb2a1064 --- /dev/null +++ b/man/ran_gamma_pois_zt.Rd @@ -0,0 +1,43 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ran.R +\name{ran_gamma_pois_zt} +\alias{ran_gamma_pois_zt} +\title{Zero-Truncated Gamma-Poisson Random Samples} +\usage{ +ran_gamma_pois_zt(n = 1, lambda = 1, theta = 0) +} +\arguments{ +\item{n}{A non-negative whole number of the number of random samples to +generate.} + +\item{lambda}{A non-negative numeric vector of means.} + +\item{theta}{A non-negative numeric vector of the dispersion for the mixture +models (student, gamma-Poisson and beta-binomial).} +} +\value{ +A numeric vector of the random samples. +} +\description{ +Drawn via rejection sampling so all returned values are at least 1. +} +\examples{ +ran_gamma_pois_zt(10, lambda = 3, theta = 1) +} +\seealso{ +Other ran_dist: +\code{\link{ran_bern}()}, +\code{\link{ran_beta_binom}()}, +\code{\link{ran_binom}()}, +\code{\link{ran_gamma}()}, +\code{\link{ran_gamma_pois}()}, +\code{\link{ran_gamma_pois_zi}()}, +\code{\link{ran_lnorm}()}, +\code{\link{ran_neg_binom}()}, +\code{\link{ran_norm}()}, +\code{\link{ran_pois}()}, +\code{\link{ran_pois_zi}()}, +\code{\link{ran_skewnorm}()}, +\code{\link{ran_student}()} +} +\concept{ran_dist} diff --git a/man/ran_lnorm.Rd b/man/ran_lnorm.Rd index 03ac9612..d1cef892 100644 --- a/man/ran_lnorm.Rd +++ b/man/ran_lnorm.Rd @@ -32,6 +32,7 @@ Other ran_dist: \code{\link{ran_gamma}()}, \code{\link{ran_gamma_pois}()}, \code{\link{ran_gamma_pois_zi}()}, +\code{\link{ran_gamma_pois_zt}()}, \code{\link{ran_neg_binom}()}, \code{\link{ran_norm}()}, \code{\link{ran_pois}()}, diff --git a/man/ran_neg_binom.Rd b/man/ran_neg_binom.Rd index 19567ec2..062c1750 100644 --- a/man/ran_neg_binom.Rd +++ b/man/ran_neg_binom.Rd @@ -32,6 +32,7 @@ Other ran_dist: \code{\link{ran_gamma}()}, \code{\link{ran_gamma_pois}()}, \code{\link{ran_gamma_pois_zi}()}, +\code{\link{ran_gamma_pois_zt}()}, \code{\link{ran_lnorm}()}, \code{\link{ran_norm}()}, \code{\link{ran_pois}()}, diff --git a/man/ran_norm.Rd b/man/ran_norm.Rd index 4bc7ce72..8ab610d8 100644 --- a/man/ran_norm.Rd +++ b/man/ran_norm.Rd @@ -31,6 +31,7 @@ Other ran_dist: \code{\link{ran_gamma}()}, \code{\link{ran_gamma_pois}()}, \code{\link{ran_gamma_pois_zi}()}, +\code{\link{ran_gamma_pois_zt}()}, \code{\link{ran_lnorm}()}, \code{\link{ran_neg_binom}()}, \code{\link{ran_pois}()}, diff --git a/man/ran_pois.Rd b/man/ran_pois.Rd index 92c70b62..716bbb10 100644 --- a/man/ran_pois.Rd +++ b/man/ran_pois.Rd @@ -29,6 +29,7 @@ Other ran_dist: \code{\link{ran_gamma}()}, \code{\link{ran_gamma_pois}()}, \code{\link{ran_gamma_pois_zi}()}, +\code{\link{ran_gamma_pois_zt}()}, \code{\link{ran_lnorm}()}, \code{\link{ran_neg_binom}()}, \code{\link{ran_norm}()}, diff --git a/man/ran_pois_zi.Rd b/man/ran_pois_zi.Rd index cb707984..de49883a 100644 --- a/man/ran_pois_zi.Rd +++ b/man/ran_pois_zi.Rd @@ -32,6 +32,7 @@ Other ran_dist: \code{\link{ran_gamma}()}, \code{\link{ran_gamma_pois}()}, \code{\link{ran_gamma_pois_zi}()}, +\code{\link{ran_gamma_pois_zt}()}, \code{\link{ran_lnorm}()}, \code{\link{ran_neg_binom}()}, \code{\link{ran_norm}()}, diff --git a/man/ran_skewnorm.Rd b/man/ran_skewnorm.Rd index 805c4d11..7d616456 100644 --- a/man/ran_skewnorm.Rd +++ b/man/ran_skewnorm.Rd @@ -37,6 +37,7 @@ Other ran_dist: \code{\link{ran_gamma}()}, \code{\link{ran_gamma_pois}()}, \code{\link{ran_gamma_pois_zi}()}, +\code{\link{ran_gamma_pois_zt}()}, \code{\link{ran_lnorm}()}, \code{\link{ran_neg_binom}()}, \code{\link{ran_norm}()}, diff --git a/man/ran_student.Rd b/man/ran_student.Rd index ffbdb0bc..4d5e8f66 100644 --- a/man/ran_student.Rd +++ b/man/ran_student.Rd @@ -34,6 +34,7 @@ Other ran_dist: \code{\link{ran_gamma}()}, \code{\link{ran_gamma_pois}()}, \code{\link{ran_gamma_pois_zi}()}, +\code{\link{ran_gamma_pois_zt}()}, \code{\link{ran_lnorm}()}, \code{\link{ran_neg_binom}()}, \code{\link{ran_norm}()}, diff --git a/man/res_bern.Rd b/man/res_bern.Rd index 13cf16c4..bff0965d 100644 --- a/man/res_bern.Rd +++ b/man/res_bern.Rd @@ -33,6 +33,7 @@ Other res_dist: \code{\link{res_gamma}()}, \code{\link{res_gamma_pois}()}, \code{\link{res_gamma_pois_zi}()}, +\code{\link{res_gamma_pois_zt}()}, \code{\link{res_lnorm}()}, \code{\link{res_neg_binom}()}, \code{\link{res_norm}()}, diff --git a/man/res_beta_binom.Rd b/man/res_beta_binom.Rd index b7fd0675..621cdf05 100644 --- a/man/res_beta_binom.Rd +++ b/man/res_beta_binom.Rd @@ -53,6 +53,7 @@ Other res_dist: \code{\link{res_gamma}()}, \code{\link{res_gamma_pois}()}, \code{\link{res_gamma_pois_zi}()}, +\code{\link{res_gamma_pois_zt}()}, \code{\link{res_lnorm}()}, \code{\link{res_neg_binom}()}, \code{\link{res_norm}()}, diff --git a/man/res_binom.Rd b/man/res_binom.Rd index ad2fc8e1..ec20d03c 100644 --- a/man/res_binom.Rd +++ b/man/res_binom.Rd @@ -35,6 +35,7 @@ Other res_dist: \code{\link{res_gamma}()}, \code{\link{res_gamma_pois}()}, \code{\link{res_gamma_pois_zi}()}, +\code{\link{res_gamma_pois_zt}()}, \code{\link{res_lnorm}()}, \code{\link{res_neg_binom}()}, \code{\link{res_norm}()}, diff --git a/man/res_gamma.Rd b/man/res_gamma.Rd index 6b98bcfd..9f1bd8cb 100644 --- a/man/res_gamma.Rd +++ b/man/res_gamma.Rd @@ -34,6 +34,7 @@ Other res_dist: \code{\link{res_binom}()}, \code{\link{res_gamma_pois}()}, \code{\link{res_gamma_pois_zi}()}, +\code{\link{res_gamma_pois_zt}()}, \code{\link{res_lnorm}()}, \code{\link{res_neg_binom}()}, \code{\link{res_norm}()}, diff --git a/man/res_gamma_pois.Rd b/man/res_gamma_pois.Rd index c8d276cd..f36e9b03 100644 --- a/man/res_gamma_pois.Rd +++ b/man/res_gamma_pois.Rd @@ -35,6 +35,7 @@ Other res_dist: \code{\link{res_binom}()}, \code{\link{res_gamma}()}, \code{\link{res_gamma_pois_zi}()}, +\code{\link{res_gamma_pois_zt}()}, \code{\link{res_lnorm}()}, \code{\link{res_neg_binom}()}, \code{\link{res_norm}()}, diff --git a/man/res_gamma_pois_zi.Rd b/man/res_gamma_pois_zi.Rd index 69aa04c0..21982b30 100644 --- a/man/res_gamma_pois_zi.Rd +++ b/man/res_gamma_pois_zi.Rd @@ -44,6 +44,7 @@ Other res_dist: \code{\link{res_binom}()}, \code{\link{res_gamma}()}, \code{\link{res_gamma_pois}()}, +\code{\link{res_gamma_pois_zt}()}, \code{\link{res_lnorm}()}, \code{\link{res_neg_binom}()}, \code{\link{res_norm}()}, diff --git a/man/res_gamma_pois_zt.Rd b/man/res_gamma_pois_zt.Rd new file mode 100644 index 00000000..ec3c4b8b --- /dev/null +++ b/man/res_gamma_pois_zt.Rd @@ -0,0 +1,47 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/res.R +\name{res_gamma_pois_zt} +\alias{res_gamma_pois_zt} +\title{Zero-Truncated Gamma-Poisson Residuals} +\usage{ +res_gamma_pois_zt(x, lambda = 1, theta = 0, type = "dev", simulate = FALSE) +} +\arguments{ +\item{x}{A whole numeric vector of values greater than or equal to 1.} + +\item{lambda}{A non-negative numeric vector of means.} + +\item{theta}{A non-negative numeric vector of the dispersion for the mixture +models (student, gamma-Poisson and beta-binomial).} + +\item{type}{A string of the residual type. 'raw' for raw residuals 'dev' for +deviance residuals and 'data' for the data.} + +\item{simulate}{A flag specifying whether to simulate residuals.} +} +\value{ +An numeric vector of the corresponding residuals. +} +\description{ +Centred on the truncated mean and standardised by the truncated variance. +} +\examples{ +res_gamma_pois_zt(c(1, 2, 5), 2, 1) +} +\seealso{ +Other res_dist: +\code{\link{res_bern}()}, +\code{\link{res_beta_binom}()}, +\code{\link{res_binom}()}, +\code{\link{res_gamma}()}, +\code{\link{res_gamma_pois}()}, +\code{\link{res_gamma_pois_zi}()}, +\code{\link{res_lnorm}()}, +\code{\link{res_neg_binom}()}, +\code{\link{res_norm}()}, +\code{\link{res_pois}()}, +\code{\link{res_pois_zi}()}, +\code{\link{res_skewnorm}()}, +\code{\link{res_student}()} +} +\concept{res_dist} diff --git a/man/res_lnorm.Rd b/man/res_lnorm.Rd index 948bbed0..bb4d559b 100644 --- a/man/res_lnorm.Rd +++ b/man/res_lnorm.Rd @@ -36,6 +36,7 @@ Other res_dist: \code{\link{res_gamma}()}, \code{\link{res_gamma_pois}()}, \code{\link{res_gamma_pois_zi}()}, +\code{\link{res_gamma_pois_zt}()}, \code{\link{res_neg_binom}()}, \code{\link{res_norm}()}, \code{\link{res_pois}()}, diff --git a/man/res_neg_binom.Rd b/man/res_neg_binom.Rd index c8680a53..6c2e8412 100644 --- a/man/res_neg_binom.Rd +++ b/man/res_neg_binom.Rd @@ -36,6 +36,7 @@ Other res_dist: \code{\link{res_gamma}()}, \code{\link{res_gamma_pois}()}, \code{\link{res_gamma_pois_zi}()}, +\code{\link{res_gamma_pois_zt}()}, \code{\link{res_lnorm}()}, \code{\link{res_norm}()}, \code{\link{res_pois}()}, diff --git a/man/res_norm.Rd b/man/res_norm.Rd index 0e964dd3..1ac62259 100644 --- a/man/res_norm.Rd +++ b/man/res_norm.Rd @@ -35,6 +35,7 @@ Other res_dist: \code{\link{res_gamma}()}, \code{\link{res_gamma_pois}()}, \code{\link{res_gamma_pois_zi}()}, +\code{\link{res_gamma_pois_zt}()}, \code{\link{res_lnorm}()}, \code{\link{res_neg_binom}()}, \code{\link{res_pois}()}, diff --git a/man/res_pois.Rd b/man/res_pois.Rd index 68e0ff4f..7d490785 100644 --- a/man/res_pois.Rd +++ b/man/res_pois.Rd @@ -33,6 +33,7 @@ Other res_dist: \code{\link{res_gamma}()}, \code{\link{res_gamma_pois}()}, \code{\link{res_gamma_pois_zi}()}, +\code{\link{res_gamma_pois_zt}()}, \code{\link{res_lnorm}()}, \code{\link{res_neg_binom}()}, \code{\link{res_norm}()}, diff --git a/man/res_pois_zi.Rd b/man/res_pois_zi.Rd index 43b3df93..fb42947d 100644 --- a/man/res_pois_zi.Rd +++ b/man/res_pois_zi.Rd @@ -35,6 +35,7 @@ Other res_dist: \code{\link{res_gamma}()}, \code{\link{res_gamma_pois}()}, \code{\link{res_gamma_pois_zi}()}, +\code{\link{res_gamma_pois_zt}()}, \code{\link{res_lnorm}()}, \code{\link{res_neg_binom}()}, \code{\link{res_norm}()}, diff --git a/man/res_skewnorm.Rd b/man/res_skewnorm.Rd index 0d2c53c5..765e35ec 100644 --- a/man/res_skewnorm.Rd +++ b/man/res_skewnorm.Rd @@ -39,6 +39,7 @@ Other res_dist: \code{\link{res_gamma}()}, \code{\link{res_gamma_pois}()}, \code{\link{res_gamma_pois_zi}()}, +\code{\link{res_gamma_pois_zt}()}, \code{\link{res_lnorm}()}, \code{\link{res_neg_binom}()}, \code{\link{res_norm}()}, diff --git a/man/res_student.Rd b/man/res_student.Rd index 52d7d783..d30a4c71 100644 --- a/man/res_student.Rd +++ b/man/res_student.Rd @@ -38,6 +38,7 @@ Other res_dist: \code{\link{res_gamma}()}, \code{\link{res_gamma_pois}()}, \code{\link{res_gamma_pois_zi}()}, +\code{\link{res_gamma_pois_zt}()}, \code{\link{res_lnorm}()}, \code{\link{res_neg_binom}()}, \code{\link{res_norm}()}, diff --git a/tests/testthat/test-dev.R b/tests/testthat/test-dev.R index 6bd2a5e0..2a620eff 100644 --- a/tests/testthat/test-dev.R +++ b/tests/testthat/test-dev.R @@ -316,6 +316,38 @@ test_that("gamma_pois_zi vectorized", { expect_equal(dev_gamma_pois_zi(0:3, 3:0, 0:3, seq(0, 1, length.out = 4)), c(6, 1.0464962875291, 2.41568518074605, Inf)) }) +test_that("gamma_pois_zt missing values", { + expect_identical(dev_gamma_pois_zt(logical(0), integer(0), numeric(0)), numeric(0)) + expect_identical(dev_gamma_pois_zt(NA, 1, 1), NA_real_) + expect_identical(dev_gamma_pois_zt(1, NA, 1), NA_real_) + expect_identical(dev_gamma_pois_zt(1, 1, NA), NA_real_) +}) + +test_that("gamma_pois_zt known values", { + expect_equal(dev_gamma_pois_zt(c(1, 3, 4), 3, 2), + c(2.69098973, 0.19520565, 0.03152929), + tolerance = 1e-6) + expect_equal(dev_gamma_pois_zt(3, 3), 0.009522421, tolerance = 1e-6) +}) + +test_that("gamma_pois_zt deviance is non-negative", { + expect_true(all(dev_gamma_pois_zt(1:20, 5, 1) >= 0)) + expect_true(all(dev_gamma_pois_zt(1:20, 10, 0.5) >= 0)) + expect_true(all(dev_gamma_pois_zt(1:20, 1, 2) >= 0)) +}) + +test_that("gamma_pois_zt deviance residuals", { + # Residuals signed by sign(x - trunc_mean) and scaled by sqrt(dev). + expect_equal(dev_gamma_pois_zt(c(1, 2, 5), 2, 1, res = TRUE), + c(-1.4823038, -0.4853515, 0.6610002), + tolerance = 1e-6) + # Sign follows sign(x - trunc_mean). For lambda = 2, theta = 1, + # trunc_mean = 2 / (1 - 1/3) = 3. + trunc_mean_at_2_1 <- 2 / (1 - dnbinom(0, mu = 2, size = 1)) + expect_equal(sign(dev_gamma_pois_zt(c(1, 5), 2, 1, res = TRUE)), + sign(c(1, 5) - trunc_mean_at_2_1)) +}) + test_that("gamma_pois_zi vectorized missing values", { expect_equal(dev_gamma_pois_zi(c(NA, 1), 0:1, 0:1, 0:1), c(NA, Inf)) expect_equal(dev_gamma_pois_zi(c(0, NA), 0:1, 0:1, 0:1), c(0, NA)) diff --git a/tests/testthat/test-log-lik.R b/tests/testthat/test-log-lik.R index aca1621c..fd354232 100644 --- a/tests/testthat/test-log-lik.R +++ b/tests/testthat/test-log-lik.R @@ -113,6 +113,35 @@ test_that("gamma_pois_zi vectorized", { expect_equal(log_lik_gamma_pois_zi(0:3, 3:0, 0:3, seq(0, 1, length.out = 4)), c(-3, -1.90954250488444, -3.43967790223022, -Inf)) }) +test_that("gamma_pois_zt missing values", { + expect_identical(log_lik_gamma_pois_zt(numeric(0), numeric(0), numeric(0)), numeric(0)) + expect_identical(log_lik_gamma_pois_zt(NA, 1, 1), NA_real_) + expect_identical(log_lik_gamma_pois_zt(1, NA, 1), NA_real_) + expect_identical(log_lik_gamma_pois_zt(1, 1, NA), NA_real_) +}) + +test_that("gamma_pois_zt known values", { + expect_equal(log_lik_gamma_pois_zt(1, 2), -1.1614393615712) + expect_equal(log_lik_gamma_pois_zt(2, 2), -1.1614393615712) + expect_equal(log_lik_gamma_pois_zt(1, 2, 0.5), -1.09861228866811) + expect_equal(log_lik_gamma_pois_zt(c(1, 3, 4), 3, 1), + c(-1.38629436111989, -1.96165850602345, -2.24934057847523)) +}) + +test_that("gamma_pois_zt vectorized", { + expect_equal(log_lik_gamma_pois_zt(1:3, 2, 0), + c(-1.1614393615712, -1.1614393615712, -1.56690446967936)) +}) + +test_that("gamma_pois_zt matches truncated NB conditional", { + # log L_zt(x; mu, theta) = log NB(x; mu, theta) - log(1 - P(0; mu, theta)) + expect_equal( + log_lik_gamma_pois_zt(2, 3, 1), + dnbinom(2, mu = 3, size = 1, log = TRUE) - + log1p(-dnbinom(0, mu = 3, size = 1)) + ) +}) + test_that("gamma missing values", { expect_identical(log_lik_gamma(NA), NA_real_) expect_identical(log_lik_gamma(1, NA), NA_real_) diff --git a/tests/testthat/test-ran.R b/tests/testthat/test-ran.R index 0c375d75..185c0b6e 100644 --- a/tests/testthat/test-ran.R +++ b/tests/testthat/test-ran.R @@ -146,6 +146,39 @@ test_that("ran_gamma_pois_zi", { expect_identical(ran_gamma_pois_zi(1, c(0.1, 100), theta = 1, prob = 0.5), 0L) }) +test_that("ran_gamma_pois_zt", { + expect_error(ran_gamma_pois_zt(NA_integer_)) + expect_error(ran_gamma_pois_zt(integer(0))) + expect_identical(ran_gamma_pois_zt(0L), integer(0)) + set.seed(101) + expect_identical(ran_gamma_pois_zt(), 1L) + expect_identical(ran_gamma_pois_zt(2), c(3L, 2L)) + expect_identical(ran_gamma_pois_zt(2, 10), c(7L, 7L)) + expect_identical(ran_gamma_pois_zt(2, c(0.1, 100)), c(1L, 95L)) + set.seed(101) + expect_identical(ran_gamma_pois_zt(2, theta = 10), c(2L, 1L)) + expect_identical(ran_gamma_pois_zt(2, 10, theta = 10), c(7L, 13L)) + expect_identical(ran_gamma_pois_zt(2, c(0.1, 100), theta = 1), c(1L, 12L)) +}) + +test_that("ran_gamma_pois_zt always >= 1", { + set.seed(42) + # Even with very small lambda where most untruncated draws would be 0, + # rejection sampling guarantees all returned values are >= 1. + y <- ran_gamma_pois_zt(2000, lambda = 0.5, theta = 1) + expect_true(all(y >= 1)) + expect_length(y, 2000) +}) + +test_that("ran_gamma_pois_zt mean approximates truncated mean", { + set.seed(42) + lambda <- 5 + theta <- 0.5 + y <- ran_gamma_pois_zt(20000, lambda = lambda, theta = theta) + expected <- lambda / (1 - dnbinom(0, mu = lambda, size = 1 / theta)) + expect_equal(mean(y), expected, tolerance = 0.05) +}) + test_that("ran_gamma", { expect_error(ran_gamma(NA_integer_)) expect_error(ran_gamma(integer(0))) diff --git a/tests/testthat/test-res.R b/tests/testthat/test-res.R index 0bbf882c..32ad2464 100644 --- a/tests/testthat/test-res.R +++ b/tests/testthat/test-res.R @@ -380,6 +380,52 @@ test_that("res_gamma_pois_zi", { expect_equal(sd(res), 0.860406154019735) }) +test_that("res_gamma_pois_zt", { + expect_identical(res_gamma_pois_zt(integer(0), integer(0), integer(0)), numeric(0)) + expect_identical(res_gamma_pois_zt(NA, 1, 1), NA_real_) + expect_identical(res_gamma_pois_zt(1, NA, 1), NA_real_) + expect_identical(res_gamma_pois_zt(1, 1, NA), NA_real_) + expect_error(res_gamma_pois_zt(1, 1, 1, type = "unknown")) + + # type = "data" + expect_identical(res_gamma_pois_zt(1:3, 2, 1, type = "data"), 1:3) + + # type = "raw" — residual is x - truncated mean + # truncated mean for lambda = 2, theta = 1 is 2 / (1 - 1/3) = 3 + expect_equal(res_gamma_pois_zt(c(1, 3, 5), 2, 1, type = "raw"), c(-2, 0, 2)) + # When theta -> 0, truncated mean approaches lambda for moderate lambda + expect_equal( + res_gamma_pois_zt(c(1, 10, 20), 10, 0, type = "raw"), + c(1, 10, 20) - 10 / (1 - exp(-10)) + ) + + # type = "standardized" — residual is (x - trunc_mean) / sqrt(trunc_var) + # For lambda = 3, theta = 1: trunc_mean = 4, trunc_var = 12 + expect_equal( + res_gamma_pois_zt(c(1, 3, 4), 3, 1, type = "standardized"), + (c(1, 3, 4) - 4) / sqrt(12) + ) + + # type = "dev" — saturated MLE found numerically per observation, + # signed by sign(x - trunc_mean) + expect_equal( + res_gamma_pois_zt(c(1, 2, 5), 2, 1), + c(-1.4823038, -0.4853515, 0.6610002), + tolerance = 1e-6 + ) + + # simulate = TRUE replaces x with random draws from ZT distribution + set.seed(101) + expect_true(all(res_gamma_pois_zt(1:5, 2, 1, simulate = TRUE, type = "data") >= 1)) + + # Standardised residuals on simulated data approximate N(0, 1) + set.seed(101) + res <- res_gamma_pois_zt(rep(2, 10000), 2, 1, simulate = TRUE, + type = "standardized") + expect_equal(mean(res), 0, tolerance = 0.05) + expect_equal(sd(res), 1, tolerance = 0.05) +}) + test_that("res_gamma", { expect_identical(res_gamma(integer(0), integer(0), integer(0)), numeric(0)) expect_identical(res_gamma(1, 1, 1), 0)