proof of concept for posterior_pit

DESCRIPTION

@@ -91,14 +91,14 @@ Description: Fit Bayesian generalized (non-)linear multivariate multilevel model
     linear, robust linear, count data, survival, response times, ordinal,
     zero-inflated, hurdle, and even self-defined mixture models all in a
     multilevel context. Further modeling options include both theory-driven and
-    data-driven non-linear terms, auto-correlation structures, censoring and 
-    truncation, meta-analytic standard errors, and quite a few more. 
-    In addition, all parameters of the response distribution can be predicted 
-    in order to perform distributional regression. Prior specifications are 
+    data-driven non-linear terms, auto-correlation structures, censoring and
+    truncation, meta-analytic standard errors, and quite a few more.
+    In addition, all parameters of the response distribution can be predicted
+    in order to perform distributional regression. Prior specifications are
     flexible and explicitly encourage users to apply prior distributions that
     actually reflect their prior knowledge. Models can easily be evaluated and
-    compared using several methods assessing posterior or prior predictions. 
-    References: Bürkner (2017) <doi:10.18637/jss.v080.i01>; 
+    compared using several methods assessing posterior or prior predictions.
+    References: Bürkner (2017) <doi:10.18637/jss.v080.i01>;
     Bürkner (2018) <doi:10.32614/RJ-2018-017>;
     Bürkner (2021) <doi:10.18637/jss.v100.i05>; Carpenter et al. (2017)
     <doi:10.18637/jss.v076.i01>.

NAMESPACE

@@ -559,11 +559,13 @@ export(pzero_inflated_binomial)
 export(pzero_inflated_negbinomial)
 export(pzero_inflated_poisson)
 export(qasym_laplace)
+export(qbeta_binomial)
 export(qfrechet)
 export(qgen_extreme_value)
 export(qshifted_lnorm)
 export(qskew_normal)
 export(qstudent_t)
+export(qzero_inflated_negbinomial)
 export(ranef)
 export(rasym_laplace)
 export(rbeta_binomial)

R/distributions.R

@@ -918,6 +918,10 @@ rinv_gaussian <- function(n, mu = 1, shape = 1) {
 #'  \item{\code{phi = (1 - mu) * beta}} precision or over-dispersion, component.
 #' }
 #'
+#' The quantile function has no known closed form for this parameterization and
+#' is therefore computed numerically via inversion of the cumulative
+#' distribution function over the finite support \code{0:size}.
+#'
 #' @name BetaBinomial
 #'
 #' @inheritParams StudentT
@@ -944,6 +948,31 @@ pbeta_binomial <- function(q, size, mu, phi, lower.tail = TRUE, log.p = FALSE) {
                       lower.tail = lower.tail, log.p = log.p)
 }
 
