diff --git a/NAMESPACE b/NAMESPACE
index 634a96b1..d3708fea 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/R/dev.R b/R/dev.R
index 7eed333f..22397440 100644
--- a/R/dev.R
+++ b/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)
+}
diff --git a/R/log-lik.R b/R/log-lik.R
index 487b739b..f14557fd 100644
--- a/R/log-lik.R
+++ b/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
+}
diff --git a/R/params.R b/R/params.R
index db18c7fd..7dfe79d4 100644
--- a/R/params.R
+++ b/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.
diff --git a/R/ran.R b/R/ran.R
index f43937b2..4ba107b7 100644
--- a/R/ran.R
+++ b/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])
+  }
+}
diff --git a/R/res.R b/R/res.R
index 84cd4a3f..0b3015a4 100644
--- a/R/res.R
+++ b/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")))
+}
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 0209b070..a8251bab 100644
--- a/_pkgdown.yml
+++ b/_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:
diff --git a/man/dev_bern.Rd b/man/dev_bern.Rd
index e920fb10..2839ce34 100644
--- a/man/dev_bern.Rd
+++ b/man/dev_bern.Rd
@@ -34,6 +34,7 @@ Other dev_dist:
 \code{\link{dev_pois_zi}()},
 \code{\link{dev_pois}()},
 \code{\link{dev_skewnorm}()},
-\code{\link{dev_student}()}
+\code{\link{dev_student}()},
+\code{\link{dev_upois}()}
 }
 \concept{dev_dist}
diff --git a/man/dev_beta_binom.Rd b/man/dev_beta_binom.Rd
index 4a14ed8f..482dc369 100644
--- a/man/dev_beta_binom.Rd
+++ b/man/dev_beta_binom.Rd
@@ -13,7 +13,7 @@ dev_beta_binom(x, size = 1, prob = 0.5, theta = 0, res = FALSE)
 
 \item{prob}{A numeric vector of values between 0 and 1 of the probability of success.}
 
-\item{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson and beta-binomial).}
+\item{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson, beta-binomial, and underdispersed Poisson).}
 
 \item{res}{A flag specifying whether to return the deviance residual as opposed to the deviance.}
 }
@@ -38,6 +38,7 @@ Other dev_dist:
 \code{\link{dev_pois_zi}()},
 \code{\link{dev_pois}()},
 \code{\link{dev_skewnorm}()},
-\code{\link{dev_student}()}
+\code{\link{dev_student}()},
+\code{\link{dev_upois}()}
 }
 \concept{dev_dist}
diff --git a/man/dev_binom.Rd b/man/dev_binom.Rd
index cbe89098..fdbc218f 100644
--- a/man/dev_binom.Rd
+++ b/man/dev_binom.Rd
@@ -36,6 +36,7 @@ Other dev_dist:
 \code{\link{dev_pois_zi}()},
 \code{\link{dev_pois}()},
 \code{\link{dev_skewnorm}()},
-\code{\link{dev_student}()}
+\code{\link{dev_student}()},
+\code{\link{dev_upois}()}
 }
 \concept{dev_dist}
diff --git a/man/dev_gamma.Rd b/man/dev_gamma.Rd
index cac20287..f1d9dfb2 100644
--- a/man/dev_gamma.Rd
+++ b/man/dev_gamma.Rd
@@ -36,6 +36,7 @@ Other dev_dist:
 \code{\link{dev_pois_zi}()},
 \code{\link{dev_pois}()},
 \code{\link{dev_skewnorm}()},
-\code{\link{dev_student}()}
+\code{\link{dev_student}()},
+\code{\link{dev_upois}()}
 }
 \concept{dev_dist}
diff --git a/man/dev_gamma_pois.Rd b/man/dev_gamma_pois.Rd
index 13758798..36dff2c5 100644
--- a/man/dev_gamma_pois.Rd
+++ b/man/dev_gamma_pois.Rd
@@ -11,7 +11,7 @@ dev_gamma_pois(x, lambda = 1, theta = 0, res = FALSE)
 
 \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{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson, beta-binomial, and underdispersed Poisson).}
 
 \item{res}{A flag specifying whether to return the deviance residual as opposed to the deviance.}
 }
@@ -36,6 +36,7 @@ Other dev_dist:
 \code{\link{dev_pois_zi}()},
 \code{\link{dev_pois}()},
 \code{\link{dev_skewnorm}()},
-\code{\link{dev_student}()}
+\code{\link{dev_student}()},
+\code{\link{dev_upois}()}
 }
 \concept{dev_dist}
diff --git a/man/dev_gamma_pois_zi.Rd b/man/dev_gamma_pois_zi.Rd
index d9b96d70..3b8835d2 100644
--- a/man/dev_gamma_pois_zi.Rd
+++ b/man/dev_gamma_pois_zi.Rd
@@ -11,7 +11,7 @@ dev_gamma_pois_zi(x, lambda = 1, theta = 0, prob = 0, res = FALSE)
 
 \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{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson, beta-binomial, and underdispersed Poisson).}
 
 \item{prob}{A numeric vector of values between 0 and 1 of the probability of success.}
 
diff --git a/man/dev_lnorm.Rd b/man/dev_lnorm.Rd
index 7a505195..1587ec78 100644
--- a/man/dev_lnorm.Rd
+++ b/man/dev_lnorm.Rd
@@ -36,6 +36,7 @@ Other dev_dist:
 \code{\link{dev_pois_zi}()},
 \code{\link{dev_pois}()},
 \code{\link{dev_skewnorm}()},
-\code{\link{dev_student}()}
+\code{\link{dev_student}()},
+\code{\link{dev_upois}()}
 }
 \concept{dev_dist}
diff --git a/man/dev_neg_binom.Rd b/man/dev_neg_binom.Rd
index 78248df6..117798a7 100644
--- a/man/dev_neg_binom.Rd
+++ b/man/dev_neg_binom.Rd
@@ -11,7 +11,7 @@ dev_neg_binom(x, lambda = 1, theta = 0, res = FALSE)
 
 \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{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson, beta-binomial, and underdispersed Poisson).}
 
 \item{res}{A flag specifying whether to return the deviance residual as opposed to the deviance.}
 }
@@ -36,6 +36,7 @@ Other dev_dist:
 \code{\link{dev_pois_zi}()},
 \code{\link{dev_pois}()},
 \code{\link{dev_skewnorm}()},
-\code{\link{dev_student}()}
+\code{\link{dev_student}()},
+\code{\link{dev_upois}()}
 }
 \concept{dev_dist}
diff --git a/man/dev_norm.Rd b/man/dev_norm.Rd
index af408cec..9971ad4f 100644
--- a/man/dev_norm.Rd
+++ b/man/dev_norm.Rd
@@ -36,6 +36,7 @@ Other dev_dist:
 \code{\link{dev_pois_zi}()},
 \code{\link{dev_pois}()},
 \code{\link{dev_skewnorm}()},
-\code{\link{dev_student}()}
+\code{\link{dev_student}()},
+\code{\link{dev_upois}()}
 }
 \concept{dev_dist}
