diff --git a/DESCRIPTION b/DESCRIPTION
index ad6aa191..3e7384ec 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -33,6 +33,7 @@ Suggests:
     covr,
     extraDistr,
     ggplot2,
+    glmnet,
     hms,
     knitr,
     memoise,
diff --git a/NAMESPACE b/NAMESPACE
index 1cb86856..0ac22bcf 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -39,6 +39,7 @@ export(dev_gamma)
 export(dev_gamma_pois)
 export(dev_gamma_pois_zi)
 export(dev_lnorm)
+export(dev_multinom)
 export(dev_neg_binom)
 export(dev_norm)
 export(dev_pois)
@@ -72,6 +73,7 @@ export(log_lik_gamma)
 export(log_lik_gamma_pois)
 export(log_lik_gamma_pois_zi)
 export(log_lik_lnorm)
+export(log_lik_multinom)
 export(log_lik_neg_binom)
 export(log_lik_norm)
 export(log_lik_pois)
@@ -150,6 +152,7 @@ export(ran_gamma)
 export(ran_gamma_pois)
 export(ran_gamma_pois_zi)
 export(ran_lnorm)
+export(ran_multinom)
 export(ran_neg_binom)
 export(ran_norm)
 export(ran_pois)
@@ -165,6 +168,7 @@ export(res_gamma)
 export(res_gamma_pois)
 export(res_gamma_pois_zi)
 export(res_lnorm)
+export(res_multinom)
 export(res_neg_binom)
 export(res_norm)
 export(res_pois)
diff --git a/R/dev.R b/R/dev.R
index 641286c5..d9fda82f 100644
--- a/R/dev.R
+++ b/R/dev.R
@@ -267,6 +267,41 @@ dev_lnorm <- function(x, meanlog = 0, sdlog = 1, res = FALSE) {
   dev_norm(log(x), mean = meanlog, sd = sdlog, res = res)
 }
 
+#' Multinomial Deviances
+#'
+#' Models the counts across two or more mutually exclusive categories from a
+#' fixed number of trials, in \emph{long} format: one row per category per
+#' trial, with `group` identifying which rows belong to the same trial.
+#'
+#' A category's deviance depends only on its own `x` and `mu = size * prob`,
+#' not on the rest of its trial, so `group` is used only to validate `size`
+#' and `prob` (see [log_lik_multinom()]), not in the calculation itself.
+#' `dev_multinom()` is the Poisson-equivalent deviance (see [dev_pois()]):
+#' summing it over a trial's rows recovers the trial's exact multinomial
+#' deviance.
+#'
+#' @inheritParams params
+#' @param x A non-negative whole numeric vector of the category counts.
+#' @param prob A numeric vector of the probability of the category. Must sum
+#'   to 1 across the rows sharing the same `group`.
+#'
+#' @return An numeric vector of the corresponding deviances or deviance residuals.
+#' @family dev_dist
+#' @export
+#'
+#' @examples
+#' dev_multinom(c(1, 3, 6), size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1))
+dev_multinom <- function(x, size = 1, prob, group, res = FALSE) {
+  chk_compatible_lengths(x, size, prob, group)
+  n <- length(x)
+  size <- rep_len(size, n)
+  prob <- rep_len(prob, n)
+  group <- rep_len(group, n)
+  chk_not_any_na(group)
+  chk_multinom_group(size, prob, group)
+  dev_pois(x, lambda = size * prob, res = res)
+}
+
 #' Negative Binomial Deviances
 #'
 #' @inheritParams params
diff --git a/R/internal.R b/R/internal.R
index a56e483b..2c000a3f 100644
--- a/R/internal.R
+++ b/R/internal.R
@@ -1,3 +1,84 @@
 dev_res <- function(x, mu, dev) {
   sign(x - mu) * sqrt(dev)
 }
+
+# Row indices for each `group`, shared across the multinom_* helpers within
+# a single call (chk_multinom_group(), multinom_row_na(), the sampling loop
+# in ran_multinom()) so `group` isn't re-split by every one of them.
+multinom_split <- function(group) {
+  split(seq_along(group), group)
+}
+
+# Checks every group shares one `size` and `prob` values summing to 1
+# (required by rmultinom() and the deviance/log-lik identities), has >= 2
+# rows (a trial needs >= 2 categories -- singletons usually mean `group`
+# was evaluated row-by-row instead of over the whole vector), and matches
+# the modal row count across groups (a short group usually means a row was
+# lost). Only non-NA values are compared, so lone NAs don't error here --
+# see multinom_row_na(). Callers must chk_not_any_na(group) first; `group`
+# itself can't be NA-tolerant since it's what identifies the trial.
+chk_multinom_group <- function(size, prob, group, groups = multinom_split(group)) {
+  for (idx in groups) {
+    if (length(idx) < 2L) {
+      stop(
+        "Each `group` must contain at least 2 rows (a multinomial trial needs at least 2 categories); found a group with only 1 row. This usually means `group`/`size`/`prob` were passed one row at a time instead of as vectors.",
+        call. = FALSE
+      )
+    }
+    known_size <- size[idx][!is.na(size[idx])]
+    if (length(unique(known_size)) > 1L) {
+      stop(
+        "`size` must be the same for every row belonging to the same `group` (multinomial trial).",
+        call. = FALSE
+      )
+    }
+    known_prob <- prob[idx][!is.na(prob[idx])]
+    known_prob_sum <- sum(known_prob)
+    # a group with a missing prob can only be validated one-sided: the known
+    # values must not already exceed 1, since a full sum-to-1 check would be
+    # (wrongly) skipped whenever any prob in the group is NA
+    prob_bad <- if (length(known_prob) < length(idx)) {
+      known_prob_sum > 1 + 1e-6
+    } else {
+      abs(known_prob_sum - 1) > 1e-6
+    }
+    if (prob_bad) {
+      stop(
+        "`prob` must sum to 1 for every `group` (multinomial trial).",
+        call. = FALSE
+      )
+    }
+  }
+  if (length(groups) > 1L) {
+    group_sizes <- lengths(groups)
+    size_counts <- table(group_sizes)
+    # ties are broken in favour of the smallest row count (table()'s names
+    # are sorted ascending, and which.max() takes the first maximum)
+    mode_size <- as.integer(names(size_counts)[which.max(size_counts)])
+    bad <- group_sizes != mode_size
+    if (any(bad)) {
+      stop(
+        sprintf(
+          "Every `group` should have the same number of rows (%d, the most common number of categories in this data); found a group (\"%s\") with %d row(s) instead. This usually means `group` lost a row that should have been part of that trial.",
+          mode_size,
+          names(groups)[bad][1],
+          group_sizes[bad][1]
+        ),
+        call. = FALSE
+      )
+    }
+  }
+}
+
+# Flags every row whose trial has an NA `size`/`prob` anywhere in the group,
+# since a trial's categories are scored/drawn jointly, not independently.
+multinom_row_na <- function(size, prob, group, groups = multinom_split(group)) {
+  bad <- is.na(size) | is.na(prob)
+  result <- rep(FALSE, length(group))
+  for (idx in groups) {
+    if (any(bad[idx])) {
+      result[idx] <- TRUE
+    }
+  }
+  result
+}
diff --git a/R/log-lik.R b/R/log-lik.R
index c579fe64..044a61a0 100644
--- a/R/log-lik.R
+++ b/R/log-lik.R
@@ -427,6 +427,54 @@ log_lik_lnorm <- function(x, meanlog = 0, sdlog = 1, tlower = 0, tupper = Inf) {
   log_lik
 }
 
+#' Multinomial Log-Likelihood
+#'
+#' Models the counts across two or more mutually exclusive categories from a
+#' fixed number of trials, in \emph{long} format: one row per category per
+#' trial, with `group` identifying which rows belong to the same trial. All
+#' rows sharing a `group` must have the same `size`, and their `prob` values
+#' must sum to 1.
+#'
+#' A trial's log-likelihood doesn't split evenly across its rows, since the
+#' multinomial coefficient belongs to the whole trial. `log_lik_multinom()`
+#' uses the multinomial-as-independent-Poissons identity: each row's value
+#' is the Poisson log-likelihood of `x` given `mu = size * prob`, minus an
+#' even share of the trial's normalizing constant, so summing over a
+#' `group` recovers the trial's exact multinomial log-likelihood.
+#'
+#' @inheritParams params
+#' @param x A non-negative whole numeric vector of the category counts.
+#' @param prob A numeric vector of the probability of the category. Must sum
+#'   to 1 across the rows sharing the same `group`. `NA` in `size` or `prob`
+#'   for any row of a trial makes the log-likelihood `NA` for every row of
+#'   that trial, since a trial's categories are scored jointly.
+#'
+#' @return An numeric vector of the corresponding log-likelihoods, one value
+#'   per row of `x`.
+#' @family log_lik_dist
+#' @export
+#'
+#' @examples
+#' log_lik_multinom(c(1, 3, 6), size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1))
+log_lik_multinom <- function(x, size = 1, prob, group) {
+  chk_compatible_lengths(x, size, prob, group)
+  n <- length(x)
+  size <- rep_len(size, n)
+  prob <- rep_len(prob, n)
+  group <- rep_len(group, n)
+  chk_not_any_na(group)
+  groups <- multinom_split(group)
+  chk_multinom_group(size, prob, group, groups)
+  mu <- size * prob
+  log_lik <- log_lik_pois(x, mu)
+  group_size <- table(group)
+  k <- as.numeric(group_size[as.character(group)])
+  const <- log_lik_pois(size, size)
+  log_lik <- log_lik - const / k
+  log_lik[multinom_row_na(size, prob, group, groups)] <- NA_real_
+  log_lik
+}
+
 #' Negative Binomial Log-Likelihood
 #'
 #' @inheritParams params
diff --git a/R/params.R b/R/params.R
index ec79aed9..357c8ab5 100644
--- a/R/params.R
+++ b/R/params.R
@@ -16,6 +16,11 @@
 #'   level.
 #' @param directional A flag specifying whether probabilities less than 0.5
 #'   should be returned as negative values.
+#' @param group A vector identifying which rows belong to the same
+#'   multinomial trial (whose `x` values sum to `size` and `prob` values sum
+#'   to 1). Every group must have at least 2 rows and the same number of
+#'   rows as the rest of the data (a fixed set of categories, as in
+#'   multinomial logistic regression), and must not contain `NA`.
 #' @param lambda A non-negative numeric vector of means.
 #' @param level A number > 0 and <= 1 specifying the probability coverage of the
 #' interval.
diff --git a/R/ran.R b/R/ran.R
index 65369c77..0d5863ce 100644
--- a/R/ran.R
+++ b/R/ran.R
@@ -113,6 +113,50 @@ ran_lnorm <- function(n = 1, meanlog = 0, sdlog = 1) {
   stats::rlnorm(n, meanlog = meanlog, sdlog = sdlog)
 }
 
