diff --git a/src/functions-reference/unbounded_discrete_distributions.qmd b/src/functions-reference/unbounded_discrete_distributions.qmd
index 9130bda00..cb49ebd83 100644
--- a/src/functions-reference/unbounded_discrete_distributions.qmd
+++ b/src/functions-reference/unbounded_discrete_distributions.qmd
@@ -586,3 +586,36 @@ The log Poisson probability mass of `y` given the log-rate `alpha + x * beta`.
The log Poisson probability mass of `y` given the log-rate `alpha + x * beta`
dropping constant additive terms.
{{< since 2.25 >}}
+
+## Beta negative binomial distribution {#beta-neg-binomial}
+
+### Probability mass function
+
+If $r \in \mathbb{R}^+$, $\alpha \in \mathbb{R}^+$, and $\beta \in \mathbb{R}^+$, then for $n \in \mathbb{N}$, \begin{equation*}
+\text{BetaNegBinomial}(n|r,\alpha,\beta) = \frac {\Gamma (n+r )}{n!\;\Gamma (r )}
+\frac {\mathrm {B} (\beta+n,\alpha +r )}{\mathrm {B} (\beta,\alpha )}. \end{equation*}
+
+### Distribution statement
+
+`n ~ ` **`beta_neg_binomial`**`(r,alpha,beta)`
+
+Increment target log probability density with `beta_neg_binomial_lupmf(n | r, alpha, beta)`.
+{{< since 2.36 >}}
+
+\index{{\tt \bfseries beta\_neg\_binomial }!sampling statement|hyperpage}
+
+### Stan functions
+
+
+\index{{\tt \bfseries beta\_neg\_binomial\_lpmf }!{\tt (ints n \textbar\ reals r, reals alpha, reals beta): real}|hyperpage}
+
+`real` **`beta_neg_binomial_lpmf`**`(ints n | reals r, reals alpha, reals beta)`
\newline
+The log beta negative binomial probability mass of `n` given parameters `r`, `alpha` and `beta`.
+{{< since 2.36 >}}
+
+
+\index{{\tt \bfseries beta\_neg\_binomial\_lupmf }!{\tt (ints n \textbar\ reals r, reals alpha, reals beta): real}|hyperpage}
+
+`real` **`beta_neg_binomial_lupmf`**`(ints n | reals r, reals alpha, reals beta)`
\newline
+The log beta negative binomial probability mass of `n` given parameters `r`, `alpha` and `beta` dropping constant additive terms.
+{{< since 2.36 >}}
\ No newline at end of file