Add underdispersed Poisson functions

NAMESPACE

@@ -45,6 +45,7 @@ export(dev_pois)
 export(dev_pois_zi)
 export(dev_skewnorm)
 export(dev_student)
+export(dev_upois)
 export(dskewnorm)
 export(exp10)
 export(exp2)
@@ -72,6 +73,7 @@ export(log_lik_pois)
 export(log_lik_pois_zi)
 export(log_lik_skewnorm)
 export(log_lik_student)
+export(log_lik_upois)
 export(log_odds)
 export(log_odds_ratio)
 export(log_odds_ratio2)
@@ -109,6 +111,7 @@ export(ran_pois)
 export(ran_pois_zi)
 export(ran_skewnorm)
 export(ran_student)
+export(ran_upois)
 export(rbern)
 export(res_bern)
 export(res_beta_binom)
@@ -123,6 +126,7 @@ export(res_pois)
 export(res_pois_zi)
 export(res_skewnorm)
 export(res_student)
+export(res_upois)
 export(rskewnorm)
 export(sextreme)
 export(skewness)

R/dev.R

@@ -287,3 +287,28 @@ dev_student <- function(x, mean = 0, sd = 1, theta = 0, res = FALSE) {
   dev_res(x, mean, dev)
 }
 
+#' Underdispersed Poisson Deviances
+#'
+#' @inheritParams params
+#' @param x A non-negative whole numeric vector of values.
+#'
+#' @return An numeric vector of the corresponding deviances or deviance residuals.
+#' @family dev_dist
+#' @export
+#'
+#' @examples
+#' dev_upois(c(1,3.5,4), 3, 2)
+dev_upois <- function(x, lambda = 1, theta = 0, res = FALSE) {
+  sat <- log_lik_upois(x = x, lambda = pmax(x - (theta / (1 + theta)), 0), theta = theta)
+  ll <- log_lik_upois(x = x, lambda = lambda, theta = theta)
+  dev <- sat - ll
+  dev <- 2 * dev
+  # Some cases where dev < 0, but all are very small diff (largest 1e-3)
+  neg <- dev < 0
+  dev[neg] <- 0
+  use_pois <- !is.na(theta) & theta == 0
+  dev_pois <- dev_pois(x = x, lambda = lambda, res = FALSE)
+  dev[use_pois] <- dev_pois[use_pois]
+  if(vld_false(res)) return(dev)
+  dev_res(x, lambda + (theta / (1 + theta)), dev)
+}

R/log-lik.R

@@ -240,3 +240,48 @@ log_lik_student <- function(x, mean = 0, sd = 1, theta = 0) {
   lstudent[use_norm] <- lnorm[use_norm]
   lstudent
 }
+
+#' Underdispersed Poisson Log-Likelihood
+#'
+#' @inheritParams params
+#' @param x A non-negative whole numeric vector of values.
+#'
+#' @return An numeric vector of the corresponding log-likelihoods.
+#' @family log_lik_dist
+#' @export
+#'
+#' @examples
+#' log_lik_upois(c(0, 1, 2), 1, 0)
+log_lik_upois <- function(x, lambda = 1, theta = 0) {
+  chk_gte(lambda)
+  chk_gte(theta)
+
+  if (length(theta) == 1) {
+    theta <- rep(theta, length(x))
+  }
+
+  if (length(lambda) == 1) {
+    lambda <- rep(lambda, length(x))
+  }
+
+  log_lik <-
+    (-lambda + (x - 1) * log(lambda) + log(lambda + theta * x)) -
+    (log(1 + theta) + lfactorial(x))
+
+  x_lambda_zero <- x == 0 & lambda == 0 & !is.na(x) & !is.na(lambda)
+  if (any(x_lambda_zero)) {
+    log_lik[x_lambda_zero] <- log(1 / (1 + theta[x_lambda_zero]))
+  }
+
+  x_one_lambda_zero <- x == 1 & lambda == 0 & !is.na(x) & !is.na(lambda)
+  if (any(x_one_lambda_zero)) {
+    log_lik[x_one_lambda_zero] <- log(theta[x_one_lambda_zero] / (1 + theta[x_one_lambda_zero]))
+  }
+
+  zero <- theta == 0 & !is.na(theta) & !is.na(lambda) & !is.na(x)
+  if (any(zero)) {
+    log_lik[zero] <- log_lik_pois(x[zero], lambda[zero])
+  }
+
+  log_lik
+}

R/params.R

@@ -18,7 +18,7 @@
 #' @param size A non-negative whole numeric vector of the number of trials.
 #' @param prob A numeric vector of values between 0 and 1 of the probability of success.
 #' @param lambda A non-negative numeric vector of means.
-#' @param theta A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson and beta-binomial).
+#' @param theta A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson, beta-binomial, and underdispersed Poisson).
 #' @param shape A non-negative numeric vector of shape.
 #' @param rate A non-negative numeric vector of rate.
 #' @param type A string of the residual type. 'raw' for raw residuals 'dev' for deviance residuals and 'data' for the data.