+#' Multinomial Random Samples
+#'
+#' Models the counts across two or more mutually exclusive categories from a
+#' fixed number of trials, in \emph{long} format: one value per category per
+#' trial, with `group` identifying which rows belong to the same trial. All
+#' rows sharing a `group` must have the same `size`, and their `prob` values
+#' must sum to 1.
+#'
+#' Unlike the other `ran_*()` functions, `ran_multinom()` has no `n`
+#' argument: the number of samples is fully determined by `length(prob)`
+#' (equivalently `length(group)`), since a trial's categories can't be
+#' generated independently of one another.
+#'
+#' @inheritParams params
+#' @param prob A numeric vector of the probability of the category. Must sum
+#'   to 1 across the rows sharing the same `group`. `NA` in `size` or `prob`
+#'   for any row of a trial makes the sample `NA` for every row of that
+#'   trial, since a trial's categories are drawn jointly.
+#' @return An integer vector of the random samples, one per row of `prob`.
+#' @family ran_dist
+#' @export
+#'
+#' @examples
+#' ran_multinom(size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1))
+ran_multinom <- function(size = 1, prob, group) {
+  chk_compatible_lengths(size, prob, group)
+  n <- length(prob)
+  size <- rep_len(size, n)
+  prob <- rep_len(prob, n)
+  group <- rep_len(group, n)
+  chk_not_any_na(group)
+  groups <- multinom_split(group)
+  chk_multinom_group(size, prob, group, groups)
+  row_na <- multinom_row_na(size, prob, group, groups)
+  x <- rep(NA_real_, n)
+  for (idx in groups) {
+    if (row_na[idx[1]]) {
+      next
+    }
+    x[idx] <- stats::rmultinom(1, size = size[idx[1]], prob = prob[idx])[, 1]
+  }
+  as.integer(x)
+}
+
 #' Negative Binomial Random Samples
 #'
 #' Identical to Gamma-Poisson Random Samples.
diff --git a/R/res.R b/R/res.R
index c260d5e5..8b9e16fe 100644
--- a/R/res.R
+++ b/R/res.R
@@ -251,6 +251,65 @@ res_lnorm <- function(
   )
 }
 
+#' Multinomial Residuals
+#'
+#' Models the counts across two or more mutually exclusive categories from a
+#' fixed number of trials, in \emph{long} format: one row per category per
+#' trial, with `group` identifying which rows belong to the same trial (see
+#' [log_lik_multinom()] for details). `res_multinom()` returns one residual
+#' per row, not one per trial, since a trial's categories aren't
+#' independent and so have no single meaningful residual as a whole; the
+#' classic per-trial deviance statistic can be recovered by summing the
+#' squared `type = "dev"` residuals within a `group`.
+#'
+#' `group` is validated (same `size`, `prob` summing to 1, no singleton or
+#' short groups, no `NA`) regardless of `simulate`, but is only otherwise
+#' used when `simulate = TRUE`, to draw a joint, correlation-preserving
+#' replicate per trial (via [ran_multinom()]) rather than simulating each
+#' category independently, which requires `res_multinom()` to see every row
+#' of a `group` in the same call.
+#'
+#' @inheritParams params
+#' @param x A non-negative whole numeric vector of the category counts.
+#' @param prob A numeric vector of the probability of the category. Must sum
+#'   to 1 across the rows sharing the same `group`.
+#'
+#' @return An numeric vector of the corresponding residuals.
+#' @family res_dist
+#' @export
+#'
+#' @examples
+#' res_multinom(c(1, 3, 6), size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1))
+res_multinom <- function(
+  x,
+  size = 1,
+  prob,
+  group,
+  type = "dev",
+  simulate = FALSE
+) {
+  chk_string(type)
+  chk_compatible_lengths(x, size, prob, group)
+  n <- length(x)
+  size <- rep_len(size, n)
+  prob <- rep_len(prob, n)
+  group <- rep_len(group, n)
+  chk_not_any_na(group)
+  chk_multinom_group(size, prob, group)
+  if (!vld_false(simulate)) {
+    x <- ran_multinom(size = size, prob = prob, group = group)
+  }
+  mu <- size * prob
+  switch(
+    type,
+    data = x,
+    raw = x - mu,
+    standardized = (x - mu) / sqrt(mu * (1 - prob)),
+    dev = dev_multinom(x, size = size, prob = prob, group = group, res = TRUE),
+    chk_subset(x, c("data", "raw", "dev", "standardized"))
+  )
+}
+
 #' Negative Binomial Residuals
 #'
 #' @inheritParams params
diff --git a/_pkgdown.yml b/_pkgdown.yml
index df964f30..630b6e77 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -106,6 +106,7 @@ reference:
   - '`dev_gamma_pois`'
   - '`dev_gamma_pois_zi`'
   - '`dev_lnorm`'
+  - '`dev_multinom`'
   - '`dev_neg_binom`'
   - '`dev_norm`'
   - '`dev_pois`'
@@ -123,6 +124,7 @@ reference:
   - '`res_gamma_pois`'
   - '`res_gamma_pois_zi`'
   - '`res_lnorm`'
+  - '`res_multinom`'
   - '`res_neg_binom`'
   - '`res_norm`'
   - '`res_pois`'
@@ -142,6 +144,7 @@ reference:
   - '`log_lik_gamma_pois_zi`'
   - '`log_lik_exp`'
   - '`log_lik_lnorm`'
+  - '`log_lik_multinom`'
   - '`log_lik_neg_binom`'
   - '`log_lik_norm`'
   - '`log_lik_pois`'
@@ -199,6 +202,7 @@ reference:
   - '`ran_gamma_pois`'
   - '`ran_gamma_pois_zi`'
   - '`ran_lnorm`'
+  - '`ran_multinom`'
   - '`ran_neg_binom`'
   - '`ran_norm`'
   - '`ran_pois`'
diff --git a/inst/WORDLIST b/inst/WORDLIST
index 127f2b03..8257c2a9 100644
--- a/inst/WORDLIST
+++ b/inst/WORDLIST
@@ -11,6 +11,7 @@ Numericise
 Numericize
 ORCID
 POSIXct
+Poissons
 Psychol
 Schaub
 Shachar
diff --git a/man/dev_bern.Rd b/man/dev_bern.Rd
index f978015f..91e829b8 100644
--- a/man/dev_bern.Rd
+++ b/man/dev_bern.Rd
@@ -31,6 +31,7 @@ Other dev_dist:
 \code{\link[=dev_gamma]{dev_gamma()}},
 \code{\link[=dev_gamma_pois]{dev_gamma_pois()}},
 \code{\link[=dev_lnorm]{dev_lnorm()}},
+\code{\link[=dev_multinom]{dev_multinom()}},
 \code{\link[=dev_neg_binom]{dev_neg_binom()}},
 \code{\link[=dev_norm]{dev_norm()}},
 \code{\link[=dev_pois]{dev_pois()}},
diff --git a/man/dev_beta_binom.Rd b/man/dev_beta_binom.Rd
index 84ca558c..7edf8110 100644
--- a/man/dev_beta_binom.Rd
+++ b/man/dev_beta_binom.Rd
@@ -45,6 +45,7 @@ Other dev_dist:
 \code{\link[=dev_gamma]{dev_gamma()}},
 \code{\link[=dev_gamma_pois]{dev_gamma_pois()}},
 \code{\link[=dev_lnorm]{dev_lnorm()}},
+\code{\link[=dev_multinom]{dev_multinom()}},
 \code{\link[=dev_neg_binom]{dev_neg_binom()}},
 \code{\link[=dev_norm]{dev_norm()}},
 \code{\link[=dev_pois]{dev_pois()}},
diff --git a/man/dev_binom.Rd b/man/dev_binom.Rd
index 4d83f953..063a09fb 100644
--- a/man/dev_binom.Rd
+++ b/man/dev_binom.Rd
@@ -33,6 +33,7 @@ Other dev_dist:
 \code{\link[=dev_gamma]{dev_gamma()}},
 \code{\link[=dev_gamma_pois]{dev_gamma_pois()}},
 \code{\link[=dev_lnorm]{dev_lnorm()}},
+\code{\link[=dev_multinom]{dev_multinom()}},
 \code{\link[=dev_neg_binom]{dev_neg_binom()}},
 \code{\link[=dev_norm]{dev_norm()}},
 \code{\link[=dev_pois]{dev_pois()}},
diff --git a/man/dev_gamma.Rd b/man/dev_gamma.Rd
index 8041ccfe..cf9ebf4d 100644
--- a/man/dev_gamma.Rd
+++ b/man/dev_gamma.Rd
@@ -32,6 +32,7 @@ Other dev_dist:
 \code{\link[=dev_binom]{dev_binom()}},
 \code{\link[=dev_gamma_pois]{dev_gamma_pois()}},
 \code{\link[=dev_lnorm]{dev_lnorm()}},
+\code{\link[=dev_multinom]{dev_multinom()}},
 \code{\link[=dev_neg_binom]{dev_neg_binom()}},
 \code{\link[=dev_norm]{dev_norm()}},
 \code{\link[=dev_pois]{dev_pois()}},
diff --git a/man/dev_gamma_pois.Rd b/man/dev_gamma_pois.Rd
index 89f97a4f..1f68251b 100644
--- a/man/dev_gamma_pois.Rd
+++ b/man/dev_gamma_pois.Rd
@@ -33,6 +33,7 @@ Other dev_dist:
 \code{\link[=dev_binom]{dev_binom()}},
 \code{\link[=dev_gamma]{dev_gamma()}},
 \code{\link[=dev_lnorm]{dev_lnorm()}},
+\code{\link[=dev_multinom]{dev_multinom()}},
 \code{\link[=dev_neg_binom]{dev_neg_binom()}},
 \code{\link[=dev_norm]{dev_norm()}},
 \code{\link[=dev_pois]{dev_pois()}},
diff --git a/man/dev_lnorm.Rd b/man/dev_lnorm.Rd
index 0f10c26c..6512b9d3 100644
--- a/man/dev_lnorm.Rd
+++ b/man/dev_lnorm.Rd
@@ -33,6 +33,7 @@ Other dev_dist:
 \code{\link[=dev_binom]{dev_binom()}},
 \code{\link[=dev_gamma]{dev_gamma()}},
 \code{\link[=dev_gamma_pois]{dev_gamma_pois()}},
+\code{\link[=dev_multinom]{dev_multinom()}},
 \code{\link[=dev_neg_binom]{dev_neg_binom()}},
 \code{\link[=dev_norm]{dev_norm()}},
 \code{\link[=dev_pois]{dev_pois()}},
diff --git a/man/dev_multinom.Rd b/man/dev_multinom.Rd
new file mode 100644
index 00000000..57dd2172
--- /dev/null
+++ b/man/dev_multinom.Rd
@@ -0,0 +1,61 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/dev.R
+\name{dev_multinom}
+\alias{dev_multinom}
+\title{Multinomial Deviances}
+\usage{
+dev_multinom(x, size = 1, prob, group, res = FALSE)
+}
+\arguments{
+\item{x}{A non-negative whole numeric vector of the category counts.}
+
+\item{size}{A non-negative whole numeric vector of the number of trials.}
+
+\item{prob}{A numeric vector of the probability of the category. Must sum
+to 1 across the rows sharing the same \code{group}.}
+
+\item{group}{A vector identifying which rows belong to the same
+multinomial trial (whose \code{x} values sum to \code{size} and \code{prob} values sum
+to 1). Every group must have at least 2 rows and the same number of
+rows as the rest of the data (a fixed set of categories, as in
+multinomial logistic regression), and must not contain \code{NA}.}
+
+\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{
+Models the counts across two or more mutually exclusive categories from a
+fixed number of trials, in \emph{long} format: one row per category per
+trial, with \code{group} identifying which rows belong to the same trial.
+}
+\details{
+A category's deviance depends only on its own \code{x} and \code{mu = size * prob},
+not on the rest of its trial, so \code{group} is used only to validate \code{size}
+and \code{prob} (see \code{\link[=log_lik_multinom]{log_lik_multinom()}}), not in the calculation itself.
+\code{dev_multinom()} is the Poisson-equivalent deviance (see \code{\link[=dev_pois]{dev_pois()}}):
+summing it over a trial's rows recovers the trial's exact multinomial
+deviance.
+}
+\examples{
+dev_multinom(c(1, 3, 6), size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1))
+}
+\seealso{
+Other dev_dist:
+\code{\link[=dev_bern]{dev_bern()}},
+\code{\link[=dev_beta_binom]{dev_beta_binom()}},
+\code{\link[=dev_binom]{dev_binom()}},
+\code{\link[=dev_gamma]{dev_gamma()}},
+\code{\link[=dev_gamma_pois]{dev_gamma_pois()}},
+\code{\link[=dev_lnorm]{dev_lnorm()}},
+\code{\link[=dev_neg_binom]{dev_neg_binom()}},
+\code{\link[=dev_norm]{dev_norm()}},
+\code{\link[=dev_pois]{dev_pois()}},
+\code{\link[=dev_pois_zi]{dev_pois_zi()}},
+\code{\link[=dev_skewlnorm]{dev_skewlnorm()}},
+\code{\link[=dev_skewnorm]{dev_skewnorm()}},
+\code{\link[=dev_student]{dev_student()}}
+}
+\concept{dev_dist}
diff --git a/man/dev_neg_binom.Rd b/man/dev_neg_binom.Rd
index 5696099c..947dcf15 100644
--- a/man/dev_neg_binom.Rd
+++ b/man/dev_neg_binom.Rd
@@ -34,6 +34,7 @@ Other dev_dist:
 \code{\link[=dev_gamma]{dev_gamma()}},
 \code{\link[=dev_gamma_pois]{dev_gamma_pois()}},
 \code{\link[=dev_lnorm]{dev_lnorm()}},