diff --git a/man/dev_pois.Rd b/man/dev_pois.Rd
index f271619d..e8b057ec 100644
--- a/man/dev_pois.Rd
+++ b/man/dev_pois.Rd
@@ -34,6 +34,7 @@ Other dev_dist:
 \code{\link{dev_norm}()},
 \code{\link{dev_pois_zi}()},
 \code{\link{dev_skewnorm}()},
-\code{\link{dev_student}()}
+\code{\link{dev_student}()},
+\code{\link{dev_upois}()}
 }
 \concept{dev_dist}
diff --git a/man/dev_pois_zi.Rd b/man/dev_pois_zi.Rd
index 88cd86b8..75ff7d94 100644
--- a/man/dev_pois_zi.Rd
+++ b/man/dev_pois_zi.Rd
@@ -36,6 +36,7 @@ Other dev_dist:
 \code{\link{dev_norm}()},
 \code{\link{dev_pois}()},
 \code{\link{dev_skewnorm}()},
-\code{\link{dev_student}()}
+\code{\link{dev_student}()},
+\code{\link{dev_upois}()}
 }
 \concept{dev_dist}
diff --git a/man/dev_skewnorm.Rd b/man/dev_skewnorm.Rd
index dd71d386..a17e438b 100644
--- a/man/dev_skewnorm.Rd
+++ b/man/dev_skewnorm.Rd
@@ -40,6 +40,7 @@ Other dev_dist:
 \code{\link{dev_norm}()},
 \code{\link{dev_pois_zi}()},
 \code{\link{dev_pois}()},
-\code{\link{dev_student}()}
+\code{\link{dev_student}()},
+\code{\link{dev_upois}()}
 }
 \concept{dev_dist}
diff --git a/man/dev_student.Rd b/man/dev_student.Rd
index 3a3fb7a5..8c5a053a 100644
--- a/man/dev_student.Rd
+++ b/man/dev_student.Rd
@@ -13,7 +13,7 @@ dev_student(x, mean = 0, sd = 1, theta = 0, res = FALSE)
 
 \item{sd}{A non-negative numeric vector of the standard deviations.}
 
-\item{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson and beta-binomial).}
+\item{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson, beta-binomial, and underdispersed Poisson).}
 
 \item{res}{A flag specifying whether to return the deviance residual as opposed to the deviance.}
 }
@@ -38,6 +38,7 @@ Other dev_dist:
 \code{\link{dev_norm}()},
 \code{\link{dev_pois_zi}()},
 \code{\link{dev_pois}()},
-\code{\link{dev_skewnorm}()}
+\code{\link{dev_skewnorm}()},
+\code{\link{dev_upois}()}
 }
 \concept{dev_dist}
diff --git a/man/dev_upois.Rd b/man/dev_upois.Rd
new file mode 100644
index 00000000..dfe6d25f
--- /dev/null
+++ b/man/dev_upois.Rd
@@ -0,0 +1,42 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/dev.R
+\name{dev_upois}
+\alias{dev_upois}
+\title{Underdispersed Poisson Deviances}
+\usage{
+dev_upois(x, lambda = 1, theta = 0, res = FALSE)
+}
+\arguments{
+\item{x}{A non-negative whole numeric vector of values.}
+
+\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, beta-binomial, and underdispersed Poisson).}
+
+\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{
+Underdispersed Poisson Deviances
+}
+\examples{
+dev_upois(c(1,3.5,4), 3, 2)
+}
+\seealso{
+Other dev_dist: 
+\code{\link{dev_bern}()},
+\code{\link{dev_beta_binom}()},
+\code{\link{dev_binom}()},
+\code{\link{dev_gamma_pois}()},
+\code{\link{dev_gamma}()},
+\code{\link{dev_lnorm}()},
+\code{\link{dev_neg_binom}()},
+\code{\link{dev_norm}()},
+\code{\link{dev_pois_zi}()},
+\code{\link{dev_pois}()},
+\code{\link{dev_skewnorm}()},
+\code{\link{dev_student}()}
+}
+\concept{dev_dist}
diff --git a/man/log_lik_bern.Rd b/man/log_lik_bern.Rd
index 1b83910c..4402c1c8 100644
--- a/man/log_lik_bern.Rd
+++ b/man/log_lik_bern.Rd
@@ -33,6 +33,7 @@ Other log_lik_dist:
 \code{\link{log_lik_pois_zi}()},
 \code{\link{log_lik_pois}()},
 \code{\link{log_lik_skewnorm}()},
-\code{\link{log_lik_student}()}
+\code{\link{log_lik_student}()},
+\code{\link{log_lik_upois}()}
 }
 \concept{log_lik_dist}
diff --git a/man/log_lik_beta_binom.Rd b/man/log_lik_beta_binom.Rd
index 3ca731df..0bd3c3b4 100644
--- a/man/log_lik_beta_binom.Rd
+++ b/man/log_lik_beta_binom.Rd
@@ -13,7 +13,7 @@ log_lik_beta_binom(x, size = 1, prob = 0.5, theta = 0)
 
 \item{prob}{A numeric vector of values between 0 and 1 of the probability of success.}
 
-\item{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson and beta-binomial).}
+\item{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson, beta-binomial, and underdispersed Poisson).}
 }
 \value{
 An numeric vector of the corresponding log-likelihoods.
@@ -37,6 +37,7 @@ Other log_lik_dist:
 \code{\link{log_lik_pois_zi}()},
 \code{\link{log_lik_pois}()},
 \code{\link{log_lik_skewnorm}()},
-\code{\link{log_lik_student}()}
+\code{\link{log_lik_student}()},
+\code{\link{log_lik_upois}()}
 }
 \concept{log_lik_dist}
diff --git a/man/log_lik_binom.Rd b/man/log_lik_binom.Rd
index 260e47f9..a9813ecc 100644
--- a/man/log_lik_binom.Rd
+++ b/man/log_lik_binom.Rd
@@ -35,6 +35,7 @@ Other log_lik_dist:
 \code{\link{log_lik_pois_zi}()},
 \code{\link{log_lik_pois}()},
 \code{\link{log_lik_skewnorm}()},
-\code{\link{log_lik_student}()}
+\code{\link{log_lik_student}()},
+\code{\link{log_lik_upois}()}
 }
 \concept{log_lik_dist}
diff --git a/man/log_lik_gamma.Rd b/man/log_lik_gamma.Rd
index 7dc2924e..56d8d32b 100644
--- a/man/log_lik_gamma.Rd
+++ b/man/log_lik_gamma.Rd
@@ -35,6 +35,7 @@ Other log_lik_dist:
 \code{\link{log_lik_pois_zi}()},
 \code{\link{log_lik_pois}()},
 \code{\link{log_lik_skewnorm}()},
-\code{\link{log_lik_student}()}
+\code{\link{log_lik_student}()},
+\code{\link{log_lik_upois}()}
 }
 \concept{log_lik_dist}