+#' @rdname BetaBinomial
+#' @export
+qbeta_binomial <- function(p, size, mu, phi, lower.tail = TRUE, log.p = FALSE) {
+  require_package("extraDistr")
+  p <- validate_p_dist(p, lower.tail = lower.tail, log.p = log.p)
+  args <- do_call(expand, nlist(p, size, mu, phi))
+  out <- vapply(seq_along(args$p), function(i) {
+    if (!is.finite(args$p[i])) {
+      return(args$p[i])
+    }
+    alpha <- args$mu[i] * args$phi[i]
+    beta <- (1 - args$mu[i]) * args$phi[i]
+    x <- 0:args$size[i]
+    cdf <- cumsum(extraDistr::dbbinom(x, args$size[i], alpha = alpha, beta = beta))
+    ind <- which(cdf >= args$p[i])[1]
+    if (is.na(ind)) {
+      args$size[i]
+    } else {
+      x[ind]
+    }
+  }, numeric(1))
+  dim(out) <- attributes(args)$max_dim
+  out
+}
+
 #' @rdname BetaBinomial
 #' @export
 rbeta_binomial <- function(n, size, mu, phi) {
@@ -1761,6 +1790,8 @@ rxbeta <- function(...) {
 #' If \eqn{x = 0} set \eqn{f(x) = \theta + (1 - \theta) * g(0)}.
 #' Else set \eqn{f(x) = (1 - \theta) * g(x)},
 #' where \eqn{g(x)} is the density of the non-zero-inflated part.
+#' For the zero-inflated negative binomial distribution, the quantile function
+#' has no known closed form and is therefore computed numerically.
 NULL
 
 #' @rdname ZeroInflated
@@ -1778,6 +1809,14 @@ pzero_inflated_poisson <- function(q, lambda, zi, lower.tail = TRUE,
   .pzero_inflated(q, "pois", zi, pars, lower.tail, log.p)
 }
 
+#' @rdname ZeroInflated
+#' @export
+qzero_inflated_poisson <- function(p, lambda, zi, lower.tail = TRUE,
+                                   log.p = FALSE) {
+  pars <- nlist(lambda)
+  .qzero_inflated(p, "pois", zi, pars, lower.tail, log.p)
+}
+
 #' @rdname ZeroInflated
 #' @export
 dzero_inflated_negbinomial <- function(x, mu, shape, zi, log = FALSE) {
@@ -1793,6 +1832,14 @@ pzero_inflated_negbinomial <- function(q, mu, shape, zi, lower.tail = TRUE,
   .pzero_inflated(q, "nbinom", zi, pars, lower.tail, log.p)
 }
 
+#' @rdname ZeroInflated
+#' @export
+qzero_inflated_negbinomial <- function(p, mu, shape, zi, lower.tail = TRUE,
+                                       log.p = FALSE) {
+  pars <- nlist(mu, size = shape)
+  .qzero_inflated(p, "nbinom", zi, pars, lower.tail, log.p)
+}
+
 #' @rdname ZeroInflated
 #' @export
 dzero_inflated_binomial <- function(x, size, prob, zi, log = FALSE) {
@@ -1808,6 +1855,14 @@ pzero_inflated_binomial <- function(q, size, prob, zi, lower.tail = TRUE,
   .pzero_inflated(q, "binom", zi, pars, lower.tail, log.p)
 }
 
+#' @rdname ZeroInflated
+#' @export
+qzero_inflated_binomial <- function(p, size, prob, zi, lower.tail = TRUE,
+                                    log.p = FALSE) {
+  pars <- nlist(size, prob)
+  .qzero_inflated(p, "binom", zi, pars, lower.tail, log.p)
+}
+
 #' @rdname ZeroInflated
 #' @export
 dzero_inflated_beta_binomial <- function(x, size, mu, phi, zi, log = FALSE) {
@@ -1823,6 +1878,14 @@ pzero_inflated_beta_binomial <- function(q, size, mu, phi, zi,
   .pzero_inflated(q, "beta_binomial", zi, pars, lower.tail, log.p)
 }
 
+#' @rdname ZeroInflated
+#' @export
+qzero_inflated_beta_binomial <- function(p, size, mu, phi, zi,
+                                         lower.tail = TRUE, log.p = FALSE) {
+  pars <- nlist(size, mu, phi)
+  .qzero_inflated(p, "beta_binomial", zi, pars, lower.tail, log.p)
+}
+
 #' @rdname ZeroInflated
 #' @export
 dzero_inflated_beta <- function(x, shape1, shape2, zi, log = FALSE) {
@@ -1840,6 +1903,15 @@ pzero_inflated_beta <- function(q, shape1, shape2, zi, lower.tail = TRUE,
   .phurdle(q, "beta", zi, pars, lower.tail, log.p, type = "real")
 }
 
+#' @rdname ZeroInflated
+#' @export
+qzero_inflated_beta <- function(p, shape1, shape2, zi, lower.tail = TRUE,
+                                log.p = FALSE) {
+  pars <- nlist(shape1, shape2)
+  # zi_beta is technically a hurdle model
+  .qzero_inflated(p, "beta", zi, pars, lower.tail, log.p)
+}
+
 # @rdname ZeroInflated
 # @export
 dzero_inflated_asym_laplace <- function(x, mu, sigma, quantile, zi,
@@ -1920,6 +1992,29 @@ pzero_inflated_asym_laplace <- function(q, mu, sigma, quantile, zi,
   out
 }
 
+# quantile function of a zero-inflated distribution
+# @param dist name of the distribution
+# @param zi bernoulli zero-inflated parameter
+# @param pars list of parameters passed to qfun
+.qzero_inflated <- function(p, dist, zi, pars, lower.tail, log.p) {
+  stopifnot(is.list(pars))
+  dist <- as_one_character(dist)
+  lower.tail <- as_one_logical(lower.tail)
+  log.p <- as_one_logical(log.p)
+  p <- validate_p_dist(p, lower.tail = lower.tail, log.p = log.p)
+  args <- expand(dots = c(nlist(p, zi), pars))
+  p <- args$p
+  zi <- args$zi
+  pars <- args[names(pars)]
+  qfun <- paste0("q", dist)
+  p_dist <- ifelse(zi == 1, 0, (p - zi) / (1 - zi))
+  p_dist <- pmin(1, pmax(0, p_dist))
+  out <- do_call(qfun, c(list(p_dist), pars))
+  out[!is.finite(p)] <- p[!is.finite(p)]
+  dim(out) <- attributes(args)$max_dim
+  out
+}
+
 #' Hurdle Distributions
 #'
 #' Density and distribution functions for hurdle distributions.