+\code{\link[=dev_multinom]{dev_multinom()}},
 \code{\link[=dev_norm]{dev_norm()}},
 \code{\link[=dev_pois]{dev_pois()}},
 \code{\link[=dev_pois_zi]{dev_pois_zi()}},
diff --git a/man/dev_norm.Rd b/man/dev_norm.Rd
index 5b556302..fc685b14 100644
--- a/man/dev_norm.Rd
+++ b/man/dev_norm.Rd
@@ -33,6 +33,7 @@ Other dev_dist:
 \code{\link[=dev_gamma]{dev_gamma()}},
 \code{\link[=dev_gamma_pois]{dev_gamma_pois()}},
 \code{\link[=dev_lnorm]{dev_lnorm()}},
+\code{\link[=dev_multinom]{dev_multinom()}},
 \code{\link[=dev_neg_binom]{dev_neg_binom()}},
 \code{\link[=dev_pois]{dev_pois()}},
 \code{\link[=dev_pois_zi]{dev_pois_zi()}},
diff --git a/man/dev_pois.Rd b/man/dev_pois.Rd
index b99db9cd..f00cd65a 100644
--- a/man/dev_pois.Rd
+++ b/man/dev_pois.Rd
@@ -31,6 +31,7 @@ Other dev_dist:
 \code{\link[=dev_gamma]{dev_gamma()}},
 \code{\link[=dev_gamma_pois]{dev_gamma_pois()}},
 \code{\link[=dev_lnorm]{dev_lnorm()}},
+\code{\link[=dev_multinom]{dev_multinom()}},
 \code{\link[=dev_neg_binom]{dev_neg_binom()}},
 \code{\link[=dev_norm]{dev_norm()}},
 \code{\link[=dev_pois_zi]{dev_pois_zi()}},
diff --git a/man/dev_pois_zi.Rd b/man/dev_pois_zi.Rd
index 23db2ccf..7aa37c09 100644
--- a/man/dev_pois_zi.Rd
+++ b/man/dev_pois_zi.Rd
@@ -34,6 +34,7 @@ Other dev_dist:
 \code{\link[=dev_gamma]{dev_gamma()}},
 \code{\link[=dev_gamma_pois]{dev_gamma_pois()}},
 \code{\link[=dev_lnorm]{dev_lnorm()}},
+\code{\link[=dev_multinom]{dev_multinom()}},
 \code{\link[=dev_neg_binom]{dev_neg_binom()}},
 \code{\link[=dev_norm]{dev_norm()}},
 \code{\link[=dev_pois]{dev_pois()}},
diff --git a/man/dev_skewlnorm.Rd b/man/dev_skewlnorm.Rd
index 5b39cef8..c3148042 100644
--- a/man/dev_skewlnorm.Rd
+++ b/man/dev_skewlnorm.Rd
@@ -40,6 +40,7 @@ Other dev_dist:
 \code{\link[=dev_gamma]{dev_gamma()}},
 \code{\link[=dev_gamma_pois]{dev_gamma_pois()}},
 \code{\link[=dev_lnorm]{dev_lnorm()}},
+\code{\link[=dev_multinom]{dev_multinom()}},
 \code{\link[=dev_neg_binom]{dev_neg_binom()}},
 \code{\link[=dev_norm]{dev_norm()}},
 \code{\link[=dev_pois]{dev_pois()}},
diff --git a/man/dev_skewnorm.Rd b/man/dev_skewnorm.Rd
index 389cea5d..c46bd9b7 100644
--- a/man/dev_skewnorm.Rd
+++ b/man/dev_skewnorm.Rd
@@ -39,6 +39,7 @@ Other dev_dist:
 \code{\link[=dev_gamma]{dev_gamma()}},
 \code{\link[=dev_gamma_pois]{dev_gamma_pois()}},
 \code{\link[=dev_lnorm]{dev_lnorm()}},
+\code{\link[=dev_multinom]{dev_multinom()}},
 \code{\link[=dev_neg_binom]{dev_neg_binom()}},
 \code{\link[=dev_norm]{dev_norm()}},
 \code{\link[=dev_pois]{dev_pois()}},
diff --git a/man/dev_student.Rd b/man/dev_student.Rd
index 3f9739d6..6ce0ec43 100644
--- a/man/dev_student.Rd
+++ b/man/dev_student.Rd
@@ -36,6 +36,7 @@ Other dev_dist:
 \code{\link[=dev_gamma]{dev_gamma()}},
 \code{\link[=dev_gamma_pois]{dev_gamma_pois()}},
 \code{\link[=dev_lnorm]{dev_lnorm()}},
+\code{\link[=dev_multinom]{dev_multinom()}},
 \code{\link[=dev_neg_binom]{dev_neg_binom()}},
 \code{\link[=dev_norm]{dev_norm()}},
 \code{\link[=dev_pois]{dev_pois()}},
diff --git a/man/log_lik_bern.Rd b/man/log_lik_bern.Rd
index 18ddb194..9212e2ae 100644
--- a/man/log_lik_bern.Rd
+++ b/man/log_lik_bern.Rd
@@ -31,6 +31,7 @@ Other log_lik_dist:
 \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}},
 \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}},
 \code{\link[=log_lik_lnorm]{log_lik_lnorm()}},
+\code{\link[=log_lik_multinom]{log_lik_multinom()}},
 \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}},
 \code{\link[=log_lik_norm]{log_lik_norm()}},
 \code{\link[=log_lik_pois]{log_lik_pois()}},
diff --git a/man/log_lik_beta.Rd b/man/log_lik_beta.Rd
index 9e6453fb..d4def7b0 100644
--- a/man/log_lik_beta.Rd
+++ b/man/log_lik_beta.Rd
@@ -36,6 +36,7 @@ Other log_lik_dist:
 \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}},
 \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}},
 \code{\link[=log_lik_lnorm]{log_lik_lnorm()}},
+\code{\link[=log_lik_multinom]{log_lik_multinom()}},
 \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}},
 \code{\link[=log_lik_norm]{log_lik_norm()}},
 \code{\link[=log_lik_pois]{log_lik_pois()}},
diff --git a/man/log_lik_beta_binom.Rd b/man/log_lik_beta_binom.Rd
index 8c737976..8528b741 100644
--- a/man/log_lik_beta_binom.Rd
+++ b/man/log_lik_beta_binom.Rd
@@ -58,6 +58,7 @@ Other log_lik_dist:
 \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}},
 \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}},
 \code{\link[=log_lik_lnorm]{log_lik_lnorm()}},
+\code{\link[=log_lik_multinom]{log_lik_multinom()}},
 \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}},
 \code{\link[=log_lik_norm]{log_lik_norm()}},
 \code{\link[=log_lik_pois]{log_lik_pois()}},
diff --git a/man/log_lik_binom.Rd b/man/log_lik_binom.Rd
index 475d74f8..5e90cf80 100644
--- a/man/log_lik_binom.Rd
+++ b/man/log_lik_binom.Rd
@@ -37,6 +37,7 @@ Other log_lik_dist:
 \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}},
 \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}},
 \code{\link[=log_lik_lnorm]{log_lik_lnorm()}},
+\code{\link[=log_lik_multinom]{log_lik_multinom()}},
 \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}},
 \code{\link[=log_lik_norm]{log_lik_norm()}},
 \code{\link[=log_lik_pois]{log_lik_pois()}},
diff --git a/man/log_lik_exp.Rd b/man/log_lik_exp.Rd
index 5192ae8f..94f16a05 100644
--- a/man/log_lik_exp.Rd
+++ b/man/log_lik_exp.Rd
@@ -34,6 +34,7 @@ Other log_lik_dist:
 \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}},
 \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}},
 \code{\link[=log_lik_lnorm]{log_lik_lnorm()}},
+\code{\link[=log_lik_multinom]{log_lik_multinom()}},
 \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}},
 \code{\link[=log_lik_norm]{log_lik_norm()}},
 \code{\link[=log_lik_pois]{log_lik_pois()}},
diff --git a/man/log_lik_gamma.Rd b/man/log_lik_gamma.Rd
index 953dfd9e..9148ff42 100644
--- a/man/log_lik_gamma.Rd
+++ b/man/log_lik_gamma.Rd
@@ -36,6 +36,7 @@ Other log_lik_dist:
 \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}},
 \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}},
 \code{\link[=log_lik_lnorm]{log_lik_lnorm()}},
+\code{\link[=log_lik_multinom]{log_lik_multinom()}},
 \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}},
 \code{\link[=log_lik_norm]{log_lik_norm()}},
 \code{\link[=log_lik_pois]{log_lik_pois()}},
diff --git a/man/log_lik_gamma_pois.Rd b/man/log_lik_gamma_pois.Rd
index 51f94e71..47e1c63b 100644
--- a/man/log_lik_gamma_pois.Rd
+++ b/man/log_lik_gamma_pois.Rd
@@ -37,6 +37,7 @@ Other log_lik_dist:
 \code{\link[=log_lik_gamma]{log_lik_gamma()}},
 \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}},
 \code{\link[=log_lik_lnorm]{log_lik_lnorm()}},
+\code{\link[=log_lik_multinom]{log_lik_multinom()}},
 \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}},
 \code{\link[=log_lik_norm]{log_lik_norm()}},
 \code{\link[=log_lik_pois]{log_lik_pois()}},
diff --git a/man/log_lik_gamma_pois_zi.Rd b/man/log_lik_gamma_pois_zi.Rd
index 1700d074..1f8e102d 100644
--- a/man/log_lik_gamma_pois_zi.Rd
+++ b/man/log_lik_gamma_pois_zi.Rd
@@ -47,6 +47,7 @@ Other log_lik_dist:
 \code{\link[=log_lik_gamma]{log_lik_gamma()}},
 \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}},
 \code{\link[=log_lik_lnorm]{log_lik_lnorm()}},
+\code{\link[=log_lik_multinom]{log_lik_multinom()}},
 \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}},
 \code{\link[=log_lik_norm]{log_lik_norm()}},
 \code{\link[=log_lik_pois]{log_lik_pois()}},
