diff --git a/src/functions-reference/bounded_discrete_distributions.Rmd b/src/functions-reference/bounded_discrete_distributions.Rmd
index bdaea5c80..152694ed0 100644
--- a/src/functions-reference/bounded_discrete_distributions.Rmd
+++ b/src/functions-reference/bounded_discrete_distributions.Rmd
@@ -155,6 +155,133 @@ The log binomial probability mass of n successes in N trials given
 logit-scaled chance of success alpha dropping constant additive terms
 `r since("2.25")`
 
+## Binomial-logit generalized linear model (Logistic Regression) {#binomial-logit-glm}
+
+Stan also supplies a single function for a generalized linear model
+with binomial likelihood and logit link function, i.e., a function
+for logistic regression with aggregated outcomes. This provides a more efficient
+implementation of logistic regression than a manually written
+regression in terms of a Binomial likelihood and matrix
+multiplication.
+
+### Probability mass function
+
+Suppose $N \in \mathbb{N}$, $x\in \mathbb{R}^{n\cdot m}, \alpha \in \mathbb{R}^n, \beta \in \mathbb{R}^m$, and $n \in
+\{0,\ldots,N\}$.  Then
+\begin{align*}
+  &\text{BinomialLogitGLM}(n~|~N, x, \alpha, \beta) = \text{Binomial}(n~|~N,\text{logit}^{-1}(\alpha_i + x_i \cdot \beta)) \\
+  &= \binom{N}{n} \left( \text{logit}^{-1}(\alpha_i + \sum_{1\leq j\leq m}x_{ij}\cdot \beta_j) \right)^{n}  \left( 1 - \text{logit}^{-1}(\alpha_i + \sum_{1\leq j\leq m}x_{ij}\cdot \beta_j) \right)^{N - n}.
+\end{align*}
+
+### Sampling statement
+
+`n ~ ` **`binomial_logit_glm`**`(N, x, alpha, beta)`
+
+Increment target log probability density with `binomial_logit_glm_lupmf(n | N, x, alpha, beta)`.
+`r since("2.34")`
+
+<!-- real; binomial_logit_glm ~; -->
+\index{{\tt \bfseries binomial\_logit\_glm }!sampling statement|hyperpage}
+
+### Stan Functions
+
+<!-- real; binomial_logit_glm_lpmf; (int n | int N, matrix x, real alpha, vector beta); -->
+\index{{\tt \bfseries binomial\_logit\_glm\_lpmf  }!{\tt (int n \textbar\ int N, matrix x, real alpha, vector beta): real}|hyperpage}
+
+`real` **`binomial_logit_glm_lpmf`**`(int n | int N, matrix x, real alpha, vector beta)`<br>\newline
+The log binomial probability mass of n given N trials and chance of success
+`inv_logit(alpha + x * beta)`.
+`r since("2.34")`
+
+<!-- real; binomial_logit_glm_lupmf; (int n | int N, matrix x, real alpha, vector beta); -->
+\index{{\tt \bfseries binomial\_logit\_glm\_lupmf  }!{\tt (int n \textbar\ int N, matrix x, real alpha, vector beta): real}|hyperpage}
+
+`real` **`binomial_logit_glm_lupmf`**`(int n | int N, matrix x, real alpha, vector beta)`<br>\newline
+The log binomial probability mass of n given N trials and chance of success
+`inv_logit(alpha + x * beta)` dropping constant additive terms.
+`r since("2.34")`
+
+<!-- real; binomial_logit_glm_lpmf; (int n | int N, matrix x, vector alpha, vector beta); -->
+\index{{\tt \bfseries binomial\_logit\_glm\_lpmf  }!{\tt (int n \textbar\ int N, matrix x, vector alpha, vector beta): real}|hyperpage}
+
+`real` **`binomial_logit_glm_lpmf`**`(int n | int N, matrix x, vector alpha, vector beta)`<br>\newline
+The log binomial probability mass of n given N trials and chance of success
+`inv_logit(alpha + x * beta)`.
+`r since("2.34")`
+
+<!-- real; binomial_logit_glm_lupmf; (int n | int N, matrix x, vector alpha, vector beta); -->
+\index{{\tt \bfseries binomial\_logit\_glm\_lupmf  }!{\tt (int n \textbar\ int N, matrix x, vector alpha, vector beta): real}|hyperpage}
+
+`real` **`binomial_logit_glm_lupmf`**`(int n | int N, matrix x, vector alpha, vector beta)`<br>\newline
+The log binomial probability mass of n given N trials and chance of success
+`inv_logit(alpha + x * beta)` dropping constant additive terms.
+`r since("2.34")`
+
+<!-- real; binomial_logit_glm_lpmf; (array[] int n | array[] int N, row_vector x, real alpha, vector beta); -->
+\index{{\tt \bfseries binomial\_logit\_glm\_lpmf  }!{\tt (array[] int n \textbar\ array[] int N, row\_vector x, real alpha, vector beta): real}|hyperpage}
+
+`real` **`binomial_logit_glm_lpmf`**`(array[] int n | array[] int N, row_vector x, real alpha, vector beta)`<br>\newline
+The log binomial probability mass of n given N trials and chance of success
+`inv_logit(alpha + x * beta)`.
+`r since("2.34")`
+
+<!-- real; binomial_logit_glm_lupmf; (array[] int n | array[] int N, row_vector x, real alpha, vector beta); -->
+\index{{\tt \bfseries binomial\_logit\_glm\_lupmf  }!{\tt (array[] int n \textbar\ array[] int N, row\_vector x, real alpha, vector beta): real}|hyperpage}
+
+`real` **`binomial_logit_glm_lupmf`**`(array[] int n | array[] int N, row_vector x, real alpha, vector beta)`<br>\newline
+The log binomial probability mass of n given N trials and chance of success
+`inv_logit(alpha + x * beta)` dropping constant additive terms.
+`r since("2.34")`
+
+<!-- real; binomial_logit_glm_lpmf; (array[] int n | array[] int N, row_vector x, vector alpha, vector beta); -->
+\index{{\tt \bfseries binomial\_logit\_glm\_lpmf  }!{\tt (array[] int n \textbar\ array[] int N, row\_vector x, vector alpha, vector beta): real}|hyperpage}
+
+`real` **`binomial_logit_glm_lpmf`**`(array[] int n | array[] int N, row_vector x, vector alpha, vector beta)`<br>\newline
+The log binomial probability mass of n given N trials and chance of success
+`inv_logit(alpha + x * beta)`.
+`r since("2.34")`
+
+<!-- real; binomial_logit_glm_lupmf; (array[] int n | array[] int N, row_vector x, vector alpha, vector beta); -->
+\index{{\tt \bfseries binomial\_logit\_glm\_lupmf  }!{\tt (array[] int n \textbar\ array[] int N, row\_vector x, vector alpha, vector beta): real}|hyperpage}
+
+`real` **`binomial_logit_glm_lupmf`**`(array[] int n | array[] int N, row_vector x, vector alpha, vector beta)`<br>\newline
+The log binomial probability mass of n given N trials and chance of success
+`inv_logit(alpha + x * beta)` dropping constant additive terms.
+`r since("2.34")`
+
+
+<!-- real; binomial_logit_glm_lpmf; (array[] int n | array[] int N, matrix x, real alpha, vector beta); -->
+\index{{\tt \bfseries binomial\_logit\_glm\_lpmf  }!{\tt (array[] int n \textbar\ array[] int N, matrix x, real alpha, vector beta): real}|hyperpage}
+
+`real` **`binomial_logit_glm_lpmf`**`(array[] int n | array[] int N, matrix x, real alpha, vector beta)`<br>\newline
+The log binomial probability mass of n given N trials and chance of success
+`inv_logit(alpha + x * beta)`.
+`r since("2.34")`
+
+<!-- real; binomial_logit_glm_lupmf; (array[] int n | array[] int N, matrix x, real alpha, vector beta); -->
+\index{{\tt \bfseries binomial\_logit\_glm\_lupmf  }!{\tt (array[] int n \textbar\ array[] int N, matrix x, real alpha, vector beta): real}|hyperpage}
+
+`real` **`binomial_logit_glm_lupmf`**`(array[] int n | array[] int N, matrix x, real alpha, vector beta)`<br>\newline
+The log binomial probability mass of n given N trials and chance of success
+`inv_logit(alpha + x * beta)` dropping constant additive terms.
+`r since("2.34")`
+
+<!-- real; binomial_logit_glm_lpmf; (array[] int n | array[] int N, matrix x, vector alpha, vector beta); -->
+\index{{\tt \bfseries binomial\_logit\_glm\_lpmf  }!{\tt (array[] int n \textbar\ array[] int N, matrix x, vector alpha, vector beta): real}|hyperpage}
+
+`real` **`binomial_logit_glm_lpmf`**`(array[] int n | array[] int N, matrix x, vector alpha, vector beta)`<br>\newline
+The log binomial probability mass of n given N trials and chance of success
+`inv_logit(alpha + x * beta)`.
+`r since("2.34")`
+
+<!-- real; binomial_logit_glm_lupmf; (array[] int n | array[] int N, matrix x, vector alpha, vector beta); -->
+\index{{\tt \bfseries binomial\_logit\_glm\_lupmf  }!{\tt (array[] int n \textbar\ array[] int N, matrix x, vector alpha, vector beta): real}|hyperpage}
+
+`real` **`binomial_logit_glm_lupmf`**`(array[] int n | array[] int N, matrix x, vector alpha, vector beta)`<br>\newline
+The log binomial probability mass of n given N trials and chance of success
+`inv_logit(alpha + x * beta)` dropping constant additive terms.
+`r since("2.34")`
+
 ## Beta-binomial distribution
 
 ### Probability mass function