diff --git a/man/log_lik_gamma_pois.Rd b/man/log_lik_gamma_pois.Rd
index 9046796e..ddf06398 100644
--- a/man/log_lik_gamma_pois.Rd
+++ b/man/log_lik_gamma_pois.Rd
@@ -11,7 +11,7 @@ log_lik_gamma_pois(x, lambda = 1, theta = 0)
 
 \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{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson, beta-binomial, and underdispersed Poisson).}
 }
 \value{
 An numeric vector of the corresponding log-likelihoods.
@@ -35,6 +35,7 @@ Other log_lik_dist:
 \code{\link{log_lik_pois_zi}()},
 \code{\link{log_lik_pois}()},
 \code{\link{log_lik_skewnorm}()},
-\code{\link{log_lik_student}()}
+\code{\link{log_lik_student}()},
+\code{\link{log_lik_upois}()}
 }
 \concept{log_lik_dist}
diff --git a/man/log_lik_gamma_pois_zi.Rd b/man/log_lik_gamma_pois_zi.Rd
index 41140716..a11686ca 100644
--- a/man/log_lik_gamma_pois_zi.Rd
+++ b/man/log_lik_gamma_pois_zi.Rd
@@ -11,7 +11,7 @@ log_lik_gamma_pois_zi(x, lambda = 1, theta = 0, prob = 0)
 
 \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{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson, beta-binomial, and underdispersed Poisson).}
 
 \item{prob}{A numeric vector of values between 0 and 1 of the probability of success.}
 }
@@ -37,6 +37,7 @@ Other log_lik_dist:
 \code{\link{log_lik_pois_zi}()},
 \code{\link{log_lik_pois}()},
 \code{\link{log_lik_skewnorm}()},
-\code{\link{log_lik_student}()}
+\code{\link{log_lik_student}()},
+\code{\link{log_lik_upois}()}
 }
 \concept{log_lik_dist}
diff --git a/man/log_lik_lnorm.Rd b/man/log_lik_lnorm.Rd
index f8fcfc2b..6f7aa152 100644
--- a/man/log_lik_lnorm.Rd
+++ b/man/log_lik_lnorm.Rd
@@ -35,6 +35,7 @@ Other log_lik_dist:
 \code{\link{log_lik_pois_zi}()},
 \code{\link{log_lik_pois}()},
 \code{\link{log_lik_skewnorm}()},
-\code{\link{log_lik_student}()}
+\code{\link{log_lik_student}()},
+\code{\link{log_lik_upois}()}
 }
 \concept{log_lik_dist}
diff --git a/man/log_lik_neg_binom.Rd b/man/log_lik_neg_binom.Rd
index 594587f1..6212690b 100644
--- a/man/log_lik_neg_binom.Rd
+++ b/man/log_lik_neg_binom.Rd
@@ -11,7 +11,7 @@ log_lik_neg_binom(x, lambda = 1, theta = 0)
 
 \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{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson, beta-binomial, and underdispersed Poisson).}
 }
 \value{
 An numeric vector of the corresponding log-likelihoods.
@@ -35,6 +35,7 @@ Other log_lik_dist:
 \code{\link{log_lik_pois_zi}()},
 \code{\link{log_lik_pois}()},
 \code{\link{log_lik_skewnorm}()},
-\code{\link{log_lik_student}()}
+\code{\link{log_lik_student}()},
+\code{\link{log_lik_upois}()}
 }
 \concept{log_lik_dist}
diff --git a/man/log_lik_norm.Rd b/man/log_lik_norm.Rd
index 658e72a9..68a1558b 100644
--- a/man/log_lik_norm.Rd
+++ b/man/log_lik_norm.Rd
@@ -35,6 +35,7 @@ Other log_lik_dist:
 \code{\link{log_lik_pois_zi}()},
 \code{\link{log_lik_pois}()},
 \code{\link{log_lik_skewnorm}()},
-\code{\link{log_lik_student}()}
+\code{\link{log_lik_student}()},
+\code{\link{log_lik_upois}()}
 }
 \concept{log_lik_dist}
diff --git a/man/log_lik_pois.Rd b/man/log_lik_pois.Rd
index 3f1702eb..be8a5ef1 100644
--- a/man/log_lik_pois.Rd
+++ b/man/log_lik_pois.Rd
@@ -33,6 +33,7 @@ Other log_lik_dist:
 \code{\link{log_lik_norm}()},
 \code{\link{log_lik_pois_zi}()},
 \code{\link{log_lik_skewnorm}()},
-\code{\link{log_lik_student}()}
+\code{\link{log_lik_student}()},
+\code{\link{log_lik_upois}()}
 }
 \concept{log_lik_dist}
diff --git a/man/log_lik_pois_zi.Rd b/man/log_lik_pois_zi.Rd
index b8c4fb6e..91e12461 100644
--- a/man/log_lik_pois_zi.Rd
+++ b/man/log_lik_pois_zi.Rd
@@ -35,6 +35,7 @@ Other log_lik_dist:
 \code{\link{log_lik_norm}()},
 \code{\link{log_lik_pois}()},
 \code{\link{log_lik_skewnorm}()},
-\code{\link{log_lik_student}()}
+\code{\link{log_lik_student}()},
+\code{\link{log_lik_upois}()}
 }
 \concept{log_lik_dist}
diff --git a/man/log_lik_skewnorm.Rd b/man/log_lik_skewnorm.Rd
index 1d68f612..2611d452 100644
--- a/man/log_lik_skewnorm.Rd
+++ b/man/log_lik_skewnorm.Rd
@@ -39,6 +39,7 @@ Other log_lik_dist:
 \code{\link{log_lik_norm}()},
 \code{\link{log_lik_pois_zi}()},
 \code{\link{log_lik_pois}()},
-\code{\link{log_lik_student}()}
+\code{\link{log_lik_student}()},
+\code{\link{log_lik_upois}()}
 }
 \concept{log_lik_dist}
diff --git a/man/log_lik_student.Rd b/man/log_lik_student.Rd
index 383b299a..4c272d68 100644
--- a/man/log_lik_student.Rd
+++ b/man/log_lik_student.Rd
@@ -13,7 +13,7 @@ log_lik_student(x, mean = 0, sd = 1, theta = 0)
 
 \item{sd}{A non-negative numeric vector of the standard deviations.}
 
-\item{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson and beta-binomial).}
+\item{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson, beta-binomial, and underdispersed Poisson).}
 }
 \value{
 An numeric vector of the corresponding log-likelihoods.
@@ -37,6 +37,7 @@ Other log_lik_dist:
 \code{\link{log_lik_norm}()},
 \code{\link{log_lik_pois_zi}()},
 \code{\link{log_lik_pois}()},
-\code{\link{log_lik_skewnorm}()}
+\code{\link{log_lik_skewnorm}()},
+\code{\link{log_lik_upois}()}
 }
 \concept{log_lik_dist}