diff --git a/man/log_lik_lnorm.Rd b/man/log_lik_lnorm.Rd
index d02d6233..ac7e102f 100644
--- a/man/log_lik_lnorm.Rd
+++ b/man/log_lik_lnorm.Rd
@@ -37,6 +37,7 @@ Other log_lik_dist:
 \code{\link[=log_lik_gamma]{log_lik_gamma()}},
 \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}},
 \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}},
+\code{\link[=log_lik_multinom]{log_lik_multinom()}},
 \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}},
 \code{\link[=log_lik_norm]{log_lik_norm()}},
 \code{\link[=log_lik_pois]{log_lik_pois()}},
diff --git a/man/log_lik_multinom.Rd b/man/log_lik_multinom.Rd
new file mode 100644
index 00000000..f7b1d8f5
--- /dev/null
+++ b/man/log_lik_multinom.Rd
@@ -0,0 +1,67 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/log-lik.R
+\name{log_lik_multinom}
+\alias{log_lik_multinom}
+\title{Multinomial Log-Likelihood}
+\usage{
+log_lik_multinom(x, size = 1, prob, group)
+}
+\arguments{
+\item{x}{A non-negative whole numeric vector of the category counts.}
+
+\item{size}{A non-negative whole numeric vector of the number of trials.}
+
+\item{prob}{A numeric vector of the probability of the category. Must sum
+to 1 across the rows sharing the same \code{group}. \code{NA} in \code{size} or \code{prob}
+for any row of a trial makes the log-likelihood \code{NA} for every row of
+that trial, since a trial's categories are scored jointly.}
+
+\item{group}{A vector identifying which rows belong to the same
+multinomial trial (whose \code{x} values sum to \code{size} and \code{prob} values sum
+to 1). Every group must have at least 2 rows and the same number of
+rows as the rest of the data (a fixed set of categories, as in
+multinomial logistic regression), and must not contain \code{NA}.}
+}
+\value{
+An numeric vector of the corresponding log-likelihoods, one value
+per row of \code{x}.
+}
+\description{
+Models the counts across two or more mutually exclusive categories from a
+fixed number of trials, in \emph{long} format: one row per category per
+trial, with \code{group} identifying which rows belong to the same trial. All
+rows sharing a \code{group} must have the same \code{size}, and their \code{prob} values
+must sum to 1.
+}
+\details{
+A trial's log-likelihood doesn't split evenly across its rows, since the
+multinomial coefficient belongs to the whole trial. \code{log_lik_multinom()}
+uses the multinomial-as-independent-Poissons identity: each row's value
+is the Poisson log-likelihood of \code{x} given \code{mu = size * prob}, minus an
+even share of the trial's normalizing constant, so summing over a
+\code{group} recovers the trial's exact multinomial log-likelihood.
+}
+\examples{
+log_lik_multinom(c(1, 3, 6), size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1))
+}
+\seealso{
+Other log_lik_dist:
+\code{\link[=log_lik_bern]{log_lik_bern()}},
+\code{\link[=log_lik_beta]{log_lik_beta()}},
+\code{\link[=log_lik_beta_binom]{log_lik_beta_binom()}},
+\code{\link[=log_lik_binom]{log_lik_binom()}},
+\code{\link[=log_lik_exp]{log_lik_exp()}},
+\code{\link[=log_lik_gamma]{log_lik_gamma()}},
+\code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}},
+\code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}},
+\code{\link[=log_lik_lnorm]{log_lik_lnorm()}},
+\code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}},
+\code{\link[=log_lik_norm]{log_lik_norm()}},
+\code{\link[=log_lik_pois]{log_lik_pois()}},
+\code{\link[=log_lik_pois_zi]{log_lik_pois_zi()}},
+\code{\link[=log_lik_skewlnorm]{log_lik_skewlnorm()}},
+\code{\link[=log_lik_skewnorm]{log_lik_skewnorm()}},
+\code{\link[=log_lik_student]{log_lik_student()}},
+\code{\link[=log_lik_unif]{log_lik_unif()}}
+}
+\concept{log_lik_dist}
diff --git a/man/log_lik_neg_binom.Rd b/man/log_lik_neg_binom.Rd
index a4e7609c..c06b5259 100644
--- a/man/log_lik_neg_binom.Rd
+++ b/man/log_lik_neg_binom.Rd
@@ -38,6 +38,7 @@ Other log_lik_dist:
 \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}},
 \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}},
 \code{\link[=log_lik_lnorm]{log_lik_lnorm()}},
+\code{\link[=log_lik_multinom]{log_lik_multinom()}},
 \code{\link[=log_lik_norm]{log_lik_norm()}},
 \code{\link[=log_lik_pois]{log_lik_pois()}},
 \code{\link[=log_lik_pois_zi]{log_lik_pois_zi()}},
diff --git a/man/log_lik_norm.Rd b/man/log_lik_norm.Rd
index 8018e102..feed9c0a 100644
--- a/man/log_lik_norm.Rd
+++ b/man/log_lik_norm.Rd
@@ -37,6 +37,7 @@ Other log_lik_dist:
 \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}},
 \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}},
 \code{\link[=log_lik_lnorm]{log_lik_lnorm()}},
+\code{\link[=log_lik_multinom]{log_lik_multinom()}},
 \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}},
 \code{\link[=log_lik_pois]{log_lik_pois()}},
 \code{\link[=log_lik_pois_zi]{log_lik_pois_zi()}},
diff --git a/man/log_lik_pois.Rd b/man/log_lik_pois.Rd
index 0eb155e4..421c67ff 100644
--- a/man/log_lik_pois.Rd
+++ b/man/log_lik_pois.Rd
@@ -35,6 +35,7 @@ Other log_lik_dist:
 \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}},
 \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}},
 \code{\link[=log_lik_lnorm]{log_lik_lnorm()}},
+\code{\link[=log_lik_multinom]{log_lik_multinom()}},
 \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}},
 \code{\link[=log_lik_norm]{log_lik_norm()}},
 \code{\link[=log_lik_pois_zi]{log_lik_pois_zi()}},
diff --git a/man/log_lik_pois_zi.Rd b/man/log_lik_pois_zi.Rd
index b1e67d4b..8138a82f 100644
--- a/man/log_lik_pois_zi.Rd
+++ b/man/log_lik_pois_zi.Rd
@@ -38,6 +38,7 @@ Other log_lik_dist:
 \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}},
 \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}},
 \code{\link[=log_lik_lnorm]{log_lik_lnorm()}},
+\code{\link[=log_lik_multinom]{log_lik_multinom()}},
 \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}},
 \code{\link[=log_lik_norm]{log_lik_norm()}},
 \code{\link[=log_lik_pois]{log_lik_pois()}},
diff --git a/man/log_lik_skewlnorm.Rd b/man/log_lik_skewlnorm.Rd
index 99660b49..d9a9fbbb 100644
--- a/man/log_lik_skewlnorm.Rd
+++ b/man/log_lik_skewlnorm.Rd
@@ -51,6 +51,7 @@ Other log_lik_dist:
 \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}},
 \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}},
 \code{\link[=log_lik_lnorm]{log_lik_lnorm()}},
+\code{\link[=log_lik_multinom]{log_lik_multinom()}},
 \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}},
 \code{\link[=log_lik_norm]{log_lik_norm()}},
 \code{\link[=log_lik_pois]{log_lik_pois()}},
diff --git a/man/log_lik_skewnorm.Rd b/man/log_lik_skewnorm.Rd
index f4acf7f7..ce0a9f2a 100644
--- a/man/log_lik_skewnorm.Rd
+++ b/man/log_lik_skewnorm.Rd
@@ -43,6 +43,7 @@ Other log_lik_dist:
 \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}},
 \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}},
 \code{\link[=log_lik_lnorm]{log_lik_lnorm()}},
+\code{\link[=log_lik_multinom]{log_lik_multinom()}},
 \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}},
 \code{\link[=log_lik_norm]{log_lik_norm()}},
 \code{\link[=log_lik_pois]{log_lik_pois()}},
diff --git a/man/log_lik_student.Rd b/man/log_lik_student.Rd
index 31e83fa1..d3700d73 100644
--- a/man/log_lik_student.Rd
+++ b/man/log_lik_student.Rd
@@ -40,6 +40,7 @@ Other log_lik_dist:
 \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}},
 \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}},
 \code{\link[=log_lik_lnorm]{log_lik_lnorm()}},
+\code{\link[=log_lik_multinom]{log_lik_multinom()}},
 \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}},
 \code{\link[=log_lik_norm]{log_lik_norm()}},
 \code{\link[=log_lik_pois]{log_lik_pois()}},
diff --git a/man/log_lik_unif.Rd b/man/log_lik_unif.Rd
index e9eef960..90b4fd54 100644
--- a/man/log_lik_unif.Rd
+++ b/man/log_lik_unif.Rd
@@ -33,6 +33,7 @@ Other log_lik_dist:
 \code{\link[=log_lik_gamma_pois]{log_lik_gamma_pois()}},
 \code{\link[=log_lik_gamma_pois_zi]{log_lik_gamma_pois_zi()}},
 \code{\link[=log_lik_lnorm]{log_lik_lnorm()}},
+\code{\link[=log_lik_multinom]{log_lik_multinom()}},
 \code{\link[=log_lik_neg_binom]{log_lik_neg_binom()}},
 \code{\link[=log_lik_norm]{log_lik_norm()}},
 \code{\link[=log_lik_pois]{log_lik_pois()}},
diff --git a/man/params.Rd b/man/params.Rd
index 47b31cdf..4105f9d4 100644
--- a/man/params.Rd
+++ b/man/params.Rd
@@ -19,6 +19,12 @@ level.}
 \item{directional}{A flag specifying whether probabilities less than 0.5
 should be returned as negative values.}
 
+\item{group}{A vector identifying which rows belong to the same
+multinomial trial (whose \code{x} values sum to \code{size} and \code{prob} values sum
+to 1). Every group must have at least 2 rows and the same number of
+rows as the rest of the data (a fixed set of categories, as in
+multinomial logistic regression), and must not contain \code{NA}.}
+
 \item{lambda}{A non-negative numeric vector of means.}
 
 \item{level}{A number > 0 and <= 1 specifying the probability coverage of the
diff --git a/man/ran_bern.Rd b/man/ran_bern.Rd
index bde299dd..20e1613b 100644
--- a/man/ran_bern.Rd
+++ b/man/ran_bern.Rd
@@ -30,6 +30,7 @@ Other ran_dist:
 \code{\link[=ran_gamma_pois]{ran_gamma_pois()}},
 \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}},
 \code{\link[=ran_lnorm]{ran_lnorm()}},
+\code{\link[=ran_multinom]{ran_multinom()}},
 \code{\link[=ran_neg_binom]{ran_neg_binom()}},
 \code{\link[=ran_norm]{ran_norm()}},
 \code{\link[=ran_pois]{ran_pois()}},
diff --git a/man/ran_beta_binom.Rd b/man/ran_beta_binom.Rd
index daa14527..da22edbb 100644
--- a/man/ran_beta_binom.Rd
+++ b/man/ran_beta_binom.Rd
@@ -43,6 +43,7 @@ Other ran_dist:
 \code{\link[=ran_gamma_pois]{ran_gamma_pois()}},
 \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}},
 \code{\link[=ran_lnorm]{ran_lnorm()}},
+\code{\link[=ran_multinom]{ran_multinom()}},
 \code{\link[=ran_neg_binom]{ran_neg_binom()}},
 \code{\link[=ran_norm]{ran_norm()}},
 \code{\link[=ran_pois]{ran_pois()}},