R/ran.R

@@ -200,3 +200,56 @@ ran_student <- function(n = 1, mean = 0, sd = 1, theta = 0) {
   r <- x * sd + mean
   r
 }
+
+# Cumulative distribution function for underdispersed poisson distribution
+pupois <- function(q, lambda, theta) {
+  sapply(q, function(x) {sum(exp(log_lik_upois(0:x, lambda, theta)))})
+}
+
+#' Underdispersed Poisson Random Samples
+#'
+#' @inheritParams params
+#' @return A numeric vector of the random samples.
+#' @family ran_dist
+#' @export
+#'
+#' @examples
+#' ran_upois(n = 10, lambda = 1, theta = 1)
+ran_upois <- function(n = 1, lambda = 1, theta = 0) {
+  chk_whole_number(n)
+  chk_gte(n)
+  chk_gte(lambda)
+  chk_gte(theta)
+
+  use_pois <- all(theta == 0 & !is.na(theta))
+  if (use_pois) {
+    stats::rpois(n, lambda = lambda)
+  }
+
+  fun <- function(n, lambda, theta) {
+    u = stats::runif(n)
+    cmf <- pupois(0:max((ceiling(lambda) * 3), 50), lambda = lambda, theta = theta)
+    last <- min(which(abs(cmf - 1) < 1e-8))
+    first = min(which(cmf > 0))
+    cmf = unique(c(0, cmf[first:last]))
+    cmf_tbl = table(cut(u, breaks = cmf, include.lowest = TRUE))
+    x = rep(1:length(cmf_tbl), as.numeric(cmf_tbl)) - 1 + ifelse(first > 1, first, 0)
+    x <- as.integer(x)
+    samp <- sample(1:n, size = n)
+    x[samp]
+  }
+
+  lambda_vec <- length(lambda) > 1L & length(unique(lambda)) != 1
+  theta_vec <- length(theta) > 1L & length(unique(theta)) != 1
+  if (lambda_vec | theta_vec) {
+    n_rep <- rep(1, n)
+    mapply(
+      fun,
+      n = n_rep,
+      lambda = lambda,
+      theta = theta
+    )
+  } else {
+    fun(n, lambda[1], theta[1])
+  }
+}

R/res.R

@@ -345,3 +345,33 @@ res_student <- function(x, mean = 0, sd = 1, theta = 0, type = "dev", simulate =
          dev = dev_student(x, mean = mean, sd = sd, theta = theta, res = TRUE),
          chk_subset(x, c("data", "raw", "dev", "standardized")))
 }
+
+#' Underdispersed Poisson Residuals
+#'
+#' @inheritParams params
+#' @param x A non-negative whole numeric vector of values.
+#'
+#' @return An numeric vector of the corresponding residuals.
+#' @family res_dist
+#' @export
+#'
+#' @examples
+#' res_upois(c(0, 1, 2), 1, 1)
+res_upois <- function(x, lambda = 1, theta = 0, type = "dev", simulate = FALSE) {
+  chk_string(type)
+  if(!vld_false(simulate)) {
+    x <- ran_upois(length(x), lambda = lambda, theta = theta)
+  }
+  if (length(lambda) == 1) {
+    lambda <- rep(lambda, length(x))
+  }
+  if (length(theta) == 1) {
+    theta <- rep(theta, length(x))
+  }
+  switch(type,
+         data = x,
+         raw = x - (lambda + (theta / (1 + theta))),
+         standardized = (x - (lambda + (theta / (1 + theta)))) / sqrt(lambda + (theta / (1 + theta)^2)),
+         dev = dev_upois(x, lambda = lambda, theta = theta, res = TRUE),
+         chk_subset(x, c("data", "raw", "dev", "standardized")))
+}

_pkgdown.yml

@@ -99,6 +99,7 @@ reference:
   - '`dev_pois_zi`'
   - '`dev_skewnorm`'
   - '`dev_student`'
+  - '`dev_upois`'
 - title: Residual
   desc: Raw and Deviance Residuals
   contents:
@@ -115,6 +116,7 @@ reference:
   - '`res_pois_zi`'
   - '`res_skewnorm`'
   - '`res_student`'
+  - '`res_upois`'
 - title: Log-Likelihood
   desc: Log-likelihood functions
   contents:
@@ -131,6 +133,7 @@ reference:
   - '`log_lik_pois_zi`'
   - '`log_lik_skewnorm`'
   - '`log_lik_student`'
+  - '`log_lik_upois`'
 - title: Random
   desc: Random sample functions
   contents:
@@ -147,6 +150,7 @@ reference:
   - '`ran_pois_zi`'
   - '`ran_skewnorm`'
   - '`ran_student`'
+  - '`ran_upois`'
 - title: Bernoulli
   desc: Bernoulli distribution functions
   contents: