From 1f33f4a93d29c6547dfb8ed1a08c5edf120c5329 Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Mon, 19 Oct 2026 08:56:01 +0300
Subject: [PATCH 1/3] closed-form derivatives for the built-in embedded Laplace
 likelihoods

---
 stan/math/mix/functor/laplace_likelihood.hpp  | 122 +++++++-
 .../laplace_marginal_bernoulli_logit_lpmf.hpp |  66 +++++
 ...place_marginal_neg_binomial_2_log_lpmf.hpp | 154 ++++++++++
 .../laplace_marginal_poisson_log_lpmf.hpp     |  63 ++++
 .../laplace_analytic_derivatives_test.cpp     | 268 ++++++++++++++++++
 5 files changed, 669 insertions(+), 4 deletions(-)
 create mode 100644 test/unit/math/laplace/laplace_analytic_derivatives_test.cpp

diff --git a/stan/math/mix/functor/laplace_likelihood.hpp b/stan/math/mix/functor/laplace_likelihood.hpp
index ec5be6b36d1..75b7d23f468 100644
--- a/stan/math/mix/functor/laplace_likelihood.hpp
+++ b/stan/math/mix/functor/laplace_likelihood.hpp
@@ -10,6 +10,7 @@
 #include <stan/math/mix/functor/hessian_times_vector.hpp>
 
 #include <tuple>