diff --git a/man/ran_binom.Rd b/man/ran_binom.Rd
index 3ec3d020..a352e434 100644
--- a/man/ran_binom.Rd
+++ b/man/ran_binom.Rd
@@ -32,6 +32,7 @@ Other ran_dist:
 \code{\link[=ran_gamma_pois]{ran_gamma_pois()}},
 \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}},
 \code{\link[=ran_lnorm]{ran_lnorm()}},
+\code{\link[=ran_multinom]{ran_multinom()}},
 \code{\link[=ran_neg_binom]{ran_neg_binom()}},
 \code{\link[=ran_norm]{ran_norm()}},
 \code{\link[=ran_pois]{ran_pois()}},
diff --git a/man/ran_gamma.Rd b/man/ran_gamma.Rd
index c348912a..9dc857c7 100644
--- a/man/ran_gamma.Rd
+++ b/man/ran_gamma.Rd
@@ -31,6 +31,7 @@ Other ran_dist:
 \code{\link[=ran_gamma_pois]{ran_gamma_pois()}},
 \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}},
 \code{\link[=ran_lnorm]{ran_lnorm()}},
+\code{\link[=ran_multinom]{ran_multinom()}},
 \code{\link[=ran_neg_binom]{ran_neg_binom()}},
 \code{\link[=ran_norm]{ran_norm()}},
 \code{\link[=ran_pois]{ran_pois()}},
diff --git a/man/ran_gamma_pois.Rd b/man/ran_gamma_pois.Rd
index 83925a90..342f4c15 100644
--- a/man/ran_gamma_pois.Rd
+++ b/man/ran_gamma_pois.Rd
@@ -32,6 +32,7 @@ Other ran_dist:
 \code{\link[=ran_gamma]{ran_gamma()}},
 \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}},
 \code{\link[=ran_lnorm]{ran_lnorm()}},
+\code{\link[=ran_multinom]{ran_multinom()}},
 \code{\link[=ran_neg_binom]{ran_neg_binom()}},
 \code{\link[=ran_norm]{ran_norm()}},
 \code{\link[=ran_pois]{ran_pois()}},
diff --git a/man/ran_gamma_pois_zi.Rd b/man/ran_gamma_pois_zi.Rd
index 44240a73..ba438271 100644
--- a/man/ran_gamma_pois_zi.Rd
+++ b/man/ran_gamma_pois_zi.Rd
@@ -35,6 +35,7 @@ Other ran_dist:
 \code{\link[=ran_gamma]{ran_gamma()}},
 \code{\link[=ran_gamma_pois]{ran_gamma_pois()}},
 \code{\link[=ran_lnorm]{ran_lnorm()}},
+\code{\link[=ran_multinom]{ran_multinom()}},
 \code{\link[=ran_neg_binom]{ran_neg_binom()}},
 \code{\link[=ran_norm]{ran_norm()}},
 \code{\link[=ran_pois]{ran_pois()}},
diff --git a/man/ran_lnorm.Rd b/man/ran_lnorm.Rd
index e852d444..58b7bb2e 100644
--- a/man/ran_lnorm.Rd
+++ b/man/ran_lnorm.Rd
@@ -32,6 +32,7 @@ Other ran_dist:
 \code{\link[=ran_gamma]{ran_gamma()}},
 \code{\link[=ran_gamma_pois]{ran_gamma_pois()}},
 \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}},
+\code{\link[=ran_multinom]{ran_multinom()}},
 \code{\link[=ran_neg_binom]{ran_neg_binom()}},
 \code{\link[=ran_norm]{ran_norm()}},
 \code{\link[=ran_pois]{ran_pois()}},
diff --git a/man/ran_multinom.Rd b/man/ran_multinom.Rd
new file mode 100644
index 00000000..de3bcc04
--- /dev/null
+++ b/man/ran_multinom.Rd
@@ -0,0 +1,59 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ran.R
+\name{ran_multinom}
+\alias{ran_multinom}
+\title{Multinomial Random Samples}
+\usage{
+ran_multinom(size = 1, prob, group)
+}
+\arguments{
+\item{size}{A non-negative whole numeric vector of the number of trials.}
+
+\item{prob}{A numeric vector of the probability of the category. Must sum
+to 1 across the rows sharing the same \code{group}. \code{NA} in \code{size} or \code{prob}
+for any row of a trial makes the sample \code{NA} for every row of that
+trial, since a trial's categories are drawn jointly.}
+
+\item{group}{A vector identifying which rows belong to the same
+multinomial trial (whose \code{x} values sum to \code{size} and \code{prob} values sum
+to 1). Every group must have at least 2 rows and the same number of
+rows as the rest of the data (a fixed set of categories, as in
+multinomial logistic regression), and must not contain \code{NA}.}
+}
+\value{
+An integer vector of the random samples, one per row of \code{prob}.
+}
+\description{
+Models the counts across two or more mutually exclusive categories from a
+fixed number of trials, in \emph{long} format: one value per category per
+trial, with \code{group} identifying which rows belong to the same trial. All
+rows sharing a \code{group} must have the same \code{size}, and their \code{prob} values
+must sum to 1.
+}
+\details{
+Unlike the other \verb{ran_*()} functions, \code{ran_multinom()} has no \code{n}
+argument: the number of samples is fully determined by \code{length(prob)}
+(equivalently \code{length(group)}), since a trial's categories can't be
+generated independently of one another.
+}
+\examples{
+ran_multinom(size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1))
+}
+\seealso{
+Other ran_dist:
+\code{\link[=ran_bern]{ran_bern()}},
+\code{\link[=ran_beta_binom]{ran_beta_binom()}},
+\code{\link[=ran_binom]{ran_binom()}},
+\code{\link[=ran_gamma]{ran_gamma()}},
+\code{\link[=ran_gamma_pois]{ran_gamma_pois()}},
+\code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}},
+\code{\link[=ran_lnorm]{ran_lnorm()}},
+\code{\link[=ran_neg_binom]{ran_neg_binom()}},
+\code{\link[=ran_norm]{ran_norm()}},
+\code{\link[=ran_pois]{ran_pois()}},
+\code{\link[=ran_pois_zi]{ran_pois_zi()}},
+\code{\link[=ran_skewlnorm]{ran_skewlnorm()}},
+\code{\link[=ran_skewnorm]{ran_skewnorm()}},
+\code{\link[=ran_student]{ran_student()}}
+}
+\concept{ran_dist}
diff --git a/man/ran_neg_binom.Rd b/man/ran_neg_binom.Rd
index b70d571e..5c66caa5 100644
--- a/man/ran_neg_binom.Rd
+++ b/man/ran_neg_binom.Rd
@@ -33,6 +33,7 @@ Other ran_dist:
 \code{\link[=ran_gamma_pois]{ran_gamma_pois()}},
 \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}},
 \code{\link[=ran_lnorm]{ran_lnorm()}},
+\code{\link[=ran_multinom]{ran_multinom()}},
 \code{\link[=ran_norm]{ran_norm()}},
 \code{\link[=ran_pois]{ran_pois()}},
 \code{\link[=ran_pois_zi]{ran_pois_zi()}},
diff --git a/man/ran_norm.Rd b/man/ran_norm.Rd
index 84564c6b..1b92e336 100644
--- a/man/ran_norm.Rd
+++ b/man/ran_norm.Rd
@@ -32,6 +32,7 @@ Other ran_dist:
 \code{\link[=ran_gamma_pois]{ran_gamma_pois()}},
 \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}},
 \code{\link[=ran_lnorm]{ran_lnorm()}},
+\code{\link[=ran_multinom]{ran_multinom()}},
 \code{\link[=ran_neg_binom]{ran_neg_binom()}},
 \code{\link[=ran_pois]{ran_pois()}},
 \code{\link[=ran_pois_zi]{ran_pois_zi()}},
diff --git a/man/ran_pois.Rd b/man/ran_pois.Rd
index 3afe4a9c..e22a041d 100644
--- a/man/ran_pois.Rd
+++ b/man/ran_pois.Rd
@@ -30,6 +30,7 @@ Other ran_dist:
 \code{\link[=ran_gamma_pois]{ran_gamma_pois()}},
 \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}},
 \code{\link[=ran_lnorm]{ran_lnorm()}},
+\code{\link[=ran_multinom]{ran_multinom()}},
 \code{\link[=ran_neg_binom]{ran_neg_binom()}},
 \code{\link[=ran_norm]{ran_norm()}},
 \code{\link[=ran_pois_zi]{ran_pois_zi()}},
diff --git a/man/ran_pois_zi.Rd b/man/ran_pois_zi.Rd
index 9668cb61..509f5fc0 100644
--- a/man/ran_pois_zi.Rd
+++ b/man/ran_pois_zi.Rd
@@ -33,6 +33,7 @@ Other ran_dist:
 \code{\link[=ran_gamma_pois]{ran_gamma_pois()}},
 \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}},
 \code{\link[=ran_lnorm]{ran_lnorm()}},
+\code{\link[=ran_multinom]{ran_multinom()}},
 \code{\link[=ran_neg_binom]{ran_neg_binom()}},
 \code{\link[=ran_norm]{ran_norm()}},
 \code{\link[=ran_pois]{ran_pois()}},
diff --git a/man/ran_skewlnorm.Rd b/man/ran_skewlnorm.Rd
index 36f925d1..c02b39d4 100644
--- a/man/ran_skewlnorm.Rd
+++ b/man/ran_skewlnorm.Rd
@@ -39,6 +39,7 @@ Other ran_dist:
 \code{\link[=ran_gamma_pois]{ran_gamma_pois()}},
 \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}},
 \code{\link[=ran_lnorm]{ran_lnorm()}},
+\code{\link[=ran_multinom]{ran_multinom()}},
 \code{\link[=ran_neg_binom]{ran_neg_binom()}},
 \code{\link[=ran_norm]{ran_norm()}},
 \code{\link[=ran_pois]{ran_pois()}},
diff --git a/man/ran_skewnorm.Rd b/man/ran_skewnorm.Rd
index e8152d03..9b1a4375 100644
--- a/man/ran_skewnorm.Rd
+++ b/man/ran_skewnorm.Rd
@@ -38,6 +38,7 @@ Other ran_dist:
 \code{\link[=ran_gamma_pois]{ran_gamma_pois()}},
 \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}},
 \code{\link[=ran_lnorm]{ran_lnorm()}},
+\code{\link[=ran_multinom]{ran_multinom()}},
 \code{\link[=ran_neg_binom]{ran_neg_binom()}},
 \code{\link[=ran_norm]{ran_norm()}},
 \code{\link[=ran_pois]{ran_pois()}},
diff --git a/man/ran_student.Rd b/man/ran_student.Rd
index 8a883907..77fdde4c 100644
--- a/man/ran_student.Rd
+++ b/man/ran_student.Rd
@@ -35,6 +35,7 @@ Other ran_dist:
 \code{\link[=ran_gamma_pois]{ran_gamma_pois()}},
 \code{\link[=ran_gamma_pois_zi]{ran_gamma_pois_zi()}},
 \code{\link[=ran_lnorm]{ran_lnorm()}},
+\code{\link[=ran_multinom]{ran_multinom()}},
 \code{\link[=ran_neg_binom]{ran_neg_binom()}},
 \code{\link[=ran_norm]{ran_norm()}},
 \code{\link[=ran_pois]{ran_pois()}},
diff --git a/man/res_bern.Rd b/man/res_bern.Rd
index c5beb276..5da2d149 100644
--- a/man/res_bern.Rd
+++ b/man/res_bern.Rd
@@ -34,6 +34,7 @@ Other res_dist:
 \code{\link[=res_gamma_pois]{res_gamma_pois()}},
 \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}},
 \code{\link[=res_lnorm]{res_lnorm()}},