diff --git a/man/log_lik_upois.Rd b/man/log_lik_upois.Rd
new file mode 100644
index 00000000..786cd9f6
--- /dev/null
+++ b/man/log_lik_upois.Rd
@@ -0,0 +1,41 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/log-lik.R
+\name{log_lik_upois}
+\alias{log_lik_upois}
+\title{Underdispersed Poisson Log-Likelihood}
+\usage{
+log_lik_upois(x, lambda = 1, theta = 0)
+}
+\arguments{
+\item{x}{A non-negative whole numeric vector of values.}
+
+\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, beta-binomial, and underdispersed Poisson).}
+}
+\value{
+An numeric vector of the corresponding log-likelihoods.
+}
+\description{
+Underdispersed Poisson Log-Likelihood
+}
+\examples{
+log_lik_upois(c(0, 1, 2), 1, 0)
+}
+\seealso{
+Other log_lik_dist: 
+\code{\link{log_lik_bern}()},
+\code{\link{log_lik_beta_binom}()},
+\code{\link{log_lik_binom}()},
+\code{\link{log_lik_gamma_pois_zi}()},
+\code{\link{log_lik_gamma_pois}()},
+\code{\link{log_lik_gamma}()},
+\code{\link{log_lik_lnorm}()},
+\code{\link{log_lik_neg_binom}()},
+\code{\link{log_lik_norm}()},
+\code{\link{log_lik_pois_zi}()},
+\code{\link{log_lik_pois}()},
+\code{\link{log_lik_skewnorm}()},
+\code{\link{log_lik_student}()}
+}
+\concept{log_lik_dist}
diff --git a/man/params.Rd b/man/params.Rd
index 1c65775b..988758d1 100644
--- a/man/params.Rd
+++ b/man/params.Rd
@@ -36,7 +36,7 @@ should be returned as negative values.}
 
 \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{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson, beta-binomial, and underdispersed Poisson).}
 
 \item{shape}{A non-negative numeric vector of shape.}
 
diff --git a/man/ran_bern.Rd b/man/ran_bern.Rd
index 1dbe3567..c8a8c2dc 100644
--- a/man/ran_bern.Rd
+++ b/man/ran_bern.Rd
@@ -33,6 +33,7 @@ Other ran_dist:
 \code{\link{ran_pois_zi}()},
 \code{\link{ran_pois}()},
 \code{\link{ran_skewnorm}()},
-\code{\link{ran_student}()}
+\code{\link{ran_student}()},
+\code{\link{ran_upois}()}
 }
 \concept{ran_dist}
diff --git a/man/ran_beta_binom.Rd b/man/ran_beta_binom.Rd
index 6f690441..1978a62a 100644
--- a/man/ran_beta_binom.Rd
+++ b/man/ran_beta_binom.Rd
@@ -13,7 +13,7 @@ ran_beta_binom(n = 1, size = 1, prob = 0.5, theta = 0)
 
 \item{prob}{A numeric vector of values between 0 and 1 of the probability of success.}
 
-\item{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson and beta-binomial).}
+\item{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson, beta-binomial, and underdispersed Poisson).}
 }
 \value{
 A numeric vector of the random samples.
@@ -37,6 +37,7 @@ Other ran_dist:
 \code{\link{ran_pois_zi}()},
 \code{\link{ran_pois}()},
 \code{\link{ran_skewnorm}()},
-\code{\link{ran_student}()}
+\code{\link{ran_student}()},
+\code{\link{ran_upois}()}
 }
 \concept{ran_dist}
diff --git a/man/ran_binom.Rd b/man/ran_binom.Rd
index 718ac1f7..b3a773c2 100644
--- a/man/ran_binom.Rd
+++ b/man/ran_binom.Rd
@@ -35,6 +35,7 @@ Other ran_dist:
 \code{\link{ran_pois_zi}()},
 \code{\link{ran_pois}()},
 \code{\link{ran_skewnorm}()},
-\code{\link{ran_student}()}
+\code{\link{ran_student}()},
+\code{\link{ran_upois}()}
 }
 \concept{ran_dist}
diff --git a/man/ran_gamma.Rd b/man/ran_gamma.Rd
index 3ad6a2d3..7c404575 100644
--- a/man/ran_gamma.Rd
+++ b/man/ran_gamma.Rd
@@ -35,6 +35,7 @@ Other ran_dist:
 \code{\link{ran_pois_zi}()},
 \code{\link{ran_pois}()},
 \code{\link{ran_skewnorm}()},
-\code{\link{ran_student}()}
+\code{\link{ran_student}()},
+\code{\link{ran_upois}()}
 }
 \concept{ran_dist}
diff --git a/man/ran_gamma_pois.Rd b/man/ran_gamma_pois.Rd
index 679839df..17658b99 100644
--- a/man/ran_gamma_pois.Rd
+++ b/man/ran_gamma_pois.Rd
@@ -11,7 +11,7 @@ ran_gamma_pois(n = 1, lambda = 1, theta = 0)
 
 \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{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson, beta-binomial, and underdispersed Poisson).}
 }
 \value{
 A numeric vector of the random samples.
@@ -35,6 +35,7 @@ Other ran_dist:
 \code{\link{ran_pois_zi}()},
 \code{\link{ran_pois}()},
 \code{\link{ran_skewnorm}()},
-\code{\link{ran_student}()}
+\code{\link{ran_student}()},
+\code{\link{ran_upois}()}
 }
 \concept{ran_dist}
diff --git a/man/ran_gamma_pois_zi.Rd b/man/ran_gamma_pois_zi.Rd
index e10d76ea..de599f29 100644
--- a/man/ran_gamma_pois_zi.Rd
+++ b/man/ran_gamma_pois_zi.Rd
@@ -11,7 +11,7 @@ ran_gamma_pois_zi(n = 1, lambda = 1, theta = 0, prob = 0)
 
 \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{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson, beta-binomial, and underdispersed Poisson).}
 
 \item{prob}{A numeric vector of values between 0 and 1 of the probability of success.}
 }
@@ -37,6 +37,7 @@ Other ran_dist:
 \code{\link{ran_pois_zi}()},
 \code{\link{ran_pois}()},
 \code{\link{ran_skewnorm}()},
-\code{\link{ran_student}()}
+\code{\link{ran_student}()},
+\code{\link{ran_upois}()}
 }
 \concept{ran_dist}
diff --git a/man/ran_lnorm.Rd b/man/ran_lnorm.Rd
index 3890f45e..8f517109 100644
--- a/man/ran_lnorm.Rd
+++ b/man/ran_lnorm.Rd
@@ -35,6 +35,7 @@ Other ran_dist:
 \code{\link{ran_pois_zi}()},
 \code{\link{ran_pois}()},
 \code{\link{ran_skewnorm}()},
-\code{\link{ran_student}()}
+\code{\link{ran_student}()},
+\code{\link{ran_upois}()}
 }
 \concept{ran_dist}
diff --git a/man/ran_neg_binom.Rd b/man/ran_neg_binom.Rd
index bf4b16d8..b07f20ea 100644
--- a/man/ran_neg_binom.Rd
+++ b/man/ran_neg_binom.Rd
@@ -11,7 +11,7 @@ ran_neg_binom(n = 1, lambda = 1, theta = 0)
 
 \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{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson, beta-binomial, and underdispersed Poisson).}
 }
 \value{
 A numeric vector of the random samples.
@@ -35,6 +35,7 @@ Other ran_dist:
 \code{\link{ran_pois_zi}()},
 \code{\link{ran_pois}()},
 \code{\link{ran_skewnorm}()},
-\code{\link{ran_student}()}
+\code{\link{ran_student}()},
+\code{\link{ran_upois}()}
 }
 \concept{ran_dist}