+#include <type_traits>
 #include <utility>
 
 namespace stan {
@@ -21,7 +22,46 @@ namespace math {
  */
 namespace laplace_likelihood {
 
+/**
+ * Specialize this trait with `value = true` for a likelihood functor `F` if
+ * the Hessian of its log likelihood with respect to `theta` is diagonal and
+ * its derivatives with respect to `theta` are known in closed form.  The
+ * specialization must have three static member functions, which take
+ * `theta` and the arguments of `F::operator()` without the message stream,
+ * all with arithmetic scalar types, and return an `Eigen::VectorXd`:
+ *  - `theta_grad(theta, args...)`: the gradient;
+ *  - `neg_hessian_diag(theta, args...)`: the diagonal of the negative
+ *    Hessian;
+ *  - `third_diff(theta, args...)`: the diagonal of the third derivative.
+ * When all arguments are arithmetic, `theta_grad()`, `diagonal_hessian()`,
+ * `block_hessian()` and `third_diff()` below then use these functions
+ * instead of nested autodiff.
+ * @tparam F Type of the likelihood functor
+ */
+template <typename F>
+struct analytic_derivatives : std::false_type {};
+
 namespace internal {
+
+/**
+ * True if the closed-form derivatives of `F` (see `analytic_derivatives`)
+ * are used for these argument types: `F` has them, and `theta` and all
+ * arguments have arithmetic scalar types.
+ */
+template <typename F, typename Theta, typename... Args>
+inline constexpr bool use_analytic_derivatives_v
+    = analytic_derivatives<std::decay_t<F>>::value&&
+        std::is_arithmetic<scalar_type_t<std::decay_t<Theta>>>::value&&
+            std::conjunction_v<
+                std::is_arithmetic<scalar_type_t<std::decay_t<Args>>>...>;
+
+template <typename F, typename Theta, typename... Args>
+using require_analytic_derivatives_t = require_t<
+    std::bool_constant<use_analytic_derivatives_v<F, Theta, Args...>>>;
+
+template <typename F, typename Theta, typename... Args>
+using require_not_analytic_derivatives_t = require_t<
+    std::bool_constant<!use_analytic_derivatives_v<F, Theta, Args...>>>;
 /**
  * @tparam F A functor with `opertor()(Args&&...)` returning a scalar
  * @tparam Theta A class assignable to an Eigen vector type
@@ -53,7 +93,8 @@ inline auto log_likelihood(F&& f, Theta&& theta, Stream* msgs, Args&&... args) {
  * @param args Variadic arguments for the likelihood function.
  */
 template <typename F, typename Theta, typename Stream, typename... Args,
-          require_eigen_vector_vt<std::is_arithmetic, Theta>* = nullptr>
+          require_eigen_vector_vt<std::is_arithmetic, Theta>* = nullptr,
+          require_not_analytic_derivatives_t<F, Theta, Args...>* = nullptr>
 inline auto theta_grad(F&& f, Theta&& theta, Stream* msgs, Args&&... args) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
@@ -104,7 +145,8 @@ inline void ll_arg_grad(F&& f, Theta&& theta, Stream* msgs, Args&&... args) {
  * @param args Variadic arguments for the likelihood function.
  */
 template <typename F, typename Theta, typename Stream, typename... Args,
-          require_eigen_vector_vt<std::is_arithmetic, Theta>* = nullptr>
+          require_eigen_vector_vt<std::is_arithmetic, Theta>* = nullptr,
+          require_not_analytic_derivatives_t<F, Theta, Args...>* = nullptr>
 inline auto diagonal_hessian(F&& f, Theta&& theta, Stream* msgs,
                              Args&&... args) {
   using Eigen::Dynamic;
@@ -134,7 +176,8 @@ inline auto diagonal_hessian(F&& f, Theta&& theta, Stream* msgs,
  * @param args Variadic arguments for the likelihood function.
  */
 template <typename F, typename Theta, typename Stream, typename... Args,
-          require_eigen_vector_vt<std::is_arithmetic, Theta>* = nullptr>
+          require_eigen_vector_vt<std::is_arithmetic, Theta>* = nullptr,
+          require_not_analytic_derivatives_t<F, Theta, Args...>* = nullptr>
 inline auto block_hessian(F&& f, Theta&& theta,
                           const Eigen::Index hessian_block_size, Stream* msgs,
                           Args&&... args) {
@@ -224,7 +267,8 @@ inline auto diff(F&& f, Theta&& theta, const Eigen::Index hessian_block_size,
  * @param args Variadic arguments for likelihood function.
  */
 template <typename F, typename Theta, typename Stream, typename... Args,
-          require_eigen_vector_t<Theta>* = nullptr>
+          require_eigen_vector_t<Theta>* = nullptr,
+          require_not_analytic_derivatives_t<F, Theta, Args...>* = nullptr>
 inline Eigen::VectorXd third_diff(F&& f, Theta&& theta, Stream&& msgs,
                                   Args&&... args) {
   nested_rev_autodiff nested;
@@ -355,6 +399,76 @@ inline auto diff_eta_implicit(F&& f, V_t&& v, Theta&& theta, Stream* msgs,
   grad(f_sum.d_.vi_);
 }
 
+/**
+ * Gradient of `f` with respect to `theta`, from the closed-form derivatives
+ * of `F` (see `analytic_derivatives`).
+ */
+template <typename F, typename Theta, typename Stream, typename... Args,
+          require_eigen_vector_vt<std::is_arithmetic, Theta>* = nullptr,
+          require_analytic_derivatives_t<F, Theta, Args...>* = nullptr>
+inline Eigen::VectorXd theta_grad(F&& f, Theta&& theta, Stream* msgs,
+                                  Args&&... args) {
+  return analytic_derivatives<std::decay_t<F>>::theta_grad(theta, args...);
+}
+
+/**
+ * Diagonal of the negative Hessian of `f` with respect to `theta`, from the
+ * closed-form derivatives of `F` (see `analytic_derivatives`).
+ */
+template <typename F, typename Theta, typename Stream, typename... Args,
+          require_eigen_vector_vt<std::is_arithmetic, Theta>* = nullptr,
+          require_analytic_derivatives_t<F, Theta, Args...>* = nullptr>
+inline Eigen::VectorXd diagonal_hessian(F&& f, Theta&& theta, Stream* msgs,
+                                        Args&&... args) {
+  return analytic_derivatives<std::decay_t<F>>::neg_hessian_diag(theta,
+                                                                 args...);
+}
+
+/**
+ * Negative Hessian of `f` with respect to `theta` as a sparse block-diagonal
+ * matrix, from the closed-form derivatives of `F` (see
+ * `analytic_derivatives`).  The Hessian of such an `F` is diagonal; the
+ * matrix has the sparsity pattern of the autodiff version
+ * (`hessian_block_diag()`): full blocks of size `hessian_block_size`, with
+ * zero off-diagonal entries.
+ */
+template <typename F, typename Theta, typename Stream, typename... Args,
+          require_eigen_vector_vt<std::is_arithmetic, Theta>* = nullptr,
+          require_analytic_derivatives_t<F, Theta, Args...>* = nullptr>
+inline Eigen::SparseMatrix<double> block_hessian(
+    F&& f, Theta&& theta, const Eigen::Index hessian_block_size, Stream* msgs,
+    Args&&... args) {
+  const Eigen::VectorXd neg_hessian_diag
+      = analytic_derivatives<std::decay_t<F>>::neg_hessian_diag(theta, args...);
+  const Eigen::Index theta_size = theta.size();
+  Eigen::SparseMatrix<double> neg_hessian(theta_size, theta_size);
+  neg_hessian.reserve(
+      Eigen::VectorXi::Constant(theta_size, hessian_block_size));
+  const Eigen::Index n_blocks = theta_size / hessian_block_size;
+  for (Eigen::Index j = 0; j < n_blocks; ++j) {
+    const Eigen::Index start = j * hessian_block_size;
+    for (Eigen::Index c = start; c < start + hessian_block_size; ++c) {
+      for (Eigen::Index r = start; r < start + hessian_block_size; ++r) {
+        neg_hessian.insert(r, c) = r == c ? neg_hessian_diag(r) : 0.0;
+      }
+    }
+  }
+  neg_hessian.makeCompressed();
+  return neg_hessian;
+}
+
+/**
+ * Diagonal of the third derivative of `f` with respect to `theta`, from the
+ * closed-form derivatives of `F` (see `analytic_derivatives`).
+ */
+template <typename F, typename Theta, typename Stream, typename... Args,
+          require_eigen_vector_t<Theta>* = nullptr,
+          require_analytic_derivatives_t<F, Theta, Args...>* = nullptr>
+inline Eigen::VectorXd third_diff(F&& f, Theta&& theta, Stream&& msgs,
+                                  Args&&... args) {
+  return analytic_derivatives<std::decay_t<F>>::third_diff(theta, args...);
+}
+
 }  // namespace internal
 
 /**
diff --git a/stan/math/mix/prob/laplace_marginal_bernoulli_logit_lpmf.hpp b/stan/math/mix/prob/laplace_marginal_bernoulli_logit_lpmf.hpp
index 79e1598486d..78e809988c0 100644
--- a/stan/math/mix/prob/laplace_marginal_bernoulli_logit_lpmf.hpp
+++ b/stan/math/mix/prob/laplace_marginal_bernoulli_logit_lpmf.hpp
@@ -3,6 +3,8 @@
 
 #include <stan/math/mix/meta.hpp>
 #include <stan/math/mix/functor/laplace_likelihood.hpp>
+#include <stan/math/prim/fun/inv_logit.hpp>
+#include <type_traits>
 #include <stan/math/mix/functor/laplace_marginal_density.hpp>
 #include <stan/math/rev/fun/elt_multiply.hpp>
 #include <stan/math/rev/fun/exp.hpp>
@@ -58,6 +60,70 @@ struct bernoulli_logit_likelihood {
   }
 };
 
+namespace laplace_likelihood {
+/**
+ * Closed-form derivatives of `bernoulli_logit_likelihood` with respect to
+ * `theta` (see `analytic_derivatives`).  With `p_g = inv_logit(theta_g +
+ * mean_g)` and `q_g = 1 - p_g = inv_logit(-theta_g - mean_g)` (computed
+ * without cancellation), the gradient is `counts_g - n_g p_g`, the negative
+ * Hessian is `n_g p_g q_g`, and the third derivative is
+ * `-n_g p_g q_g (q_g - p_g)`.
+ */
+template <>
+struct analytic_derivatives<bernoulli_logit_likelihood> : std::true_type {
+  // sum of the counts and number of observations of each group
+  template <typename YVec>
+  static void group_sums(Eigen::Index theta_size, const YVec& y,
+                         const std::vector<int>& y_index,
+                         Eigen::VectorXd& counts, Eigen::VectorXd& n) {
+    counts = Eigen::VectorXd::Zero(theta_size);
+    n = Eigen::VectorXd::Zero(theta_size);
+    for (size_t i = 0; i < y_index.size(); i++) {
+      counts(y_index[i] - 1) += y[i];
+      n(y_index[i] - 1) += 1;
+    }
+  }
+
+  template <typename Theta, typename YVec, typename Mean>
+  static void probabilities(const Theta& theta, const YVec& y,
+                            const std::vector<int>& y_index, const Mean& mean,
+                            Eigen::VectorXd& counts, Eigen::VectorXd& n,
+                            Eigen::VectorXd& p, Eigen::VectorXd& q) {
+    group_sums(theta.size(), y, y_index, counts, n);
+    Eigen::VectorXd theta_offset = add(theta, mean);
+    p = inv_logit(theta_offset);
+    q = inv_logit(-theta_offset);
+  }
+
+  template <typename Theta, typename YVec, typename Mean>
+  static Eigen::VectorXd theta_grad(const Theta& theta, const YVec& y,
+                                    const std::vector<int>& y_index,
+                                    const Mean& mean) {
+    Eigen::VectorXd counts, n, p, q;
+    probabilities(theta, y, y_index, mean, counts, n, p, q);
+    return counts - n.cwiseProduct(p);
+  }
+
+  template <typename Theta, typename YVec, typename Mean>
+  static Eigen::VectorXd neg_hessian_diag(const Theta& theta, const YVec& y,
+                                          const std::vector<int>& y_index,
+                                          const Mean& mean) {
+    Eigen::VectorXd counts, n, p, q;
+    probabilities(theta, y, y_index, mean, counts, n, p, q);
+    return n.cwiseProduct(p).cwiseProduct(q);
+  }
+
+  template <typename Theta, typename YVec, typename Mean>
+  static Eigen::VectorXd third_diff(const Theta& theta, const YVec& y,
+                                    const std::vector<int>& y_index,
+                                    const Mean& mean) {
+    Eigen::VectorXd counts, n, p, q;
+    probabilities(theta, y, y_index, mean, counts, n, p, q);
+    return -n.cwiseProduct(p).cwiseProduct(q).cwiseProduct(q - p);
+  }
+};
+}  // namespace laplace_likelihood
+
 /**
  * Wrapper function around the laplace_marginal function for
  * a logistic Bernoulli likelihood. Returns the marginal density
diff --git a/stan/math/mix/prob/laplace_marginal_neg_binomial_2_log_lpmf.hpp b/stan/math/mix/prob/laplace_marginal_neg_binomial_2_log_lpmf.hpp
index 699a54b914a..bc56d3444e3 100644
--- a/stan/math/mix/prob/laplace_marginal_neg_binomial_2_log_lpmf.hpp
+++ b/stan/math/mix/prob/laplace_marginal_neg_binomial_2_log_lpmf.hpp
@@ -2,6 +2,9 @@
 #define STAN_MATH_MIX_PROB_LAPLACE_MARGINAL_NEG_BINOMIAL_2_LOG_LPMF_HPP
 
 #include <stan/math/mix/functor/laplace_likelihood.hpp>
+#include <stan/math/prim/fun/inv_logit.hpp>
+#include <cmath>
+#include <type_traits>
 #include <stan/math/mix/functor/laplace_marginal_density.hpp>
 
 #include <stan/math/rev/core/operator_addition.hpp>
@@ -193,6 +196,157 @@ struct neg_binomial_2_log_likelihood_summary {
   }
 };
 
+namespace laplace_likelihood {
+/**
+ * Closed-form derivatives of `neg_binomial_2_log_likelihood` with respect
+ * to `theta` (see `analytic_derivatives`).  With the dispersion `eta_g` of
+ * group `g` (a scalar `eta` or one entry per group),
+ * `s_g = counts_g + eta_g n_g`, `p_g = inv_logit(theta_g + mean_g -
+ * log(eta_g))` and `q_g = 1 - p_g` (computed without cancellation), the
+ * gradient is `counts_g - s_g p_g`, the negative Hessian is `s_g p_g q_g`,
+ * and the third derivative is `-s_g p_g q_g (q_g - p_g)`.
+ */
+template <>
+struct analytic_derivatives<neg_binomial_2_log_likelihood> : std::true_type {
+  // gradient, negative Hessian and third derivative from the sums of the
+  // counts `counts`, `s = counts + eta n`, and x = theta + mean - log(eta)
+  static Eigen::VectorXd grad(const Eigen::VectorXd& counts,
+                              const Eigen::VectorXd& s,
+                              const Eigen::VectorXd& x) {
+    return counts - s.cwiseProduct(inv_logit(x));
+  }
+  static Eigen::VectorXd neg_hess(const Eigen::VectorXd& s,
+                                  const Eigen::VectorXd& x) {
+    return s.cwiseProduct(inv_logit(x)).cwiseProduct(inv_logit(-x));
+  }
+  static Eigen::VectorXd third(const Eigen::VectorXd& s,
+                               const Eigen::VectorXd& x) {
+    const Eigen::VectorXd p = inv_logit(x);
+    const Eigen::VectorXd q = inv_logit(-x);
+    return -s.cwiseProduct(p).cwiseProduct(q).cwiseProduct(q - p);
+  }
+
+  template <typename Theta, typename Eta, typename Mean>
+  static void parts(const Theta& theta, const Eta& eta,
+                    const std::vector<int>& y, const std::vector<int>& y_index,
+                    const Mean& mean, Eigen::VectorXd& counts,
+                    Eigen::VectorXd& s, Eigen::VectorXd& x) {
+    const Eigen::Index theta_size = theta.size();
+    counts = Eigen::VectorXd::Zero(theta_size);
+    Eigen::VectorXd n = Eigen::VectorXd::Zero(theta_size);
+    for (size_t i = 0; i < y.size(); i++) {
+      counts(y_index[i] - 1) += y[i];
+      n(y_index[i] - 1) += 1;
+    }
+    Eigen::VectorXd eta_g(theta_size);
+    if constexpr (is_stan_scalar<Eta>::value) {
+      eta_g.setConstant(eta);
+    } else {
+      check_size_match("neg_binomial_2_log_likelihood", "eta", eta.size(),
+                       "theta", theta_size);
+      for (Eigen::Index g = 0; g < theta_size; ++g) {
+        eta_g(g) = eta[g];
+      }
+    }
+    s = counts + eta_g.cwiseProduct(n);
+    x = add(theta, mean) - eta_g.array().log().matrix();
+  }
+
+  template <typename Theta, typename Eta, typename Mean>
+  static Eigen::VectorXd theta_grad(const Theta& theta, const Eta& eta,
+                                    const std::vector<int>& y,
+                                    const std::vector<int>& y_index,
+                                    const Mean& mean) {
+    Eigen::VectorXd counts, s, x;
+    parts(theta, eta, y, y_index, mean, counts, s, x);
+    return grad(counts, s, x);
+  }
+
+  template <typename Theta, typename Eta, typename Mean>
+  static Eigen::VectorXd neg_hessian_diag(const Theta& theta, const Eta& eta,
+                                          const std::vector<int>& y,
+                                          const std::vector<int>& y_index,
+                                          const Mean& mean) {
+    Eigen::VectorXd counts, s, x;
+    parts(theta, eta, y, y_index, mean, counts, s, x);
+    return neg_hess(s, x);
+  }
+
+  template <typename Theta, typename Eta, typename Mean>
+  static Eigen::VectorXd third_diff(const Theta& theta, const Eta& eta,
+                                    const std::vector<int>& y,
+                                    const std::vector<int>& y_index,
+                                    const Mean& mean) {
+    Eigen::VectorXd counts, s, x;
+    parts(theta, eta, y, y_index, mean, counts, s, x);
+    return third(s, x);
+  }
+};
+
+/**
+ * Closed-form derivatives of `neg_binomial_2_log_likelihood_summary` with
+ * respect to `theta`; the same formulas as for
+ * `neg_binomial_2_log_likelihood`, with a scalar `eta`.
+ */
+template <>
+struct analytic_derivatives<neg_binomial_2_log_likelihood_summary>
+    : std::true_type {
+  template <typename Theta, typename Eta, typename Mean>
+  static void parts(const Theta& theta, const Eta& eta,
+                    const std::vector<int>& n_per_group,
+                    const std::vector<int>& counts_per_group, const Mean& mean,
+                    Eigen::VectorXd& counts, Eigen::VectorXd& s,
+                    Eigen::VectorXd& x) {
+    const Eigen::Index theta_size = theta.size();
+    counts = Eigen::Map<const Eigen::VectorXi>(counts_per_group.data(),
+                                               counts_per_group.size())
+                 .cast<double>();
+    const Eigen::VectorXd n = Eigen::Map<const Eigen::VectorXi>(
+                                  n_per_group.data(), n_per_group.size())
+                                  .cast<double>();
+    check_size_match("neg_binomial_2_log_likelihood_summary",
+                     "counts_per_group", counts.size(), "theta", theta_size);
+    check_size_match("neg_binomial_2_log_likelihood_summary", "n_per_group",
+                     n.size(), "theta", theta_size);
+    s = counts + eta * n;
+    x = add(theta, mean) - Eigen::VectorXd::Constant(theta_size, std::log(eta));
+  }
+
+  template <typename Theta, typename Eta, typename Mean>
+  static Eigen::VectorXd theta_grad(const Theta& theta, const Eta& eta,
+                                    const std::vector<int>& y,
+                                    const std::vector<int>& n_per_group,
+                                    const std::vector<int>& counts_per_group,
+                                    const Mean& mean) {
+    Eigen::VectorXd counts, s, x;
+    parts(theta, eta, n_per_group, counts_per_group, mean, counts, s, x);
+    return analytic_derivatives<neg_binomial_2_log_likelihood>::grad(counts, s,
+                                                                     x);
+  }
+
+  template <typename Theta, typename Eta, typename Mean>
+  static Eigen::VectorXd neg_hessian_diag(
+      const Theta& theta, const Eta& eta, const std::vector<int>& y,
+      const std::vector<int>& n_per_group,
+      const std::vector<int>& counts_per_group, const Mean& mean) {
+    Eigen::VectorXd counts, s, x;
+    parts(theta, eta, n_per_group, counts_per_group, mean, counts, s, x);
+    return analytic_derivatives<neg_binomial_2_log_likelihood>::neg_hess(s, x);
+  }
+
+  template <typename Theta, typename Eta, typename Mean>
+  static Eigen::VectorXd third_diff(const Theta& theta, const Eta& eta,
+                                    const std::vector<int>& y,
+                                    const std::vector<int>& n_per_group,
+                                    const std::vector<int>& counts_per_group,
+                                    const Mean& mean) {
+    Eigen::VectorXd counts, s, x;
+    parts(theta, eta, n_per_group, counts_per_group, mean, counts, s, x);
+    return analytic_derivatives<neg_binomial_2_log_likelihood>::third(s, x);
+  }
+};
+}  // namespace laplace_likelihood
+
 /**
  * Wrapper function around the laplace_marginal function for
  * a negative binomial likelihood. Uses the 2nd parameterization.
diff --git a/stan/math/mix/prob/laplace_marginal_poisson_log_lpmf.hpp b/stan/math/mix/prob/laplace_marginal_poisson_log_lpmf.hpp
index c9cb99fa629..7bc279a9291 100644
--- a/stan/math/mix/prob/laplace_marginal_poisson_log_lpmf.hpp
+++ b/stan/math/mix/prob/laplace_marginal_poisson_log_lpmf.hpp
@@ -2,6 +2,7 @@
 #define STAN_MATH_MIX_PROB_LAPLACE_MARGINAL_POISSON_LOG_LPMF_HPP
 
 #include <stan/math/mix/functor/laplace_likelihood.hpp>
+#include <type_traits>
 #include <stan/math/mix/functor/laplace_marginal_density.hpp>
 #include <stan/math/mix/meta.hpp>
 #include <stan/math/rev/fun/dot_product.hpp>
@@ -53,6 +54,68 @@ struct poisson_log_likelihood {
   }
 };
 
+namespace laplace_likelihood {
+/**
+ * Closed-form derivatives of `poisson_log_likelihood` with respect to
+ * `theta` (see `analytic_derivatives`).  With `m_g = n_g exp(theta_g +
+ * mean_g)`, where `n_g` is the number of observations of group `g`, the
+ * gradient is `counts_g - m_g`, the negative Hessian is `m_g`, and the third
+ * derivative is `-m_g`.
+ */
+template <>
+struct analytic_derivatives<poisson_log_likelihood> : std::true_type {
+  // sum of the counts and number of observations of each group
+  template <typename YVec>
+  static void group_sums(Eigen::Index theta_size, const YVec& y,
+                         const std::vector<int>& y_index,
+                         Eigen::VectorXd& counts, Eigen::VectorXd& n) {
+    counts = Eigen::VectorXd::Zero(theta_size);
+    n = Eigen::VectorXd::Zero(theta_size);
+    for (size_t i = 0; i < y_index.size(); i++) {
+      counts(y_index[i] - 1) += y[i];
+      n(y_index[i] - 1) += 1;
+    }
+  }
+
+  template <typename Theta, typename YVec, typename Mean>
+  static void rates(const Theta& theta, const YVec& y,
+                    const std::vector<int>& y_index, const Mean& mean,
+                    Eigen::VectorXd& counts, Eigen::VectorXd& m) {
+    Eigen::VectorXd n;
+    group_sums(theta.size(), y, y_index, counts, n);
+    Eigen::VectorXd theta_offset = add(theta, mean);
+    m = n.cwiseProduct(theta_offset.array().exp().matrix());
+  }
+
+  template <typename Theta, typename YVec, typename Mean>
+  static Eigen::VectorXd theta_grad(const Theta& theta, const YVec& y,
+                                    const std::vector<int>& y_index,
+                                    const Mean& mean) {
+    Eigen::VectorXd counts, m;
+    rates(theta, y, y_index, mean, counts, m);
+    return counts - m;
+  }
+
+  template <typename Theta, typename YVec, typename Mean>
+  static Eigen::VectorXd neg_hessian_diag(const Theta& theta, const YVec& y,
+                                          const std::vector<int>& y_index,
+                                          const Mean& mean) {
+    Eigen::VectorXd counts, m;
+    rates(theta, y, y_index, mean, counts, m);
+    return m;
+  }
+
+  template <typename Theta, typename YVec, typename Mean>
+  static Eigen::VectorXd third_diff(const Theta& theta, const YVec& y,
+                                    const std::vector<int>& y_index,
+                                    const Mean& mean) {
+    Eigen::VectorXd counts, m;
+    rates(theta, y, y_index, mean, counts, m);
+    return -m;
+  }
+};
+}  // namespace laplace_likelihood
+
 /**
  * Wrapper function around the laplace_marginal function for
  * a log poisson likelihood. Returns the marginal density
diff --git a/test/unit/math/laplace/laplace_analytic_derivatives_test.cpp b/test/unit/math/laplace/laplace_analytic_derivatives_test.cpp
new file mode 100644
index 00000000000..5ed15dcca23
--- /dev/null
+++ b/test/unit/math/laplace/laplace_analytic_derivatives_test.cpp
@@ -0,0 +1,268 @@
+#include <stan/math/mix.hpp>
+#include <gtest/gtest.h>
+#include <cmath>
+#include <ostream>
+#include <tuple>
+#include <type_traits>
+#include <vector>
+
+namespace {
+
+using stan::math::bernoulli_logit_likelihood;
+using stan::math::neg_binomial_2_log_likelihood;
+using stan::math::neg_binomial_2_log_likelihood_summary;
+using stan::math::poisson_log_likelihood;
+using stan::math::laplace_likelihood::analytic_derivatives;
+namespace lli = stan::math::laplace_likelihood::internal;
+
+// The likelihood `LL` behind a lambda: no `analytic_derivatives`
+// specialization, so the derivative functions use nested autodiff.
+template <typename LL>
+auto autodiff_version() {
+  return [](auto&&... args) { return LL{}(args...); };
+}
+
+void expect_near_rel(const Eigen::VectorXd& expected,
+                     const Eigen::VectorXd& actual, double tol,
+                     const char* what) {
+  ASSERT_EQ(expected.size(), actual.size()) << what;
+  for (Eigen::Index i = 0; i < expected.size(); ++i) {
+    EXPECT_TRUE(std::isfinite(actual(i))) << what << " [" << i << "]";
+    EXPECT_NEAR(expected(i), actual(i),
+                tol * std::fmax(1.0, std::fabs(expected(i))))
+        << what << " [" << i << "]";
+  }
+}
+
+void expect_same_sparse(const Eigen::SparseMatrix<double>& expected,
+                        const Eigen::SparseMatrix<double>& actual, double tol,
+                        const char* what) {
+  ASSERT_EQ(expected.rows(), actual.rows()) << what;
+  ASSERT_EQ(expected.cols(), actual.cols()) << what;
+  EXPECT_EQ(expected.nonZeros(), actual.nonZeros()) << what;
+  const Eigen::MatrixXd e = expected;
+  const Eigen::MatrixXd a = actual;
+  for (Eigen::Index c = 0; c < e.cols(); ++c) {
+    for (Eigen::Index r = 0; r < e.rows(); ++r) {
+      EXPECT_NEAR(e(r, c), a(r, c), tol * std::fmax(1.0, std::fabs(e(r, c))))
+          << what << " (" << r << ", " << c << ")";
+    }
+  }
+}
+
+// Each derivative function of laplace_likelihood::internal, for the
+// analytic version (`LL{}`) and for the autodiff version, with
+// hessian_block_size 1 and 2.
+template <typename LL, typename... Args>
+void expect_analytic_equals_autodiff(const Eigen::VectorXd& theta, double tol,
+                                     const Args&... args) {
+  static_assert(lli::use_analytic_derivatives_v<LL, const Eigen::VectorXd&,
+                                                const Args&...>,
+                "the analytic derivatives must be selected");
+  auto ad = autodiff_version<LL>();
+  std::ostream* msgs = nullptr;
+  expect_near_rel(lli::theta_grad(ad, theta, msgs, args...),
+                  lli::theta_grad(LL{}, theta, msgs, args...), tol,
+                  "theta_grad");
+  expect_near_rel(lli::diagonal_hessian(ad, theta, msgs, args...),
+                  lli::diagonal_hessian(LL{}, theta, msgs, args...), tol,
+                  "diagonal_hessian");
+  for (int block_size : {1, 2}) {
+    expect_same_sparse(
+        lli::block_hessian(ad, theta, block_size, msgs, args...),
+        lli::block_hessian(LL{}, theta, block_size, msgs, args...), tol,
+        "block_hessian");
+  }
+  expect_near_rel(lli::third_diff(ad, theta, msgs, args...),
+                  lli::third_diff(LL{}, theta, msgs, args...), tol,
+                  "third_diff");
+}
+
+// 4 groups, 1 to 3 observations per group
+const std::vector<int> y_counts{0, 3, 7, 2, 1, 4, 12, 9};
+const std::vector<int> y_binary{0, 1, 1, 0, 1, 1, 0, 1};
+const std::vector<int> y_index{1, 1, 1, 2, 3, 3, 4, 4};
+
+// a likelihood with the formulas of poisson_log_likelihood that counts the
+// calls of its analytic derivatives
+struct counted_poisson_likelihood {
+  template <typename... Args>
+  auto operator()(Args&&... args) const {
+    return poisson_log_likelihood{}(std::forward<Args>(args)...);
+  }
+};
+int n_analytic_calls = 0;
+
+}  // namespace
+
+namespace stan {
+namespace math {
+namespace laplace_likelihood {
+template <>
+struct analytic_derivatives<counted_poisson_likelihood> : std::true_type {
+  template <typename... Args>
+  static Eigen::VectorXd theta_grad(const Args&... args) {
+    ++n_analytic_calls;
+    return analytic_derivatives<poisson_log_likelihood>::theta_grad(args...);
+  }
+  template <typename... Args>
+  static Eigen::VectorXd neg_hessian_diag(const Args&... args) {
+    ++n_analytic_calls;
+    return analytic_derivatives<poisson_log_likelihood>::neg_hessian_diag(
+        args...);
+  }
+  template <typename... Args>
+  static Eigen::VectorXd third_diff(const Args&... args) {
+    ++n_analytic_calls;
+    return analytic_derivatives<poisson_log_likelihood>::third_diff(args...);
+  }
+};
+}  // namespace laplace_likelihood
+}  // namespace math
+}  // namespace stan
+
+namespace {
+
+struct se_covariance {
+  template <typename T1, typename T2>
+  auto operator()(const std::vector<double>& x, const T1& sigma, const T2& rho,
+                  std::ostream*) const {
+    return stan::math::add_diag(stan::math::gp_exp_quad_cov(x, sigma, rho),
+                                1e-8);
+  }
+};
+const std::vector<double> x_gp{0.0, 1.3, 2.1, 3.7};
+
+}  // namespace
+
+TEST(laplace_analytic_derivatives, trait) {
+  EXPECT_TRUE(analytic_derivatives<poisson_log_likelihood>::value);
+  EXPECT_TRUE(analytic_derivatives<bernoulli_logit_likelihood>::value);
+  EXPECT_TRUE(analytic_derivatives<neg_binomial_2_log_likelihood>::value);
+  EXPECT_TRUE(
+      analytic_derivatives<neg_binomial_2_log_likelihood_summary>::value);
+  auto ad = autodiff_version<poisson_log_likelihood>();
+  EXPECT_FALSE(analytic_derivatives<decltype(ad)>::value);
+  // a var argument selects the autodiff version
+  EXPECT_FALSE((lli::use_analytic_derivatives_v<
+                neg_binomial_2_log_likelihood, Eigen::VectorXd, stan::math::var,
+                std::vector<int>, std::vector<int>, Eigen::VectorXd>));
+}
+
+TEST(laplace_analytic_derivatives, used_by_laplace_marginal) {
+  const Eigen::VectorXd mean = Eigen::VectorXd::Constant(4, 0.5);
+  std::ostream* msgs = nullptr;
+  n_analytic_calls = 0;
+  double lp = stan::math::laplace_marginal<false>(
+      counted_poisson_likelihood{},
+      std::forward_as_tuple(y_counts, y_index, mean), 1, se_covariance{},
+      std::forward_as_tuple(x_gp, 1.0, 1.5), msgs);
+  EXPECT_TRUE(std::isfinite(lp));
+  EXPECT_GT(n_analytic_calls, 0);
+}
+
+TEST(laplace_analytic_derivatives, equal_to_autodiff) {
+  Eigen::VectorXd theta(4);
+  theta << 0.3, -1.2, 1.7, 0.05;
+  const Eigen::VectorXd mean = Eigen::VectorXd::Constant(4, 0.4);
+  const double tol = 1e-12;
+  for (double shift : {0.0, -30.0, 30.0}) {
+    const Eigen::VectorXd th = theta.array() + shift;
+    expect_analytic_equals_autodiff<poisson_log_likelihood>(th, tol, y_counts,
+                                                            y_index, mean);
+    expect_analytic_equals_autodiff<bernoulli_logit_likelihood>(
+        th, tol, y_binary, y_index, mean);
+    expect_analytic_equals_autodiff<neg_binomial_2_log_likelihood>(
+        th, tol, 2.5, y_counts, y_index, mean);
+    Eigen::VectorXd eta(4);
+    eta << 0.7, 2.5, 11.0, 1.3;
+    expect_analytic_equals_autodiff<neg_binomial_2_log_likelihood>(
+        th, tol, eta, y_counts, y_index, mean);
+    const std::vector<int> n_per_group{3, 1, 2, 2};
+    const std::vector<int> counts_per_group{10, 2, 5, 21};
+    expect_analytic_equals_autodiff<neg_binomial_2_log_likelihood_summary>(
+        th, tol, 2.5, y_counts, n_per_group, counts_per_group, mean);
+  }
+}
+
+TEST(laplace_analytic_derivatives, finite_at_large_theta) {
+  // exp(theta) overflows here; the inv_logit forms stay finite
+  Eigen::VectorXd theta(4);
+  theta << -800, 800, -750, 760;
+  const Eigen::VectorXd mean = Eigen::VectorXd::Zero(4);
+  expect_analytic_equals_autodiff<bernoulli_logit_likelihood>(
+      theta, 1e-12, y_binary, y_index, mean);
+  expect_analytic_equals_autodiff<neg_binomial_2_log_likelihood>(
+      theta, 1e-12, 2.5, y_counts, y_index, mean);
+}
+
+TEST(laplace_analytic_derivatives, end_to_end_equal_to_autodiff) {
+  using stan::math::var;
+  const Eigen::VectorXd mean = Eigen::VectorXd::Constant(4, 0.5);
+  std::ostream* msgs = nullptr;
+  const double tol = 1e-7;
+  for (int solver : {1, 2, 3}) {
+    for (int block_size : {1, 2}) {
+      auto ops = std::make_tuple(Eigen::VectorXd::Zero(4).eval(), 1e-10, 200,
+                                 solver, 1000, false);
+      // value and gradient w.r.t. sigma, rho (and eta) of `lpmf`
+      auto run = [&](auto&& lpmf) {
+        var sigma = 1.0, rho = 1.5, eta = 2.5;
+        var lp = lpmf(sigma, rho, eta);
+        stan::math::grad(lp.vi_);
+        std::vector<double> out{lp.val(), sigma.adj(), rho.adj(), eta.adj()};
+        stan::math::recover_memory();
+        return out;
+      };
+      auto compare = [&](const std::vector<double>& a,
+                         const std::vector<double>& b, const char* what) {
+        for (size_t i = 0; i < a.size(); ++i) {
+          EXPECT_NEAR(a[i], b[i], tol * std::fmax(1.0, std::fabs(a[i])))
+              << what << " solver " << solver << " block size " << block_size
+              << " output " << i;
+        }
+      };
+      compare(run([&](auto s, auto r, auto) {
+                return stan::math::laplace_marginal_tol<false>(
+                    autodiff_version<poisson_log_likelihood>(),
+                    std::forward_as_tuple(y_counts, y_index, mean), block_size,
+                    se_covariance{}, std::forward_as_tuple(x_gp, s, r), ops,
+                    msgs);
+              }),
+              run([&](auto s, auto r, auto) {
+                return stan::math::laplace_marginal_tol_poisson_log_lpmf(
+                    y_counts, y_index, mean, block_size, se_covariance{},
+                    std::forward_as_tuple(x_gp, s, r), ops, msgs);
+              }),
+              "poisson");
+      compare(run([&](auto s, auto r, auto) {
+                return stan::math::laplace_marginal_tol<false>(
+                    autodiff_version<bernoulli_logit_likelihood>(),
+                    std::forward_as_tuple(y_binary, y_index, mean), block_size,
+                    se_covariance{}, std::forward_as_tuple(x_gp, s, r), ops,
+                    msgs);
+              }),
+              run([&](auto s, auto r, auto) {
+                return stan::math::laplace_marginal_tol_bernoulli_logit_lpmf(
+                    y_binary, y_index, mean, block_size, se_covariance{},
+                    std::forward_as_tuple(x_gp, s, r), ops, msgs);
+              }),
+              "bernoulli");
+      // eta as a var: the Newton iterations use the analytic derivatives,
+      // the gradient w.r.t. eta the autodiff ones
+      compare(run([&](auto s, auto r, auto e) {
+                return stan::math::laplace_marginal_tol<false>(
+                    autodiff_version<neg_binomial_2_log_likelihood>(),
+                    std::forward_as_tuple(e, y_counts, y_index, mean),
+                    block_size, se_covariance{},
+                    std::forward_as_tuple(x_gp, s, r), ops, msgs);
+              }),
+              run([&](auto s, auto r, auto e) {
+                return stan::math::laplace_marginal_tol_neg_binomial_2_log_lpmf(
+                    y_counts, y_index, e, mean, block_size, se_covariance{},
+                    std::forward_as_tuple(x_gp, s, r), ops, msgs);
+              }),
+              "neg_binomial_2 (eta var)");
+    }
+  }
+}

From d7d750235b4ff79ef996a714b56ee2e8769bcc92 Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Mon, 19 Oct 2026 16:54:49 +0300
Subject: [PATCH 2/3] fix cpplint whitespace/braces in the analytic-derivatives
 test

---
 test/unit/math/laplace/laplace_analytic_derivatives_test.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/unit/math/laplace/laplace_analytic_derivatives_test.cpp b/test/unit/math/laplace/laplace_analytic_derivatives_test.cpp
index 5ed15dcca23..1e7c81558aa 100644
--- a/test/unit/math/laplace/laplace_analytic_derivatives_test.cpp
+++ b/test/unit/math/laplace/laplace_analytic_derivatives_test.cpp
@@ -88,7 +88,7 @@ const std::vector<int> y_index{1, 1, 1, 2, 3, 3, 4, 4};
 struct counted_poisson_likelihood {
   template <typename... Args>
   auto operator()(Args&&... args) const {
-    return poisson_log_likelihood{}(std::forward<Args>(args)...);
+    return poisson_log_likelihood()(std::forward<Args>(args)...);
   }
 };
 int n_analytic_calls = 0;

From 31cd42671202290ad215a46ae2ad35b49789ef91 Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Mon, 26 Oct 2026 23:48:49 +0200
Subject: [PATCH 3/3] share group_sums, test tail values against references,
 remove unused internal::diff

---
 stan/math/mix/functor/laplace_likelihood.hpp  | 103 ++++--------------
 .../laplace_marginal_bernoulli_logit_lpmf.hpp |  18 +--
 ...place_marginal_neg_binomial_2_log_lpmf.hpp |  18 +--
 .../laplace_marginal_poisson_log_lpmf.hpp     |  18 +--
 .../laplace_analytic_derivatives_test.cpp     |  41 +++++++
 5 files changed, 76 insertions(+), 122 deletions(-)

diff --git a/stan/math/mix/functor/laplace_likelihood.hpp b/stan/math/mix/functor/laplace_likelihood.hpp
index 75b7d23f468..07958076a5c 100644
--- a/stan/math/mix/functor/laplace_likelihood.hpp
+++ b/stan/math/mix/functor/laplace_likelihood.hpp
@@ -12,6 +12,7 @@
 #include <tuple>
 #include <type_traits>
 #include <utility>
+#include <vector>
 
 namespace stan {
 namespace math {
@@ -43,6 +44,28 @@ struct analytic_derivatives : std::false_type {};
 
 namespace internal {
 
+/**
+ * Sum of the observations and number of observations of each group, for the
+ * closed-form derivatives of the built-in likelihoods.
+ * @tparam YVec Type of the observations
+ * @param[in] theta_size number of groups
+ * @param[in] y observations
+ * @param[in] y_index group of each observation (1-based)
+ * @param[out] counts sum of `y` in each group
+ * @param[out] n number of observations in each group
+ */
+template <typename YVec>
+inline void group_sums(Eigen::Index theta_size, const YVec& y,
+                       const std::vector<int>& y_index, Eigen::VectorXd& counts,
+                       Eigen::VectorXd& n) {
+  counts = Eigen::VectorXd::Zero(theta_size);
+  n = Eigen::VectorXd::Zero(theta_size);
+  for (size_t i = 0; i < y_index.size(); i++) {
+    counts(y_index[i] - 1) += y[i];
+    n(y_index[i] - 1) += 1;
+  }
+}
+
 /**
  * True if the closed-form derivatives of `F` (see `analytic_derivatives`)
  * are used for these argument types: `F` has them, and `theta` and all
@@ -202,57 +225,6 @@ inline auto block_hessian(F&& f, Theta&& theta,
   }
 }
 
-/**
- * Computes theta gradient and negative block diagonal Hessian of `f` wrt
- * `theta` and `args...`
- * @note If `Args` contains \ref var types then their adjoints will be
- * calculated as a side effect.
- * @tparam F A functor with `opertor()(Args&&...)` returning a scalar
- * @tparam Theta A class assignable to an Eigen vector type
- * @tparam Stream Type of stream for messages.
- * @tparam Args Type of variadic arguments.
- * @param f Log likelihood function.
- * @param theta Latent Gaussian model.
- * @param hessian_block_size If the Hessian of the log likelihood function w.r.t
- *                           the latent Gaussian variable is block-diagonal,
- *                           size of each block.
- * @param msgs Stream for messages.
- * @param args Variadic arguments for the likelihood function.
- */
-template <typename F, typename Theta, typename Stream, typename... Args,
-          require_eigen_vector_vt<std::is_arithmetic, Theta>* = nullptr>
-inline auto diff(F&& f, Theta&& theta, const Eigen::Index hessian_block_size,
-                 Stream* msgs, Args&&... args) {
-  using Eigen::Dynamic;
-  using Eigen::Matrix;
-  const Eigen::Index theta_size = theta.size();
-  auto theta_gradient = [&theta, &f, &msgs](auto&&... args) {
-    nested_rev_autodiff nested;
-    Matrix<var, Dynamic, 1> theta_var = theta;
-    var f_var = f(theta_var, args..., msgs);
-    grad(f_var.vi_);
-    return theta_var.adj().eval();
-  }(args...);
-  if (hessian_block_size == 1) {
-    auto v = Eigen::VectorXd::Ones(theta_size);
-    Eigen::VectorXd hessian_v = Eigen::VectorXd::Zero(theta_size);
-    hessian_times_vector(f, hessian_v, std::forward<Theta>(theta), std::move(v),
-                         value_of(args)..., msgs);
-    Eigen::SparseMatrix<double> hessian_theta(theta_size, theta_size);
-    hessian_theta.reserve(Eigen::VectorXi::Constant(theta_size, 1));
-    for (Eigen::Index i = 0; i < theta_size; i++) {
-      hessian_theta.insert(i, i) = hessian_v(i);
-    }
-    return std::make_pair(std::move(theta_gradient), (-hessian_theta).eval());
-  } else {
-    return std::make_pair(
-        std::move(theta_gradient),
-        (-hessian_block_diag(f, std::forward<Theta>(theta), hessian_block_size,
-                             value_of(args)..., msgs))
-            .eval());
-  }
-}
-
 /**
  * Compute third order derivative of `f` wrt `theta` and `args...`
  * @note If `Args` contains \ref var types then their adjoints will be
@@ -570,35 +542,6 @@ inline auto log_likelihood(F&& f, Theta&& theta, TupleArgs&& ll_tup,
       std::forward<Theta>(theta), msgs);
 }
 
-/**
- * A wrapper that accepts a tuple as arguments.
- * @tparam F A functor with `opertor()(Args&&...)` returning a scalar
- * @tparam Theta A class assignable to an Eigen vector type
- * @tparam TupleArgs Type of arguments for covariance function.
- * @tparam Stream Type of stream for messages.
- * @param f Log likelihood function.
- * @param theta Latent Gaussian model.
- * @param hessian_block_size If Hessian of log likelihood w.r.t theta is
- *                           block diagonal, size of block.
- * @param ll_tuple Arguments for likelihood function
- * @param msgs Stream messages.
- */
-template <typename F, typename Theta, typename TupleArgs, typename Stream,
-          require_eigen_vector_t<Theta>* = nullptr,
-          require_tuple_t<TupleArgs>* = nullptr>
-inline auto diff(F&& f, Theta&& theta, const Eigen::Index hessian_block_size,
-                 TupleArgs&& ll_tuple, Stream* msgs) {
-  return stan::math::apply(
-      [](auto&& f, auto&& theta, auto hessian_block_size, auto* msgs,
-         auto&&... args) {
-        return internal::diff(
-            std::forward<decltype(f)>(f), std::forward<decltype(theta)>(theta),
-            hessian_block_size, msgs, std::forward<decltype(args)>(args)...);
-      },
-      std::forward<TupleArgs>(ll_tuple), std::forward<F>(f),
-      std::forward<Theta>(theta), hessian_block_size, msgs);
-}
-
 /**
  * A wrapper that accepts a tuple as arguments.
  * @tparam F Type of log likelhood function.
diff --git a/stan/math/mix/prob/laplace_marginal_bernoulli_logit_lpmf.hpp b/stan/math/mix/prob/laplace_marginal_bernoulli_logit_lpmf.hpp
index 78e809988c0..e2a85fc9b8d 100644
--- a/stan/math/mix/prob/laplace_marginal_bernoulli_logit_lpmf.hpp
+++ b/stan/math/mix/prob/laplace_marginal_bernoulli_logit_lpmf.hpp
@@ -4,7 +4,6 @@
 #include <stan/math/mix/meta.hpp>
 #include <stan/math/mix/functor/laplace_likelihood.hpp>
 #include <stan/math/prim/fun/inv_logit.hpp>
-#include <type_traits>
 #include <stan/math/mix/functor/laplace_marginal_density.hpp>
 #include <stan/math/rev/fun/elt_multiply.hpp>
 #include <stan/math/rev/fun/exp.hpp>
@@ -22,6 +21,8 @@
 #include <stan/math/prim/fun/binomial_coefficient_log.hpp>
 #include <stan/math/rev/fun/to_vector.hpp>
 #include <stan/math/prim/fun/to_vector.hpp>
+#include <type_traits>
+#include <vector>
 
 namespace stan {
 namespace math {
@@ -71,25 +72,12 @@ namespace laplace_likelihood {
  */
 template <>
 struct analytic_derivatives<bernoulli_logit_likelihood> : std::true_type {
-  // sum of the counts and number of observations of each group
-  template <typename YVec>
-  static void group_sums(Eigen::Index theta_size, const YVec& y,
-                         const std::vector<int>& y_index,
-                         Eigen::VectorXd& counts, Eigen::VectorXd& n) {
-    counts = Eigen::VectorXd::Zero(theta_size);
-    n = Eigen::VectorXd::Zero(theta_size);
-    for (size_t i = 0; i < y_index.size(); i++) {
-      counts(y_index[i] - 1) += y[i];
-      n(y_index[i] - 1) += 1;
-    }
-  }
-
   template <typename Theta, typename YVec, typename Mean>
   static void probabilities(const Theta& theta, const YVec& y,
                             const std::vector<int>& y_index, const Mean& mean,
                             Eigen::VectorXd& counts, Eigen::VectorXd& n,
                             Eigen::VectorXd& p, Eigen::VectorXd& q) {
-    group_sums(theta.size(), y, y_index, counts, n);
+    internal::group_sums(theta.size(), y, y_index, counts, n);
     Eigen::VectorXd theta_offset = add(theta, mean);
     p = inv_logit(theta_offset);
     q = inv_logit(-theta_offset);
diff --git a/stan/math/mix/prob/laplace_marginal_neg_binomial_2_log_lpmf.hpp b/stan/math/mix/prob/laplace_marginal_neg_binomial_2_log_lpmf.hpp
index bc56d3444e3..7e9232f1688 100644
--- a/stan/math/mix/prob/laplace_marginal_neg_binomial_2_log_lpmf.hpp
+++ b/stan/math/mix/prob/laplace_marginal_neg_binomial_2_log_lpmf.hpp
@@ -3,8 +3,6 @@
 
 #include <stan/math/mix/functor/laplace_likelihood.hpp>
 #include <stan/math/prim/fun/inv_logit.hpp>
-#include <cmath>
-#include <type_traits>
 #include <stan/math/mix/functor/laplace_marginal_density.hpp>
 
 #include <stan/math/rev/core/operator_addition.hpp>
@@ -25,6 +23,8 @@
 #include <stan/math/fwd/fun/sum.hpp>
 #include <stan/math/prim/err/check_size_match.hpp>
 #include <stan/math/prim/fun/binomial_coefficient_log.hpp>
+#include <type_traits>
+#include <vector>
 
 namespace stan {
 namespace math {
@@ -232,21 +232,15 @@ struct analytic_derivatives<neg_binomial_2_log_likelihood> : std::true_type {
                     const Mean& mean, Eigen::VectorXd& counts,
                     Eigen::VectorXd& s, Eigen::VectorXd& x) {
     const Eigen::Index theta_size = theta.size();
-    counts = Eigen::VectorXd::Zero(theta_size);
-    Eigen::VectorXd n = Eigen::VectorXd::Zero(theta_size);
-    for (size_t i = 0; i < y.size(); i++) {
-      counts(y_index[i] - 1) += y[i];
-      n(y_index[i] - 1) += 1;
-    }
+    Eigen::VectorXd n;
+    internal::group_sums(theta_size, y, y_index, counts, n);
     Eigen::VectorXd eta_g(theta_size);
     if constexpr (is_stan_scalar<Eta>::value) {
       eta_g.setConstant(eta);
     } else {
       check_size_match("neg_binomial_2_log_likelihood", "eta", eta.size(),
                        "theta", theta_size);
-      for (Eigen::Index g = 0; g < theta_size; ++g) {
-        eta_g(g) = eta[g];
-      }
+      eta_g = eta;
     }
     s = counts + eta_g.cwiseProduct(n);
     x = add(theta, mean) - eta_g.array().log().matrix();
@@ -309,7 +303,7 @@ struct analytic_derivatives<neg_binomial_2_log_likelihood_summary>
     check_size_match("neg_binomial_2_log_likelihood_summary", "n_per_group",
                      n.size(), "theta", theta_size);
     s = counts + eta * n;
-    x = add(theta, mean) - Eigen::VectorXd::Constant(theta_size, std::log(eta));
+    x = (add(theta, mean).array() - log(eta)).matrix();
   }
 
   template <typename Theta, typename Eta, typename Mean>
diff --git a/stan/math/mix/prob/laplace_marginal_poisson_log_lpmf.hpp b/stan/math/mix/prob/laplace_marginal_poisson_log_lpmf.hpp
index 7bc279a9291..0ac2f43f0fa 100644
--- a/stan/math/mix/prob/laplace_marginal_poisson_log_lpmf.hpp
+++ b/stan/math/mix/prob/laplace_marginal_poisson_log_lpmf.hpp
@@ -2,7 +2,6 @@
 #define STAN_MATH_MIX_PROB_LAPLACE_MARGINAL_POISSON_LOG_LPMF_HPP
 
 #include <stan/math/mix/functor/laplace_likelihood.hpp>
-#include <type_traits>
 #include <stan/math/mix/functor/laplace_marginal_density.hpp>
 #include <stan/math/mix/meta.hpp>
 #include <stan/math/rev/fun/dot_product.hpp>
@@ -13,6 +12,8 @@
 #include <stan/math/fwd/fun/lgamma.hpp>
 #include <stan/math/rev/fun/lgamma.hpp>
 #include <stan/math/prim/fun/add.hpp>
+#include <type_traits>
+#include <vector>
 
 namespace stan {
 namespace math {
@@ -64,25 +65,12 @@ namespace laplace_likelihood {
  */
 template <>
 struct analytic_derivatives<poisson_log_likelihood> : std::true_type {
-  // sum of the counts and number of observations of each group
-  template <typename YVec>
-  static void group_sums(Eigen::Index theta_size, const YVec& y,
-                         const std::vector<int>& y_index,
-                         Eigen::VectorXd& counts, Eigen::VectorXd& n) {
-    counts = Eigen::VectorXd::Zero(theta_size);
-    n = Eigen::VectorXd::Zero(theta_size);
-    for (size_t i = 0; i < y_index.size(); i++) {
-      counts(y_index[i] - 1) += y[i];
-      n(y_index[i] - 1) += 1;
-    }
-  }
-
   template <typename Theta, typename YVec, typename Mean>
   static void rates(const Theta& theta, const YVec& y,
                     const std::vector<int>& y_index, const Mean& mean,
                     Eigen::VectorXd& counts, Eigen::VectorXd& m) {
     Eigen::VectorXd n;
-    group_sums(theta.size(), y, y_index, counts, n);
+    internal::group_sums(theta.size(), y, y_index, counts, n);
     Eigen::VectorXd theta_offset = add(theta, mean);
     m = n.cwiseProduct(theta_offset.array().exp().matrix());
   }
diff --git a/test/unit/math/laplace/laplace_analytic_derivatives_test.cpp b/test/unit/math/laplace/laplace_analytic_derivatives_test.cpp
index 1e7c81558aa..f10b785d641 100644
--- a/test/unit/math/laplace/laplace_analytic_derivatives_test.cpp
+++ b/test/unit/math/laplace/laplace_analytic_derivatives_test.cpp
@@ -22,6 +22,8 @@ auto autodiff_version() {
   return [](auto&&... args) { return LL{}(args...); };
 }
 
+// relative tolerance for entries above 1 in absolute value,
+// absolute tolerance below 1
 void expect_near_rel(const Eigen::VectorXd& expected,
                      const Eigen::VectorXd& actual, double tol,
                      const char* what) {
@@ -185,6 +187,45 @@ TEST(laplace_analytic_derivatives, equal_to_autodiff) {
   }
 }
 
+// The shifted cases of equal_to_autodiff compare entries of about 1e-13 with
+// an absolute tolerance of 1e-12, so they do not see the cancellation in
+// `1 - p`.  Reference values (mpmath, 50 digits, from the double values of
+// `theta + mean` and `theta + mean - log(eta)`) check the stable forms.
+TEST(laplace_analytic_derivatives, tail_reference_values) {
+  Eigen::VectorXd theta(4);
+  theta << 30.3, 28.8, 31.7, 30.05;
+  const Eigen::VectorXd mean = Eigen::VectorXd::Constant(4, 0.4);
+  std::ostream* msgs = nullptr;
+  auto expect_rel = [](const std::vector<double>& ref,
+                       const Eigen::VectorXd& actual, const char* what) {
+    ASSERT_EQ(static_cast<Eigen::Index>(ref.size()), actual.size()) << what;
+    for (int i = 0; i < 4; ++i) {
+      EXPECT_NEAR(ref[i], actual(i), 1e-14 * std::fabs(ref[i]))
+          << what << " [" << i << "]";
+    }
+  };
+  expect_rel({1.3940574134239623e-13, 2.0825772910546344e-13,
+              2.2918021714044091e-14, 1.1933367674388607e-13},
+             lli::diagonal_hessian(bernoulli_logit_likelihood{}, theta, msgs,
+                                   y_binary, y_index, mean),
+             "bernoulli neg_hessian");
+  expect_rel({1.3940574134238327e-13, 2.082577291053767e-13,
+              2.2918021714043565e-14, 1.1933367674387183e-13},
+             lli::third_diff(bernoulli_logit_likelihood{}, theta, msgs,
+                             y_binary, y_index, mean),
+             "bernoulli third_diff");
+  expect_rel({2.0330003945763301e-12, 2.3428994524350021e-12,
+              2.8647527142554155e-13, 3.8783444941756066e-12},
+             lli::diagonal_hessian(neg_binomial_2_log_likelihood{}, theta, msgs,
+                                   2.5, y_counts, y_index, mean),
+             "neg_binomial_2 neg_hessian");
+  expect_rel({2.0330003945758578e-12, 2.3428994524325625e-12,
+              2.8647527142552514e-13, 3.8783444941744496e-12},
+             lli::third_diff(neg_binomial_2_log_likelihood{}, theta, msgs, 2.5,
+                             y_counts, y_index, mean),
+             "neg_binomial_2 third_diff");
+}
+
 TEST(laplace_analytic_derivatives, finite_at_large_theta) {
   // exp(theta) overflows here; the inv_logit forms stay finite
   Eigen::VectorXd theta(4);