+\code{\link[=res_multinom]{res_multinom()}},
 \code{\link[=res_neg_binom]{res_neg_binom()}},
 \code{\link[=res_norm]{res_norm()}},
 \code{\link[=res_pois]{res_pois()}},
diff --git a/man/res_beta_binom.Rd b/man/res_beta_binom.Rd
index 479167cd..5452b575 100644
--- a/man/res_beta_binom.Rd
+++ b/man/res_beta_binom.Rd
@@ -54,6 +54,7 @@ Other res_dist:
 \code{\link[=res_gamma_pois]{res_gamma_pois()}},
 \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}},
 \code{\link[=res_lnorm]{res_lnorm()}},
+\code{\link[=res_multinom]{res_multinom()}},
 \code{\link[=res_neg_binom]{res_neg_binom()}},
 \code{\link[=res_norm]{res_norm()}},
 \code{\link[=res_pois]{res_pois()}},
diff --git a/man/res_binom.Rd b/man/res_binom.Rd
index b8026c4e..a00ed641 100644
--- a/man/res_binom.Rd
+++ b/man/res_binom.Rd
@@ -36,6 +36,7 @@ Other res_dist:
 \code{\link[=res_gamma_pois]{res_gamma_pois()}},
 \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}},
 \code{\link[=res_lnorm]{res_lnorm()}},
+\code{\link[=res_multinom]{res_multinom()}},
 \code{\link[=res_neg_binom]{res_neg_binom()}},
 \code{\link[=res_norm]{res_norm()}},
 \code{\link[=res_pois]{res_pois()}},
diff --git a/man/res_gamma.Rd b/man/res_gamma.Rd
index 58dff315..6902641e 100644
--- a/man/res_gamma.Rd
+++ b/man/res_gamma.Rd
@@ -35,6 +35,7 @@ Other res_dist:
 \code{\link[=res_gamma_pois]{res_gamma_pois()}},
 \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}},
 \code{\link[=res_lnorm]{res_lnorm()}},
+\code{\link[=res_multinom]{res_multinom()}},
 \code{\link[=res_neg_binom]{res_neg_binom()}},
 \code{\link[=res_norm]{res_norm()}},
 \code{\link[=res_pois]{res_pois()}},
diff --git a/man/res_gamma_pois.Rd b/man/res_gamma_pois.Rd
index ae824418..c2e7e481 100644
--- a/man/res_gamma_pois.Rd
+++ b/man/res_gamma_pois.Rd
@@ -36,6 +36,7 @@ Other res_dist:
 \code{\link[=res_gamma]{res_gamma()}},
 \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}},
 \code{\link[=res_lnorm]{res_lnorm()}},
+\code{\link[=res_multinom]{res_multinom()}},
 \code{\link[=res_neg_binom]{res_neg_binom()}},
 \code{\link[=res_norm]{res_norm()}},
 \code{\link[=res_pois]{res_pois()}},
diff --git a/man/res_gamma_pois_zi.Rd b/man/res_gamma_pois_zi.Rd
index 5c3080bd..3d4a80d2 100644
--- a/man/res_gamma_pois_zi.Rd
+++ b/man/res_gamma_pois_zi.Rd
@@ -45,6 +45,7 @@ Other res_dist:
 \code{\link[=res_gamma]{res_gamma()}},
 \code{\link[=res_gamma_pois]{res_gamma_pois()}},
 \code{\link[=res_lnorm]{res_lnorm()}},
+\code{\link[=res_multinom]{res_multinom()}},
 \code{\link[=res_neg_binom]{res_neg_binom()}},
 \code{\link[=res_norm]{res_norm()}},
 \code{\link[=res_pois]{res_pois()}},
diff --git a/man/res_lnorm.Rd b/man/res_lnorm.Rd
index 3fbfd1a5..1a674850 100644
--- a/man/res_lnorm.Rd
+++ b/man/res_lnorm.Rd
@@ -36,6 +36,7 @@ Other res_dist:
 \code{\link[=res_gamma]{res_gamma()}},
 \code{\link[=res_gamma_pois]{res_gamma_pois()}},
 \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}},
+\code{\link[=res_multinom]{res_multinom()}},
 \code{\link[=res_neg_binom]{res_neg_binom()}},
 \code{\link[=res_norm]{res_norm()}},
 \code{\link[=res_pois]{res_pois()}},
diff --git a/man/res_multinom.Rd b/man/res_multinom.Rd
new file mode 100644
index 00000000..ce90ef38
--- /dev/null
+++ b/man/res_multinom.Rd
@@ -0,0 +1,69 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/res.R
+\name{res_multinom}
+\alias{res_multinom}
+\title{Multinomial Residuals}
+\usage{
+res_multinom(x, size = 1, prob, group, type = "dev", simulate = FALSE)
+}
+\arguments{
+\item{x}{A non-negative whole numeric vector of the category counts.}
+
+\item{size}{A non-negative whole numeric vector of the number of trials.}
+
+\item{prob}{A numeric vector of the probability of the category. Must sum
+to 1 across the rows sharing the same \code{group}.}
+
+\item{group}{A vector identifying which rows belong to the same
+multinomial trial (whose \code{x} values sum to \code{size} and \code{prob} values sum
+to 1). Every group must have at least 2 rows and the same number of
+rows as the rest of the data (a fixed set of categories, as in
+multinomial logistic regression), and must not contain \code{NA}.}
+
+\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{
+Models the counts across two or more mutually exclusive categories from a
+fixed number of trials, in \emph{long} format: one row per category per
+trial, with \code{group} identifying which rows belong to the same trial (see
+\code{\link[=log_lik_multinom]{log_lik_multinom()}} for details). \code{res_multinom()} returns one residual
+per row, not one per trial, since a trial's categories aren't
+independent and so have no single meaningful residual as a whole; the
+classic per-trial deviance statistic can be recovered by summing the
+squared \code{type = "dev"} residuals within a \code{group}.
+}
+\details{
+\code{group} is validated (same \code{size}, \code{prob} summing to 1, no singleton or
+short groups, no \code{NA}) regardless of \code{simulate}, but is only otherwise
+used when \code{simulate = TRUE}, to draw a joint, correlation-preserving
+replicate per trial (via \code{\link[=ran_multinom]{ran_multinom()}}) rather than simulating each
+category independently, which requires \code{res_multinom()} to see every row
+of a \code{group} in the same call.
+}
+\examples{
+res_multinom(c(1, 3, 6), size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1))
+}
+\seealso{
+Other res_dist:
+\code{\link[=res_bern]{res_bern()}},
+\code{\link[=res_beta_binom]{res_beta_binom()}},
+\code{\link[=res_binom]{res_binom()}},
+\code{\link[=res_gamma]{res_gamma()}},
+\code{\link[=res_gamma_pois]{res_gamma_pois()}},
+\code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}},
+\code{\link[=res_lnorm]{res_lnorm()}},
+\code{\link[=res_neg_binom]{res_neg_binom()}},
+\code{\link[=res_norm]{res_norm()}},
+\code{\link[=res_pois]{res_pois()}},
+\code{\link[=res_pois_zi]{res_pois_zi()}},
+\code{\link[=res_skewlnorm]{res_skewlnorm()}},
+\code{\link[=res_skewnorm]{res_skewnorm()}},
+\code{\link[=res_student]{res_student()}}
+}
+\concept{res_dist}
diff --git a/man/res_neg_binom.Rd b/man/res_neg_binom.Rd
index 50cf1915..f92338ba 100644
--- a/man/res_neg_binom.Rd
+++ b/man/res_neg_binom.Rd
@@ -37,6 +37,7 @@ Other res_dist:
 \code{\link[=res_gamma_pois]{res_gamma_pois()}},
 \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}},
 \code{\link[=res_lnorm]{res_lnorm()}},
+\code{\link[=res_multinom]{res_multinom()}},
 \code{\link[=res_norm]{res_norm()}},
 \code{\link[=res_pois]{res_pois()}},
 \code{\link[=res_pois_zi]{res_pois_zi()}},
diff --git a/man/res_norm.Rd b/man/res_norm.Rd
index ebabf927..3c6eff54 100644
--- a/man/res_norm.Rd
+++ b/man/res_norm.Rd
@@ -36,6 +36,7 @@ Other res_dist:
 \code{\link[=res_gamma_pois]{res_gamma_pois()}},
 \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}},
 \code{\link[=res_lnorm]{res_lnorm()}},
+\code{\link[=res_multinom]{res_multinom()}},
 \code{\link[=res_neg_binom]{res_neg_binom()}},
 \code{\link[=res_pois]{res_pois()}},
 \code{\link[=res_pois_zi]{res_pois_zi()}},
diff --git a/man/res_pois.Rd b/man/res_pois.Rd
index 6da9d015..193bbea0 100644
--- a/man/res_pois.Rd
+++ b/man/res_pois.Rd
@@ -34,6 +34,7 @@ Other res_dist:
 \code{\link[=res_gamma_pois]{res_gamma_pois()}},
 \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}},
 \code{\link[=res_lnorm]{res_lnorm()}},
+\code{\link[=res_multinom]{res_multinom()}},
 \code{\link[=res_neg_binom]{res_neg_binom()}},
 \code{\link[=res_norm]{res_norm()}},
 \code{\link[=res_pois_zi]{res_pois_zi()}},
diff --git a/man/res_pois_zi.Rd b/man/res_pois_zi.Rd
index 8b2aaea2..c8706f20 100644
--- a/man/res_pois_zi.Rd
+++ b/man/res_pois_zi.Rd
@@ -36,6 +36,7 @@ Other res_dist:
 \code{\link[=res_gamma_pois]{res_gamma_pois()}},
 \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}},
 \code{\link[=res_lnorm]{res_lnorm()}},
+\code{\link[=res_multinom]{res_multinom()}},
 \code{\link[=res_neg_binom]{res_neg_binom()}},
 \code{\link[=res_norm]{res_norm()}},
 \code{\link[=res_pois]{res_pois()}},
diff --git a/man/res_skewlnorm.Rd b/man/res_skewlnorm.Rd
index 98364d56..535b3bee 100644
--- a/man/res_skewlnorm.Rd
+++ b/man/res_skewlnorm.Rd
@@ -48,6 +48,7 @@ Other res_dist:
 \code{\link[=res_gamma_pois]{res_gamma_pois()}},
 \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}},
 \code{\link[=res_lnorm]{res_lnorm()}},
+\code{\link[=res_multinom]{res_multinom()}},
 \code{\link[=res_neg_binom]{res_neg_binom()}},
 \code{\link[=res_norm]{res_norm()}},
 \code{\link[=res_pois]{res_pois()}},
diff --git a/man/res_skewnorm.Rd b/man/res_skewnorm.Rd
index f9ec34cf..7af499f1 100644
--- a/man/res_skewnorm.Rd
+++ b/man/res_skewnorm.Rd
@@ -40,6 +40,7 @@ Other res_dist:
 \code{\link[=res_gamma_pois]{res_gamma_pois()}},
 \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}},
 \code{\link[=res_lnorm]{res_lnorm()}},
+\code{\link[=res_multinom]{res_multinom()}},
 \code{\link[=res_neg_binom]{res_neg_binom()}},
 \code{\link[=res_norm]{res_norm()}},
 \code{\link[=res_pois]{res_pois()}},