diff --git a/man/ran_norm.Rd b/man/ran_norm.Rd
index 54dd58ca..8a90d589 100644
--- a/man/ran_norm.Rd
+++ b/man/ran_norm.Rd
@@ -35,6 +35,7 @@ Other ran_dist:
 \code{\link{ran_pois_zi}()},
 \code{\link{ran_pois}()},
 \code{\link{ran_skewnorm}()},
-\code{\link{ran_student}()}
+\code{\link{ran_student}()},
+\code{\link{ran_upois}()}
 }
 \concept{ran_dist}
diff --git a/man/ran_pois.Rd b/man/ran_pois.Rd
index 08bff4da..1f7633cf 100644
--- a/man/ran_pois.Rd
+++ b/man/ran_pois.Rd
@@ -33,6 +33,7 @@ Other ran_dist:
 \code{\link{ran_norm}()},
 \code{\link{ran_pois_zi}()},
 \code{\link{ran_skewnorm}()},
-\code{\link{ran_student}()}
+\code{\link{ran_student}()},
+\code{\link{ran_upois}()}
 }
 \concept{ran_dist}
diff --git a/man/ran_pois_zi.Rd b/man/ran_pois_zi.Rd
index d5e810f9..05e42234 100644
--- a/man/ran_pois_zi.Rd
+++ b/man/ran_pois_zi.Rd
@@ -35,6 +35,7 @@ Other ran_dist:
 \code{\link{ran_norm}()},
 \code{\link{ran_pois}()},
 \code{\link{ran_skewnorm}()},
-\code{\link{ran_student}()}
+\code{\link{ran_student}()},
+\code{\link{ran_upois}()}
 }
 \concept{ran_dist}
diff --git a/man/ran_skewnorm.Rd b/man/ran_skewnorm.Rd
index ea47611a..edfdc440 100644
--- a/man/ran_skewnorm.Rd
+++ b/man/ran_skewnorm.Rd
@@ -39,6 +39,7 @@ Other ran_dist:
 \code{\link{ran_norm}()},
 \code{\link{ran_pois_zi}()},
 \code{\link{ran_pois}()},
-\code{\link{ran_student}()}
+\code{\link{ran_student}()},
+\code{\link{ran_upois}()}
 }
 \concept{ran_dist}
diff --git a/man/ran_student.Rd b/man/ran_student.Rd
index 03b046f6..9a5da01b 100644
--- a/man/ran_student.Rd
+++ b/man/ran_student.Rd
@@ -13,7 +13,7 @@ ran_student(n = 1, mean = 0, sd = 1, theta = 0)
 
 \item{sd}{A non-negative numeric vector of the standard deviations.}
 
-\item{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson and beta-binomial).}
+\item{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson, beta-binomial, and underdispersed Poisson).}
 }
 \value{
 A numeric vector of the random samples.
@@ -37,6 +37,7 @@ Other ran_dist:
 \code{\link{ran_norm}()},
 \code{\link{ran_pois_zi}()},
 \code{\link{ran_pois}()},
-\code{\link{ran_skewnorm}()}
+\code{\link{ran_skewnorm}()},
+\code{\link{ran_upois}()}
 }
 \concept{ran_dist}
diff --git a/man/ran_upois.Rd b/man/ran_upois.Rd
new file mode 100644
index 00000000..1d7d0d1e
--- /dev/null
+++ b/man/ran_upois.Rd
@@ -0,0 +1,41 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ran.R
+\name{ran_upois}
+\alias{ran_upois}
+\title{Underdispersed Poisson Random Samples}
+\usage{
+ran_upois(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, beta-binomial, and underdispersed Poisson).}
+}
+\value{
+A numeric vector of the random samples.
+}
+\description{
+Underdispersed Poisson Random Samples
+}
+\examples{
+ran_upois(n = 10, lambda = 1, theta = 1)
+}
+\seealso{
+Other ran_dist: 
+\code{\link{ran_bern}()},
+\code{\link{ran_beta_binom}()},
+\code{\link{ran_binom}()},
+\code{\link{ran_gamma_pois_zi}()},
+\code{\link{ran_gamma_pois}()},
+\code{\link{ran_gamma}()},
+\code{\link{ran_lnorm}()},
+\code{\link{ran_neg_binom}()},
+\code{\link{ran_norm}()},
+\code{\link{ran_pois_zi}()},
+\code{\link{ran_pois}()},
+\code{\link{ran_skewnorm}()},
+\code{\link{ran_student}()}
+}
+\concept{ran_dist}
diff --git a/man/res_bern.Rd b/man/res_bern.Rd
index 8d37377a..fa862d6a 100644
--- a/man/res_bern.Rd
+++ b/man/res_bern.Rd
@@ -37,6 +37,7 @@ Other res_dist:
 \code{\link{res_pois_zi}()},
 \code{\link{res_pois}()},
 \code{\link{res_skewnorm}()},
-\code{\link{res_student}()}
+\code{\link{res_student}()},
+\code{\link{res_upois}()}
 }
 \concept{res_dist}
diff --git a/man/res_beta_binom.Rd b/man/res_beta_binom.Rd
index 74331467..bdc00fb0 100644
--- a/man/res_beta_binom.Rd
+++ b/man/res_beta_binom.Rd
@@ -20,7 +20,7 @@ res_beta_binom(
 
 \item{prob}{A numeric vector of values between 0 and 1 of the probability of success.}
 
-\item{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson and beta-binomial).}
+\item{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson, beta-binomial, and underdispersed Poisson).}
 
 \item{type}{A string of the residual type. 'raw' for raw residuals 'dev' for deviance residuals and 'data' for the data.}
 
@@ -48,6 +48,7 @@ Other res_dist:
 \code{\link{res_pois_zi}()},
 \code{\link{res_pois}()},
 \code{\link{res_skewnorm}()},
-\code{\link{res_student}()}
+\code{\link{res_student}()},
+\code{\link{res_upois}()}
 }
 \concept{res_dist}
diff --git a/man/res_binom.Rd b/man/res_binom.Rd
index acb7b19b..2e66c045 100644
--- a/man/res_binom.Rd
+++ b/man/res_binom.Rd
@@ -39,6 +39,7 @@ Other res_dist:
 \code{\link{res_pois_zi}()},
 \code{\link{res_pois}()},
 \code{\link{res_skewnorm}()},
-\code{\link{res_student}()}
+\code{\link{res_student}()},
+\code{\link{res_upois}()}
 }
 \concept{res_dist}
diff --git a/man/res_gamma.Rd b/man/res_gamma.Rd
index 15fb6a97..2d57768f 100644
--- a/man/res_gamma.Rd
+++ b/man/res_gamma.Rd
@@ -39,6 +39,7 @@ Other res_dist:
 \code{\link{res_pois_zi}()},
 \code{\link{res_pois}()},
 \code{\link{res_skewnorm}()},
-\code{\link{res_student}()}
+\code{\link{res_student}()},
+\code{\link{res_upois}()}
 }
 \concept{res_dist}