diff --git a/man/res_student.Rd b/man/res_student.Rd
index 9f9bf7b8..13d5c8cd 100644
--- a/man/res_student.Rd
+++ b/man/res_student.Rd
@@ -39,6 +39,7 @@ Other res_dist:
 \code{\link[=res_gamma_pois]{res_gamma_pois()}},
 \code{\link[=res_gamma_pois_zi]{res_gamma_pois_zi()}},
 \code{\link[=res_lnorm]{res_lnorm()}},
+\code{\link[=res_multinom]{res_multinom()}},
 \code{\link[=res_neg_binom]{res_neg_binom()}},
 \code{\link[=res_norm]{res_norm()}},
 \code{\link[=res_pois]{res_pois()}},
diff --git a/tests/testthat/test-dev.R b/tests/testthat/test-dev.R
index 1afb5c5e..b4fcb0cc 100644
--- a/tests/testthat/test-dev.R
+++ b/tests/testthat/test-dev.R
@@ -473,6 +473,103 @@ test_that("dev_lnorm deviance", {
   expect_equal(deviance, deviance(mod))
 })
 
+test_that("dev_multinom", {
+  expect_identical(
+    dev_multinom(numeric(0), numeric(0), numeric(0), group = numeric(0)),
+    numeric(0)
+  )
+  # a mismatched, non-recyclable length errors clearly instead of being
+  # silently (and wrongly) recycled/NA-padded by rep_len()
+  expect_error(
+    dev_multinom(1:3, 10, numeric(0), group = c(1, 1, 1)),
+    "must be all zero length or the same length"
+  )
+  expect_error(
+    dev_multinom(1:5, c(10, 10, 10), c(0.2, 0.3, 0.5), group = c(1, 1, 1, 2, 2)),
+    "must be all zero length or the same length"
+  )
+  expect_identical(
+    dev_multinom(c(2, 3, 5), size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1)),
+    c(0, 0, 0)
+  )
+  expect_identical(
+    dev_multinom(c(1, 4), c(10, 10), c(0.2, 0.8), group = c(1, 1)),
+    dev_pois(c(1, 4), c(10, 10) * c(0.2, 0.8))
+  )
+  # row-level deviance identity: summing the per-cell deviances within a
+  # trial recovers the classic multinomial saturated-model deviance
+  x <- c(1, 3, 6)
+  size <- 10
+  prob <- c(0.2, 0.3, 0.5)
+  group <- c(1, 1, 1)
+  expect_equal(
+    sum(dev_multinom(x, size, prob, group)),
+    2 * sum(x * log(x / (size * prob)))
+  )
+  # group is validated: prob must sum to 1 per group, and every group must
+  # have at least 2 rows and match the modal row count
+  expect_error(
+    dev_multinom(4, 10, 1, group = 1),
+    "must contain at least 2 rows"
+  )
+  expect_error(
+    dev_multinom(c(1, 2), c(10, 10), c(0.2, 0.2), group = c(1, 1)),
+    "`prob` must sum to 1"
+  )
+})
+
+test_that("dev_multinom matches glmnet deviance", {
+  skip_if_not_installed("glmnet")
+  withr::with_seed(7, {
+    n <- 150
+    K <- 4
+    x <- cbind(rnorm(n), rnorm(n), rnorm(n))
+    eta <- cbind(
+      0,
+      0.9 * x[, 1] - 0.4 * x[, 2],
+      -0.3 * x[, 1] + 0.6 * x[, 3],
+      0.5 * x[, 2] - 0.2 * x[, 3]
+    )
+    prob_true <- exp(eta) / rowSums(exp(eta))
+    y_wide <- t(sapply(seq_len(n), function(i) {
+      stats::rmultinom(1, size = 1, prob = prob_true[i, ])
+    }))
+    colnames(y_wide) <- paste0("cat", seq_len(K))
+
+    group <- rep(seq_len(n), each = K)
+    for (lambda in c(0, 0.01, 0.05)) {
+      fit <- glmnet::glmnet(
+        x,
+        y_wide,
+        family = "multinomial",
+        lambda = lambda,
+        intercept = TRUE
+      )
+      dev_glmnet <- deviance(fit)
+      prob_hat <- predict(fit, newx = x, type = "response")[, , 1]
+      expect_equal(
+        sum(dev_multinom(
+          as.vector(t(y_wide)),
+          size = 1,
+          prob = as.vector(t(prob_hat)),
+          group = group
+        )),
+        as.numeric(dev_glmnet)
+      )
+    }
+  })
+})
+
+test_that("dev_multinom res", {
+  x <- c(1, 3, 6)
+  size <- 10
+  prob <- c(0.2, 0.3, 0.5)
+  group <- c(1, 1, 1)
+  res <- dev_multinom(x, size, prob, group, res = TRUE)
+  expect_equal(sign(res), sign(x - size * prob))
+  expect_equal(sum(res^2), sum(dev_multinom(x, size, prob, group)))
+})
+
 test_that("dev_neg_binom", {
   expect_identical(
     dev_neg_binom(integer(0), integer(0), integer(0)),
diff --git a/tests/testthat/test-internal.R b/tests/testthat/test-internal.R
new file mode 100644
index 00000000..c6fa4034
--- /dev/null
+++ b/tests/testthat/test-internal.R
@@ -0,0 +1,72 @@
+test_that("chk_multinom_group singleton group errors", {
+  expect_error(
+    chk_multinom_group(size = 10, prob = 1, group = 1),
+    "must contain at least 2 rows"
+  )
+})
+
+test_that("chk_multinom_group size check ignores NA but catches known mismatches", {
+  expect_no_error(chk_multinom_group(c(10, 10, NA), c(0.2, 0.3, 0.5), c(1, 1, 1)))
+  expect_error(
+    chk_multinom_group(c(10, 20, NA), c(0.2, 0.3, 0.5), c(1, 1, 1)),
+    "`size` must be the same for every row belonging to the same `group`"
+  )
+})
+
+test_that("chk_multinom_group prob check ignores NA but still catches known values that already break the sum", {
+  # a lone NA prob shouldn't skip validation of the other, already-known values
+  expect_error(
+    chk_multinom_group(c(10, 10, 10), c(0.9, 0.9, NA), c(1, 1, 1)),
+    "`prob` must sum to 1 for every `group`"
+  )
+  # known values that don't yet exceed 1 are fine to leave for the NA to complete
+  expect_no_error(chk_multinom_group(c(10, 10, 10), c(0.4, 0.3, NA), c(1, 1, 1)))
+  # a fully-known group must still sum to exactly 1
+  expect_no_error(chk_multinom_group(c(10, 10, 10), c(0.2, 0.3, 0.5), c(1, 1, 1)))
+  expect_error(
+    chk_multinom_group(c(10, 10, 10), c(0.2, 0.3, 0.4), c(1, 1, 1)),
+    "`prob` must sum to 1 for every `group`"
+  )
+})
+
+test_that("chk_multinom_group modal row count check", {
+  expect_no_error(
+    chk_multinom_group(
+      c(10, 10, 10, 6, 6, 6),
+      c(0.2, 0.3, 0.5, 0.2, 0.3, 0.5),
+      c(1, 1, 1, 2, 2, 2)
+    )
+  )
+  expect_error(
+    chk_multinom_group(
+      c(10, 10, 10, 6, 6),
+      c(0.2, 0.3, 0.5, 0.5, 0.5),
+      c(1, 1, 1, 2, 2)
+    ),
+    "Every `group` should have the same number of rows"
+  )
+})
+
+test_that("chk_multinom_group breaks an exact row-count tie toward the smaller size", {
+  # 2 groups of 2 rows, 2 groups of 3 rows -- documented to prefer the
+  # smaller count (table()'s ascending sort + which.max()'s first-match)
+  expect_error(
+    chk_multinom_group(
+      c(10, 10, 10, 10, 10, 10, 10, 10, 10, 10),
+      c(0.5, 0.5, 0.2, 0.3, 0.5, 0.5, 0.5, 0.2, 0.3, 0.5),
+      c(1, 1, 2, 2, 2, 3, 3, 4, 4, 4)
+    ),
+    "should have the same number of rows \\(2,"
+  )
+})
+
+test_that("multinom_row_na flags a whole group when any size/prob in it is NA", {
+  expect_identical(
+    multinom_row_na(c(10, NA, 10, 10), c(0.2, 0.8, 0.4, 0.6), c(1, 1, 2, 2)),
+    c(TRUE, TRUE, FALSE, FALSE)
+  )
+  expect_identical(
+    multinom_row_na(c(10, 10), c(0.2, 0.8), c(1, 1)),
+    c(FALSE, FALSE)
+  )
+})
diff --git a/tests/testthat/test-log-lik.R b/tests/testthat/test-log-lik.R
index 863441ca..37abfd1d 100644
--- a/tests/testthat/test-log-lik.R
+++ b/tests/testthat/test-log-lik.R
@@ -332,6 +332,102 @@ test_that("log_lik_lnorm truncated", {
   )
 })
 
+test_that("log_lik_multinom", {
+  expect_identical(
+    log_lik_multinom(numeric(0), numeric(0), numeric(0), numeric(0)),
+    numeric(0)
+  )
+  # a mismatched, non-recyclable length errors clearly instead of being
+  # silently (and wrongly) recycled by rep_len(), or silently returning
+  # numeric(0) just because one argument happened to be empty
+  expect_error(
+    log_lik_multinom(1:3, 10, numeric(0), c(1, 1, 1)),
+    "must be all zero length or the same length"
+  )
+  expect_error(
+    log_lik_multinom(1:5, c(10, 10, 10), c(0.2, 0.3, 0.5), c(1, 1, 1, 2, 2)),
+    "must be all zero length or the same length"
+  )
+  expect_error(
+    log_lik_multinom(c(1, 3, 6), c(10, 10, 5), c(0.2, 0.3, 0.5), c(1, 1, 1)),
+    "`size` must be the same for every row belonging to the same `group`"
+  )
+  expect_error(
+    log_lik_multinom(c(1, 3, 6), 10, c(0.2, 0.3, 0.4), c(1, 1, 1)),
+    "`prob` must sum to 1 for every `group`"
+  )
+  expect_error(
+    log_lik_multinom(4, 10, 1, 1),
+    "must contain at least 2 rows"
+  )
+  # group must not contain NA -- there's no way to know which trial an
+  # unlabelled row belongs to
+  expect_error(
+    log_lik_multinom(c(4, 6, 10), c(10, 10, 10), c(0.4, 0.6, 1), c(1, 1, NA)),
+    "must not have any missing values"
+  )
+  expect_error(
+    log_lik_multinom(c(4, 6), c(10, 10), c(0.4, 0.6), c(NA, NA)),
+    "must not have any missing values"
+  )
+  # NA in size or prob is a structural input shared by the whole trial, so
+  # it makes the whole trial's result NA, not just the row where it appears
+  expect_identical(
+    log_lik_multinom(c(4, 6), c(10, NA), c(0.4, 0.6), c(1, 1)),
+    c(NA_real_, NA_real_)
+  )
+  # NA in x is just a missing observation for that one category -- it
+  # doesn't need the rest of the trial to compute its own log_lik, so it
+  # doesn't taint sibling rows the way NA size/prob/group does
+  res <- log_lik_multinom(c(4, NA), c(10, 10), c(0.4, 0.6), c(1, 1))
+  expect_false(is.na(res[1]))
+  expect_true(is.na(res[2]))
+  # an NA elsewhere doesn't leak into an unrelated, fully-known group (same
+  # number of categories in both groups, since groups must match on that)
+  res <- log_lik_multinom(
+    c(4, 3, 3, 2, 2, 2),
+    c(10, 10, NA, 6, 6, 6),
+    c(0.2, 0.3, 0.5, 0.2, 0.3, 0.5),
+    c(1, 1, 1, 2, 2, 2)
+  )
+  expect_identical(res[1:3], c(NA_real_, NA_real_, NA_real_))
+  expect_false(anyNA(res[4:6]))
+  # every group must have the same number of rows (categories) as the most
+  # common number of rows per group in the data
+  expect_error(
+    log_lik_multinom(
+      c(4, 3, 3, 5, 5, 3, 3),
+      c(10, 10, 10, 10, 10, 6, 6),
+      c(0.2, 0.3, 0.5, 0.5, 0.5, 0.5, 0.5),
+      c(1, 1, 1, 2, 2, 3, 3)
+    ),
+    "Every `group` should have the same number of rows"
+  )
+  expect_equal(
+    sum(log_lik_multinom(
+      c(1, 3, 6),
+      size = 10,
+      prob = c(0.2, 0.3, 0.5),
+      group = c(1, 1, 1)
+    )),
+    dmultinom(c(1, 3, 6), size = 10, prob = c(0.2, 0.3, 0.5), log = TRUE)
+  )
+  # multiple trials in long format (same number of categories per group)
+  x <- c(1, 3, 6, 2, 1, 1)
+  size <- c(10, 10, 10, 4, 4, 4)
+  prob <- c(0.2, 0.3, 0.5, 0.5, 0.25, 0.25)
+  group <- c(1, 1, 1, 2, 2, 2)
+  ll <- log_lik_multinom(x, size, prob, group)
+  expect_equal(
+    sum(ll[group == 1]),
+    dmultinom(x[group == 1], size = 10, prob = prob[group == 1], log = TRUE)
+  )
+  expect_equal(
+    sum(ll[group == 2]),
+    dmultinom(x[group == 2], size = 4, prob = prob[group == 2], log = TRUE)
+  )
+})
+
 test_that("log_lik_neg_binom", {
   expect_identical(
     log_lik_neg_binom(0, 2, 1),
diff --git a/tests/testthat/test-ran.R b/tests/testthat/test-ran.R
index 4fcc5e00..68b11318 100644
--- a/tests/testthat/test-ran.R
+++ b/tests/testthat/test-ran.R
@@ -107,6 +107,96 @@ test_that("ran_gamma_pois", {
   })
 })
 
+test_that("ran_multinom", {
+  expect_identical(
+    ran_multinom(size = numeric(0), prob = numeric(0), group = numeric(0)),
+    integer(0)
+  )
+  # a mismatched, non-recyclable length errors clearly instead of being
+  # silently (and wrongly) recycled by rep_len()
+  expect_error(
+    ran_multinom(size = c(10, 10, 10), prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1, 2, 2)),
+    "must be all zero length or the same length"
+  )
+  expect_error(
+    ran_multinom(size = c(10, 5), prob = c(0.5, 0.5), group = c(1, 1)),
+    "`size` must be the same for every row belonging to the same `group`"
+  )
+  expect_error(
+    ran_multinom(size = 10, prob = c(0.5, 0.4), group = c(1, 1)),
+    "`prob` must sum to 1 for every `group`"
+  )
+  expect_error(
+    ran_multinom(size = 10, prob = 1, group = 1),
+    "must contain at least 2 rows"
+  )
+  expect_error(
+    ran_multinom(
+      size = c(10, 10, 10, 10),
+      prob = c(0.2, 0.3, 0.5, 1),
+      group = c(1, 1, 1, 2)
+    ),
+    "must contain at least 2 rows"
+  )
+  # every group must have the same number of rows (categories) as the most
+  # common number of rows per group in the data (ordinary multinomial
+  # logistic regression assumes a fixed set of categories for every trial)
+  expect_error(
+    ran_multinom(
+      size = c(10, 10, 10, 10, 10, 10, 6, 6),
+      prob = c(0.2, 0.3, 0.5, 0.2, 0.3, 0.5, 0.5, 0.5),
+      group = c(1, 1, 1, 2, 2, 2, 3, 3)
+    ),
+    "Every `group` should have the same number of rows"
+  )
+  # group must not contain NA -- there's no way to know which trial an
+  # unlabelled row belongs to
+  expect_error(
+    ran_multinom(size = c(10, 10, 10), prob = c(0.4, 0.6, 1), group = c(1, 1, NA)),
+    "must not have any missing values"
+  )
+  expect_error(
+    ran_multinom(size = c(10, 10), prob = c(0.4, 0.6), group = c(NA, NA)),
+    "must not have any missing values"
+  )
+  # NA in size or prob is a structural input shared by the whole trial (the
+  # categories aren't independent draws), so it makes the whole trial's
+  # result NA, not just the row where the NA appears
+  expect_identical(
+    ran_multinom(size = c(10, NA), prob = c(0.4, 0.6), group = c(1, 1)),
+    c(NA_integer_, NA_integer_)
+  )
+  expect_identical(
+    ran_multinom(size = c(10, 10), prob = c(0.4, NA), group = c(1, 1)),
+    c(NA_integer_, NA_integer_)
+  )
+  # an NA elsewhere doesn't leak into an unrelated, fully-known group (same
+  # number of categories in both groups, since groups must match on that)
+  withr::with_seed(101, {
+    x <- ran_multinom(
+      size = c(10, 10, NA, 6, 6, 6),
+      prob = c(0.2, 0.3, 0.5, 0.2, 0.3, 0.5),
+      group = c(1, 1, 1, 2, 2, 2)
+    )
+    expect_identical(x[1:3], c(NA_integer_, NA_integer_, NA_integer_))
+    expect_identical(sum(x[4:6]), 6L)
+  })
+  withr::with_seed(101, {
+    x <- ran_multinom(size = 10, prob = c(0.2, 0.3, 0.5), group = c(1, 1, 1))
+    expect_identical(sum(x), 10L)
+    expect_length(x, 3L)
+  })
+  withr::with_seed(101, {
+    x <- ran_multinom(
+      size = c(10, 10, 6, 6),
+      prob = c(0.2, 0.8, 0.5, 0.5),
+      group = c(1, 1, 2, 2)
+    )
+    expect_identical(x[1] + x[2], 10L)
+    expect_identical(x[3] + x[4], 6L)
+  })
+})
+
 test_that("ran_neg_binom", {
   expect_error(ran_neg_binom(NA_integer_))
   expect_error(ran_neg_binom(integer(0)))
diff --git a/tests/testthat/test-res.R b/tests/testthat/test-res.R
index e7506af7..0fbf9c5b 100644
--- a/tests/testthat/test-res.R
+++ b/tests/testthat/test-res.R
@@ -507,6 +507,149 @@ test_that("res_gamma_pois simulate", {
   })
 })
 
+test_that("res_multinom", {
+  x <- c(1, 3, 6)
+  size <- 10
+  prob <- c(0.2, 0.3, 0.5)
+  group <- c(1, 1, 1)
+
+  expect_identical(
+    res_multinom(numeric(0), numeric(0), numeric(0), numeric(0)),
+    numeric(0)
+  )
+  # a mismatched, non-recyclable length errors clearly instead of being
+  # silently (and wrongly) recycled/NA-padded by rep_len()
+  expect_error(
+    res_multinom(1:3, 10, numeric(0), c(1, 1, 1)),
+    "must be all zero length or the same length"
+  )
+  expect_error(
+    res_multinom(1:5, c(10, 10, 10), c(0.2, 0.3, 0.5), c(1, 1, 1, 2, 2)),
+    "must be all zero length or the same length"
+  )
+
+  expect_equal(
+    res_multinom(x, size, prob, group, type = "dev"),
+    dev_multinom(x, size, prob, group, res = TRUE)
+  )
+  expect_equal(
+    res_multinom(x, size, prob, group, type = "raw"),
+    x - size * prob
+  )
+  # a category count is marginally binomial, so the standardized residual
+  # matches res_binom exactly
+  expect_equal(
+    res_multinom(x, size, prob, group, type = "standardized"),
+    res_binom(x, size, prob, type = "standardized")
+  )
+  expect_equal(
+    res_multinom(x, size, prob, group, type = "data"),
+    x
+  )
+  expect_error(res_multinom(x, size, prob, group, type = "unknown"))
+  expect_error(
+    res_multinom(4, 10, 1, 1, simulate = TRUE),
+    "must contain at least 2 rows"
+  )
+  # NA in size or prob is a structural input shared by the whole trial, so
+  # it makes the whole trial's result NA, not just the row where it appears
+  expect_identical(
+    res_multinom(c(4, 6), c(10, NA), c(0.4, 0.6), c(1, 1), type = "data", simulate = TRUE),
+    c(NA_integer_, NA_integer_)
+  )
+  # every group must have the same number of rows (categories) as the most
+  # common number of rows per group in the data
+  expect_error(
+    res_multinom(
+      c(4, 3, 3, 5, 5),
+      c(10, 10, 10, 6, 6),
+      c(0.2, 0.3, 0.5, 0.5, 0.5),
+      c(1, 1, 1, 2, 2),
+      type = "data",
+      simulate = TRUE
+    ),
+    "Every `group` should have the same number of rows"
+  )
+  # group must not contain NA -- there's no way to know which trial an
+  # unlabelled row belongs to
+  expect_error(
+    res_multinom(
+      c(4, 6, 10),
+      c(10, 10, 10),
+      c(0.4, 0.6, 1),
+      c(1, 1, NA),
+      type = "data",
+      simulate = TRUE
+    ),
+    "must not have any missing values"
+  )
+  expect_error(
+    res_multinom(
+      c(4, 6),
+      c(10, 10),
+      c(0.4, 0.6),
+      c(NA, NA),
+      type = "data",
+      simulate = TRUE
+    ),
+    "must not have any missing values"
+  )
+  # group is validated even when simulate = FALSE, so a bad prob sum can't
+  # silently fall through to a wrong (rather than an errored) residual
+  expect_error(
+    res_multinom(4, 10, 1, 1, simulate = FALSE),
+    "must contain at least 2 rows"
+  )
+  expect_error(
+    res_multinom(c(4, 6), c(10, 10), c(0.4, 0.4), c(1, 1), simulate = FALSE),
+    "`prob` must sum to 1"
+  )
+
+  # sum of squared deviance residuals recovers the row-level deviance
+  expect_equal(
+    sum(res_multinom(x, size, prob, group, type = "dev")^2),
+    2 * sum(x * log(x / (size * prob)))
+  )
+
+  withr::with_seed(101, {
+    sim <- res_multinom(x, size, prob, group, type = "data", simulate = TRUE)
+    expect_identical(sum(sim), as.integer(size))
+    expect_length(sim, 3L)
+  })
+})
+
+test_that("res_multinom simulate", {
+  n_group <- 10000
+  size <- rep(10, n_group * 3)
+  prob <- rep(c(0.2, 0.3, 0.5), n_group)
+  group <- rep(seq_len(n_group), each = 3)
+
+  withr::with_seed(101, {
+    res <- res_multinom(
+      rep(0, n_group * 3),
+      size = size,
+      prob = prob,
+      group = group,
+      simulate = TRUE,
+      type = "dev"
+    )
+    expect_equal(mean(res), -0.0794000903313601)
+    expect_equal(sd(res), 0.859943455555756)
+  })
+  withr::with_seed(101, {
+    res <- res_multinom(
+      rep(0, n_group * 3),
+      size = size,
+      prob = prob,
+      group = group,
+      simulate = TRUE,
+      type = "standardized"
+    )
+    expect_equal(mean(res), 0.000255549022088425)
+    expect_equal(sd(res), 1.00164264126585)
+  })
+})
+
 test_that("res_neg_binom", {
   expect_identical(
     res_neg_binom(integer(0), integer(0), integer(0)),