diff --git a/man/res_gamma_pois.Rd b/man/res_gamma_pois.Rd
index 37684a98..83c6caed 100644
--- a/man/res_gamma_pois.Rd
+++ b/man/res_gamma_pois.Rd
@@ -11,7 +11,7 @@ res_gamma_pois(x, lambda = 1, theta = 0, type = "dev", simulate = FALSE)
 
 \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{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson, beta-binomial, and underdispersed Poisson).}
 
 \item{type}{A string of the residual type. 'raw' for raw residuals 'dev' for deviance residuals and 'data' for the data.}
 
@@ -39,6 +39,7 @@ Other res_dist:
 \code{\link{res_pois_zi}()},
 \code{\link{res_pois}()},
 \code{\link{res_skewnorm}()},
-\code{\link{res_student}()}
+\code{\link{res_student}()},
+\code{\link{res_upois}()}
 }
 \concept{res_dist}
diff --git a/man/res_gamma_pois_zi.Rd b/man/res_gamma_pois_zi.Rd
index 6efeb624..6f1b467d 100644
--- a/man/res_gamma_pois_zi.Rd
+++ b/man/res_gamma_pois_zi.Rd
@@ -18,7 +18,7 @@ res_gamma_pois_zi(
 
 \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{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson, beta-binomial, and underdispersed Poisson).}
 
 \item{prob}{A numeric vector of values between 0 and 1 of the probability of success.}
 
@@ -48,6 +48,7 @@ Other res_dist:
 \code{\link{res_pois_zi}()},
 \code{\link{res_pois}()},
 \code{\link{res_skewnorm}()},
-\code{\link{res_student}()}
+\code{\link{res_student}()},
+\code{\link{res_upois}()}
 }
 \concept{res_dist}
diff --git a/man/res_lnorm.Rd b/man/res_lnorm.Rd
index 1e72b5a0..6bcc08ac 100644
--- a/man/res_lnorm.Rd
+++ b/man/res_lnorm.Rd
@@ -39,6 +39,7 @@ Other res_dist:
 \code{\link{res_pois_zi}()},
 \code{\link{res_pois}()},
 \code{\link{res_skewnorm}()},
-\code{\link{res_student}()}
+\code{\link{res_student}()},
+\code{\link{res_upois}()}
 }
 \concept{res_dist}
diff --git a/man/res_neg_binom.Rd b/man/res_neg_binom.Rd
index 5d6f38f7..0e1d4be7 100644
--- a/man/res_neg_binom.Rd
+++ b/man/res_neg_binom.Rd
@@ -11,7 +11,7 @@ res_neg_binom(x, lambda = 1, theta = 0, type = "dev", simulate = FALSE)
 
 \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{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson, beta-binomial, and underdispersed Poisson).}
 
 \item{type}{A string of the residual type. 'raw' for raw residuals 'dev' for deviance residuals and 'data' for the data.}
 
@@ -39,6 +39,7 @@ Other res_dist:
 \code{\link{res_pois_zi}()},
 \code{\link{res_pois}()},
 \code{\link{res_skewnorm}()},
-\code{\link{res_student}()}
+\code{\link{res_student}()},
+\code{\link{res_upois}()}
 }
 \concept{res_dist}
diff --git a/man/res_norm.Rd b/man/res_norm.Rd
index 513dc4d4..1d683177 100644
--- a/man/res_norm.Rd
+++ b/man/res_norm.Rd
@@ -39,6 +39,7 @@ Other res_dist:
 \code{\link{res_pois_zi}()},
 \code{\link{res_pois}()},
 \code{\link{res_skewnorm}()},
-\code{\link{res_student}()}
+\code{\link{res_student}()},
+\code{\link{res_upois}()}
 }
 \concept{res_dist}
diff --git a/man/res_pois.Rd b/man/res_pois.Rd
index 54e1e845..efff57d9 100644
--- a/man/res_pois.Rd
+++ b/man/res_pois.Rd
@@ -37,6 +37,7 @@ Other res_dist:
 \code{\link{res_norm}()},
 \code{\link{res_pois_zi}()},
 \code{\link{res_skewnorm}()},
-\code{\link{res_student}()}
+\code{\link{res_student}()},
+\code{\link{res_upois}()}
 }
 \concept{res_dist}
diff --git a/man/res_pois_zi.Rd b/man/res_pois_zi.Rd
index 68c907df..010cc935 100644
--- a/man/res_pois_zi.Rd
+++ b/man/res_pois_zi.Rd
@@ -39,6 +39,7 @@ Other res_dist:
 \code{\link{res_norm}()},
 \code{\link{res_pois}()},
 \code{\link{res_skewnorm}()},
-\code{\link{res_student}()}
+\code{\link{res_student}()},
+\code{\link{res_upois}()}
 }
 \concept{res_dist}
diff --git a/man/res_skewnorm.Rd b/man/res_skewnorm.Rd
index 88527684..010594ce 100644
--- a/man/res_skewnorm.Rd
+++ b/man/res_skewnorm.Rd
@@ -41,6 +41,7 @@ Other res_dist:
 \code{\link{res_norm}()},
 \code{\link{res_pois_zi}()},
 \code{\link{res_pois}()},
-\code{\link{res_student}()}
+\code{\link{res_student}()},
+\code{\link{res_upois}()}
 }
 \concept{res_dist}
diff --git a/man/res_student.Rd b/man/res_student.Rd
index e7c430d5..435b76e5 100644
--- a/man/res_student.Rd
+++ b/man/res_student.Rd
@@ -13,7 +13,7 @@ res_student(x, mean = 0, sd = 1, theta = 0, type = "dev", simulate = FALSE)
 
 \item{sd}{A non-negative numeric vector of the standard deviations.}
 
-\item{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson and beta-binomial).}
+\item{theta}{A non-negative numeric vector of the dispersion for the mixture models (student, gamma-Poisson, beta-binomial, and underdispersed Poisson).}
 
 \item{type}{A string of the residual type. 'raw' for raw residuals 'dev' for deviance residuals and 'data' for the data.}
 
@@ -41,6 +41,7 @@ Other res_dist:
 \code{\link{res_norm}()},
 \code{\link{res_pois_zi}()},
 \code{\link{res_pois}()},
-\code{\link{res_skewnorm}()}
+\code{\link{res_skewnorm}()},
+\code{\link{res_upois}()}
 }
 \concept{res_dist}
diff --git a/man/res_upois.Rd b/man/res_upois.Rd
new file mode 100644
index 00000000..0b32c5fd
--- /dev/null
+++ b/man/res_upois.Rd
@@ -0,0 +1,45 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/res.R
+\name{res_upois}
+\alias{res_upois}
+\title{Underdispersed Poisson Residuals}
+\usage{
+res_upois(x, lambda = 1, theta = 0, type = "dev", simulate = FALSE)
+}
+\arguments{
+\item{x}{A non-negative whole numeric vector of values.}
+
+\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, beta-binomial, and underdispersed Poisson).}
+
+\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{
+Underdispersed Poisson Residuals
+}
+\examples{
+res_upois(c(0, 1, 2), 1, 1)
+}
+\seealso{
+Other res_dist: 
+\code{\link{res_bern}()},
+\code{\link{res_beta_binom}()},
+\code{\link{res_binom}()},
+\code{\link{res_gamma_pois_zi}()},
+\code{\link{res_gamma_pois}()},
+\code{\link{res_gamma}()},
+\code{\link{res_lnorm}()},
+\code{\link{res_neg_binom}()},
+\code{\link{res_norm}()},
+\code{\link{res_pois_zi}()},
+\code{\link{res_pois}()},
+\code{\link{res_skewnorm}()},
+\code{\link{res_student}()}
+}
+\concept{res_dist}
diff --git a/tests/testthat/test-dev.R b/tests/testthat/test-dev.R
index 415c3349..795de1ba 100644
--- a/tests/testthat/test-dev.R
+++ b/tests/testthat/test-dev.R
@@ -768,3 +768,77 @@ test_that("skewnorm deviance", {
   deviance <- sum(dev_skewnorm(samples, coef(mod)[1]))
   expect_equal(deviance, deviance(mod))
 })
+
+test_that("upois known values", {
+  expect_identical(dev_upois(1, 1), 0)
+  expect_equal(dev_upois(1, 1, 1), 0.424635855096438)
+  expect_identical(dev_upois(0, 0), 0)
+  expect_identical(dev_upois(0, 0, 1), 0)
+  expect_identical(dev_upois(0, 0, 0), dev_pois(0, 0))
+  expect_equal(dev_upois(0, 1, 0), dev_pois(0, 1))
+  expect_equal(dev_upois(1, 0), Inf)
+  expect_equal(dev_upois(1, 0, 1), 0)
+  expect_equal(dev_upois(1, 0, 5), 0)
+  expect_equal(dev_upois(0, 1, 5), 2)
+  expect_equal(dev_upois(0, 1, 1), 2)
+  expect_equal(dev_upois(20, 1, 1), 77.139306160016)
+  expect_equal(dev_upois(10, 1, 5), 23.8440589297922)
+  expect_equal(dev_upois(0, 2, 5), 4)
+  expect_equal(dev_upois(0, 2, 1), 4)
+  expect_equal(dev_upois(3, 2, 1), 0.0831945648654893)
+  expect_equal(dev_upois(13, 2, 1), 24.0432116320838)
+  expect_equal(dev_upois(1, 2, 2), 2.25534033186796)
+  expect_equal(dev_upois(10, 2, 2), 13.6367082152856)
+  expect_equal(dev_upois(100, 2, 2), 579.376048918714)
+  expect_equal(dev_upois(10, 2, 0), dev_pois(10, 2))
+  expect_equal(dev_upois(100, 2, 0), dev_pois(100, 2))
+})
+
+test_that("upois vectorized", {
+  expect_equal(dev_upois(0:3, 2, 0.1),
+               c(4, 0.716043163399263, 0.00437744156764275, 0.346749115995378))
+  expect_equal(dev_upois(c(0, 1, 3, 0), 3, 0.5),
+               c(6, 2.46944208933045, 0.0416124417688755, 6))
+  expect_equal(dev_upois(0:3, 0:3, rep(1, 4)),
+               c(0, 0.424635855096438, 0.157573069847393, 0.0966910188449219
+  ))
+  expect_equal(dev_upois(0:3, 3:0, 0:3), c(6, 1.61370563888011, 0.0377745205120377, Inf))
+})
+
+test_that("upois vectorized missing values", {
+  expect_equal(dev_upois(c(NA,1), 0:1, 0:1), c(NA,0.424635855096438))
+  expect_equal(dev_upois(c(0,NA), 0:1, 1:2), c(0,NA))
+  expect_equal(dev_upois(c(0:1), c(NA,1), 1:2), c(NA, 0.830704476771521))
+  expect_equal(dev_upois(c(0:1), c(0,NA), 1:2), c(0,NA))
+  expect_equal(dev_upois(c(0:1), c(0:1), c(NA,1)), c(NA, 0.424635855096438))
+  expect_equal(dev_upois(c(0:1), c(0:1), c(1,NA)), c(0,NA))
+})
+
+test_that("upois res", {
+  expect_equal(dev_upois(10, 0.5, 0.5), dev_upois(10, 0.5, 0.5, res = TRUE)^2)
+  expect_equal(dev_upois(0:1, c(0.3,0.6), 0.5), dev_upois(0:1, c(0.3,0.6), 0.5, res = TRUE)^2)
+})
+
+test_that("upois log_lik", {
+  expect_equal(dev_upois(5:6, 1:2, 2:3),
+               2 * (log_lik_upois(5:6, 5:6 - (2:3 / (1 + 2:3)), 2:3) -
+                      log_lik_upois(5:6, 1:2, 2:3)))
+})
+
+test_that("upois ran", {
+  set.seed(101)
+  samples <- ran_upois(100000, 3, 0.5)
+  expect_equal(mean(samples), 3.33127)
+  expect_equal(var(samples), 3.18490203611947)
+  res <- dev_upois(samples, 3, 0.5, res = TRUE)
+  expect_equal(mean(res), -0.113602120823609)
+  expect_equal(sd(res), 1.02371893617533)
+})
+
+test_that("upois deviance", {
+  # Poisson case
+  samples <- ran_upois(1000, 3, 0)
+  mod <- glm(samples~1, family = poisson)
+  deviance <- sum(dev_upois(samples, exp(coef(mod)[1])))
+  expect_equal(deviance, deviance(mod))
+})
diff --git a/tests/testthat/test-log-lik.R b/tests/testthat/test-log-lik.R
index f64fb18e..601a2fb6 100644
--- a/tests/testthat/test-log-lik.R
+++ b/tests/testthat/test-log-lik.R
@@ -243,3 +243,43 @@ test_that("skewnorm vectorized", {
     c(-111.116353440211, -11.8095006207706, -11.1163537268622, -111.116353440211)
   )
 })
+
+test_that("upois missing values", {
+  expect_identical(log_lik_upois(numeric(0), numeric(0), numeric(0)), numeric(0))
+  expect_identical(log_lik_upois(1, lambda = numeric(0)), numeric(0))
+  expect_identical(log_lik_upois(1, 1, theta = numeric(0)), numeric(0))
+  expect_identical(log_lik_upois(NA, 1, 1), NA_real_)
+  expect_identical(log_lik_upois(1, NA, 1), NA_real_)
+  expect_identical(log_lik_upois(1, 1, NA), NA_real_)
+})
+
+test_that("upois known values", {
+  expect_identical(log_lik_upois(1, 2), dpois(1, 2, log = TRUE))
+  expect_identical(log_lik_upois(0, 2), dpois(0, 2, log = TRUE))
+  expect_equal(log_lik_upois(0, 2, 1), -2.69314718055995)
+  expect_identical(log_lik_upois(1, 0, 0), -Inf)
+  expect_identical(log_lik_upois(0, 0, 0), 0)
+  expect_identical(log_lik_upois(1, 0, 10), log(10 / (1 + 10)))
+  expect_identical(log_lik_upois(0, 0, 10), log(1 / (1 + 10)))
+  expect_identical(log_lik_upois(100, 1000, 0), dpois(100, 1000, log = TRUE))
+  expect_equal(log_lik_upois(3, 3.5, 0), log_lik_pois(3, 3.5))
+  expect_equal(log_lik_upois(3, 3.5, 0), -1.53347056374195)
+  expect_equal(log_lik_upois(3, 3.5, 0.1), -1.5465426453093)
+  expect_equal(log_lik_upois(3, 3.5, 0.2), -1.55756811532101)
+  expect_equal(log_lik_upois(3, 3.5, 1), -1.60757853589567)
+})
+
+test_that("upois vectorized", {
+  expect_identical(log_lik_upois(0:1, 0:1, 0:1), c(0, -1))
+  expect_equal(log_lik_upois(0:1, 0:1, 1:2), c(-0.693147180559945, -1))
+  expect_identical(log_lik_upois(0:1, 0:1, 0), dpois(0:1, 0:1, log = TRUE))
+  expect_identical(log_lik_upois(0:1, 1000, 0), dpois(0:1, 1000, log = TRUE))
+  expect_equal(log_lik_upois(c(0, 2), 2, 0.5), c(-2.40546510810816, -1.30685281944005))
+  expect_equal(log_lik_upois(c(0, 1, 3), c(0, 0, 1), 1),
+                   c(-0.693147180559945, -0.693147180559945, -2.09861228866811))
+  expect_equal(log_lik_upois(c(0, 1, 3), c(0, 0, 1), 2),
+                   c(-1.09861228866811, -0.405465108108164, -1.94446160884085))
+  expect_equal(log_lik_upois(c(19, 1000, 21), c(10, 2000, 10), c(0, 1, 10)),
+               c(-5.59076742031263, -311.513401018534, -6.33270476504199))
+})
+
diff --git a/tests/testthat/test-ran.R b/tests/testthat/test-ran.R
index 952cce80..e6e1966e 100644
--- a/tests/testthat/test-ran.R
+++ b/tests/testthat/test-ran.R
@@ -278,3 +278,41 @@ test_that("ran_skewnorm", {
   expect_equal(ran_skewnorm(2, c(10,100), sd = 0, shape = -1), c(10, 100))
   expect_equal(ran_skewnorm(1, c(10,100), sd = 10, shape = -1), c(-2.24639864962439))
 })
+
+
+test_that("ran_upois", {
+  expect_error(ran_upois(NA_integer_))
+  expect_error(ran_upois(NULL))
+  expect_error(ran_upois(integer(0)))
+  expect_error(ran_upois(-1))
+  expect_error(ran_upois(1, -1, 1))
+  expect_error(ran_upois(1, 1, -1))
+  expect_identical(ran_upois(0L), integer(0))
+  expect_warning(ran_upois(1L, 1:3, 1:2), c("longer argument not a multiple of length of shorter"))
+  set.seed(101)
+  expect_identical(ran_upois(1, 1, 0), 0L)
+  expect_identical(ran_upois(1, 2, 0), 1L)
+  expect_identical(ran_upois(1, 0, 0), 0L)
+  expect_identical(ran_upois(10, 1, 0), c(0L, 1L, 2L, 0L, 1L, 1L, 1L, 3L, 3L, 0L))
+  expect_identical(length(ran_upois(10, 1, 0)), 10L)
+  expect_identical(ran_upois(1, 1000, 0), 1056L)
+  # Vectorized
+  set.seed(101)
+  expect_identical(ran_upois(2, c(1, 2), 0), c(1L, 1L))
+  expect_identical(ran_upois(2, c(1, 1000), 0), c(2L, 1020L))
+  expect_identical(ran_upois(10, 10, seq(0, 9, 1)), c(9L, 13L, 10L, 12L, 8L, 16L, 6L, 10L, 10L, 8L))
+  # Correct mean
+  set.seed(101)
+  x <- ran_upois(10000, 10, 0)
+  expect_equal(mean(x), 10, tolerance = 0.01)
+  set.seed(101)
+  x <- ran_upois(10000, 10, 1)
+  expect_equal(mean(x), 10 + (1 / (1 + 1)), tolerance = 0.01)
+  # Correct variance
+  set.seed(101)
+  x <- ran_upois(10000, 10, 0)
+  expect_equal(var(x), 10, tolerance = 0.01)
+  set.seed(101)
+  x <- ran_upois(10000, 10, 1)
+  expect_equal(var(x), 10 + (1 / (1 + 1)^2), tolerance = 0.01)
+})
diff --git a/tests/testthat/test-res.R b/tests/testthat/test-res.R
index d7b02dfb..1f2c4cd8 100644
--- a/tests/testthat/test-res.R
+++ b/tests/testthat/test-res.R
@@ -570,3 +570,36 @@ test_that("res_skewnorm", {
   expect_equal(mean(res), 0.0236707225149864)
   expect_equal(sd(res), 3.34458775450197)
 })
+
+test_that("res_upois", {
+  expect_identical(res_upois(integer(0), integer(0), integer(0)), numeric(0))
+  expect_equal(res_upois(1, 1, 1), -0.651640894278772)
+  expect_identical(res_upois(1, 1, 1, type = "raw"), -0.5)
+  expect_equal(res_upois(1, 1, 0.5, type = "raw"), -1/3)
+  expect_equal(res_upois(1, 1, 0.5), -0.405015814635496)
+  expect_equal(res_upois(2, 2, 2, type = "raw"), -2/3)
+  expect_equal(res_upois(0, 1, 0, type = "raw"), -1)
+  expect_identical(res_upois(NA, 1, 1), NA_real_)
+  expect_identical(res_upois(1, NA, 1), NA_real_)
+  expect_identical(res_upois(1, 1, NA), NA_real_)
+  expect_error(res_upois(1, 3, 1, type = "unknown"))
+  expect_equal(res_upois(1, 3, 1, type = "raw"), -2.5)
+  expect_equal(res_upois(c(1,3.5,4), 3, 10, type = "raw"),
+               c(-2.90909090909091, -0.409090909090909, 0.0909090909090908))
+  expect_equal(res_upois(c(1,3.5,4), 3, 5, type = "raw"),
+               c(-2.83333333333333, -0.333333333333333, 0.166666666666667))
+  set.seed(101)
+  expect_equal(res_upois(1:2, 2, 2, simulate = TRUE, type = "raw"), c(-0.666666666666667, -2.66666666666667))
+  expect_equal(res_upois(1:2, 2, 2, simulate = TRUE), c(-0.535571699965782, -0.535571699965782))
+  set.seed(101)
+  expect_equal(res_upois(1:2, 2, 0.5, simulate = TRUE, type = "raw"), c(-0.333333333333333, -2.33333333333333))
+  expect_equal(res_upois(1:2, 2, 0.5, simulate = TRUE), c(-1.06882485121411, -1.06882485121411))
+  set.seed(101)
+  expect_equal(res_upois(1:2, 2, 1, simulate = TRUE, type = "data"), c(2L, 0L))
+  set.seed(101)
+  expect_equal(res_upois(1:2, 2, 1, simulate = TRUE, type = "standardized"), c(-0.333333333333333, -1.66666666666667))
+  set.seed(101)
+  res <- res_upois(rep(2, 10000), 2, 1, simulate = TRUE, type = "standardized")
+  expect_equal(mean(res), 0.00306666666666664)
+  expect_equal(sd(res), 1.00071178233097)
+})