From dcc52403bf33d08c4fb260585b6d2054231509a1 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 20 Sep 2023 09:06:52 +0300 Subject: [PATCH 1/6] Add implementation and tests --- stan/math/prim/prob.hpp | 1 + .../prim/prob/binomial_logit_glm_lpmf.hpp | 174 ++++++++++++++++++ .../mix/prob/binomial_logit_glm_lpmf_test.cpp | 35 ++++ .../prob/binomial_logit_glm_lpmf_test.cpp | 82 +++++++++ 4 files changed, 292 insertions(+) create mode 100644 stan/math/prim/prob/binomial_logit_glm_lpmf.hpp create mode 100644 test/unit/math/mix/prob/binomial_logit_glm_lpmf_test.cpp create mode 100644 test/unit/math/prim/prob/binomial_logit_glm_lpmf_test.cpp diff --git a/stan/math/prim/prob.hpp b/stan/math/prim/prob.hpp index bfd8cfdc07e..3af4c5763f6 100644 --- a/stan/math/prim/prob.hpp +++ b/stan/math/prim/prob.hpp @@ -46,6 +46,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/prob/binomial_logit_glm_lpmf.hpp b/stan/math/prim/prob/binomial_logit_glm_lpmf.hpp new file mode 100644 index 00000000000..c46a17d998c --- /dev/null +++ b/stan/math/prim/prob/binomial_logit_glm_lpmf.hpp @@ -0,0 +1,174 @@ +#ifndef STAN_MATH_PRIM_PROB_BINOMIAL_LOGIT_GLM_LPMF_HPP +#define STAN_MATH_PRIM_PROB_BINOMIAL_LOGIT_GLM_LPMF_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** \ingroup multivar_dists + * Returns the log PMF of the Generalized Linear Model (GLM) + * with Binomial distribution and logit link function. + * The idea is that binomial_logit_glm_lpmf(n | N, x, alpha, beta) should + * compute a more efficient version of + * binomial_logit_lpmf(y | N, alpha + x * beta) by using analytically + * simplified gradients. + * If containers are supplied, returns the log sum of the probabilities. + * + * @tparam T_y type of binary vector of successes variables; + * this can also be a single binary value; + * @tparam T_y type of binary vector of population size variables; + * this can also be a single binary value; + * @tparam T_x type of the matrix of independent variables (features) + * @tparam T_alpha type of the intercept(s); + * this can be a vector (of the same length as y) of intercepts or a single + * value (for models with constant intercept); + * @tparam T_beta type of the weight vector + * + * @param n binary scalar or vector parameter. If it is a scalar it will be + * broadcast - used for all instances. + * @param N binary scalar or vector parameter. If it is a scalar it will be + * broadcast - used for all instances. + * @param x design matrix or row vector. If it is a row vector it will be + * broadcast - used for all instances. + * @param alpha intercept (in log odds) + * @param beta weight vector + * @return log probability or log sum of probabilities + * @throw std::domain_error if x, beta or alpha is infinite. + * @throw std::domain_error if n is negative or greater than N + * @throw std::domain_error if N is negative + * @throw std::invalid_argument if container sizes mismatch. + */ +template * = nullptr> +return_type_t binomial_logit_glm_lpmf( + const T_n& n, const T_N& N, const T_x& x, + const T_alpha& alpha, const T_beta& beta) { + constexpr int T_x_rows = T_x::RowsAtCompileTime; + using T_xbeta_partials = partials_return_t; + using T_partials_return = partials_return_t; + using T_theta_tmp = + typename std::conditional_t>; + using T_xbeta_tmp = + typename std::conditional_t>; + using T_n_ref = ref_type_if_t::value, T_n>; + using T_N_ref = ref_type_if_t::value, T_N>; + using T_x_ref = ref_type_if_t::value, T_x>; + using T_alpha_ref = ref_type_if_t::value, T_alpha>; + using T_beta_ref = ref_type_if_t::value, T_beta>; + + T_n_ref n_ref = n; + T_N_ref N_ref = N; + T_x_ref x_ref = x; + T_alpha_ref alpha_ref = alpha; + T_beta_ref beta_ref = beta; + + if (size_zero(n, N, alpha, beta, x)) { + return 0; + } + + if (!include_summand::value) { + return 0; + } + + const size_t N_instances = max_size(n, N, x.col(0), alpha); + const size_t N_attributes = x.cols(); + + static const char* function = "binomial_logit_glm_lpmf"; + check_consistent_sizes(function, "Successes variable", n, + "Population size parameter", N); + check_consistent_size(function, "Successes variable", n, N_instances); + check_consistent_size(function, "Population size parameter", N, N_instances); + check_consistent_size(function, "Weight vector", beta, N_attributes); + check_consistent_size(function, "Vector of intercepts", alpha, N_instances); + + const auto& n_val = as_value_column_array_or_scalar(n_ref); + const auto& N_val = as_value_column_array_or_scalar(N_ref); + + check_bounded(function, "Successes variable", n_val, 0, N_val); + check_nonnegative(function, "Population size parameter", N_val); + + const auto& alpha_val = as_value_column_array_or_scalar(alpha_ref); + const auto& beta_val = as_value_column_vector_or_scalar(beta_ref); + const auto& x_val = value_of(x_ref); + Eigen::Array theta(N_instances); + if (T_x_rows == 1) { + theta = forward_as((x_val * beta_val)(0, 0)) + alpha_val; + } else { + theta = (x_val * beta_val).array() + alpha_val; + } + + const auto& log_inv_logit_theta = log_inv_logit(theta); + const auto& log1m_inv_logit_theta = log1m_inv_logit(theta); + + T_partials_return logp = sum(n_val * log_inv_logit_theta + + (N_val - n_val) * log1m_inv_logit_theta); + + using std::isfinite; + if (!isfinite(logp)) { + check_finite(function, "Weight vector", beta); + check_finite(function, "Intercept", alpha); + check_finite(function, "Matrix of independent variables", x); + } + + if (include_summand::value) { + size_t broadcast_n = max_size(N, n) == N_instances ? 1 : N_instances; + logp += sum(binomial_coefficient_log(N_val, n_val)) * broadcast_n; + } + + auto ops_partials = make_partials_propagator(x_ref, alpha_ref, beta_ref); + + if (!is_constant_all::value) { + Eigen::Matrix theta_derivative + = n_val - N_val * exp(log_inv_logit_theta); + + if (!is_constant_all::value) { + if (T_x_rows == 1) { + edge<2>(ops_partials).partials_ + = forward_as>( + theta_derivative.sum() * x_val); + } else { + partials<2>(ops_partials) = x_val.transpose() * theta_derivative; + } + } + + if (!is_constant_all::value) { + if (T_x_rows == 1) { + edge<0>(ops_partials).partials_ + = forward_as>( + beta_val * theta_derivative.sum()); + } else { + edge<0>(ops_partials).partials_ + = (beta_val * theta_derivative.transpose()).transpose(); + } + } + if (!is_constant_all::value) { + partials<1>(ops_partials) = theta_derivative; + } + } + return ops_partials.build(logp); +} + +template +inline return_type_t binomial_logit_glm_lpmf( + const T_n& n, const T_N& N, const T_x& x, + const T_alpha& alpha, const T_beta& beta) { + return binomial_logit_glm_lpmf(n, N, x, alpha, beta); +} +} // namespace math +} // namespace stan +#endif diff --git a/test/unit/math/mix/prob/binomial_logit_glm_lpmf_test.cpp b/test/unit/math/mix/prob/binomial_logit_glm_lpmf_test.cpp new file mode 100644 index 00000000000..f71124b696a --- /dev/null +++ b/test/unit/math/mix/prob/binomial_logit_glm_lpmf_test.cpp @@ -0,0 +1,35 @@ +#include +#include + +TEST(mathMixScalFun, binomial_logit_glm_lpmf) { + auto f = [](const auto n, const auto N) { + return [=](const auto& x, const auto& alpha, const auto& beta) { + return stan::math::binomial_logit_glm_lpmf(n, N, x, alpha, beta); + }; + }; + + std::vector n_arr{1, 4}; + std::vector N_arr{10, 45}; + Eigen::MatrixXd x = Eigen::MatrixXd::Random(2, 2); + Eigen::RowVectorXd x_rowvec = x.row(0); + Eigen::VectorXd alpha = Eigen::VectorXd::Random(2); + Eigen::VectorXd beta = Eigen::VectorXd::Random(2); + + stan::test::expect_ad(f(n_arr[0], N_arr[0]), x, alpha, beta); + stan::test::expect_ad(f(n_arr[0], N_arr), x, alpha, beta); + stan::test::expect_ad(f(n_arr, N_arr[0]), x, alpha, beta); + stan::test::expect_ad(f(n_arr, N_arr), x, alpha, beta); + stan::test::expect_ad(f(n_arr[0], N_arr[0]), x, alpha[0], beta); + stan::test::expect_ad(f(n_arr[0], N_arr), x, alpha[0], beta); + stan::test::expect_ad(f(n_arr, N_arr[0]), x, alpha[0], beta); + stan::test::expect_ad(f(n_arr, N_arr), x, alpha[0], beta); + stan::test::expect_ad(f(n_arr[0], N_arr[0]), x_rowvec, alpha, beta); + stan::test::expect_ad(f(n_arr[0], N_arr), x_rowvec, alpha, beta); + stan::test::expect_ad(f(n_arr, N_arr[0]), x_rowvec, alpha, beta); + stan::test::expect_ad(f(n_arr, N_arr), x_rowvec, alpha, beta); + stan::test::expect_ad(f(n_arr[0], N_arr[0]), x_rowvec, alpha[0], beta); + stan::test::expect_ad(f(n_arr[0], N_arr), x_rowvec, alpha[0], beta); + stan::test::expect_ad(f(n_arr, N_arr[0]), x_rowvec, alpha[0], beta); + stan::test::expect_ad(f(n_arr, N_arr), x_rowvec, alpha[0], beta); + +} diff --git a/test/unit/math/prim/prob/binomial_logit_glm_lpmf_test.cpp b/test/unit/math/prim/prob/binomial_logit_glm_lpmf_test.cpp new file mode 100644 index 00000000000..f23e78b890c --- /dev/null +++ b/test/unit/math/prim/prob/binomial_logit_glm_lpmf_test.cpp @@ -0,0 +1,82 @@ +#include +#include + +TEST(ProbBinomialLogitGLM, matchesNonGLM) { + using stan::math::binomial_logit_lpmf; + using stan::math::binomial_logit_glm_lpmf; + + std::vector n{1, 2}; + std::vector N{5, 4}; + Eigen::MatrixXd x = Eigen::MatrixXd::Random(2, 2); + Eigen::RowVectorXd x_row = x.row(0); + Eigen::VectorXd alpha = Eigen::VectorXd::Random(2); + Eigen::VectorXd beta = Eigen::VectorXd::Random(2); + + Eigen::VectorXd theta = alpha + x * beta; + + EXPECT_FLOAT_EQ(binomial_logit_lpmf(n, N, theta), + binomial_logit_glm_lpmf(n, N, x, alpha, beta)); + EXPECT_FLOAT_EQ(binomial_logit_lpmf(n[0], N, theta), + binomial_logit_glm_lpmf(n[0], N, x, alpha, beta)); + EXPECT_FLOAT_EQ(binomial_logit_lpmf(n, N[0], theta), + binomial_logit_glm_lpmf(n, N[0], x, alpha, beta)); + EXPECT_FLOAT_EQ(binomial_logit_lpmf(n[0], N[0], theta), + binomial_logit_glm_lpmf(n[0], N[0], x, alpha, beta)); + + theta = (alpha[0] + (x * beta).array()).matrix(); + + EXPECT_FLOAT_EQ(binomial_logit_lpmf(n, N, theta), + binomial_logit_glm_lpmf(n, N, x, alpha[0], beta)); + EXPECT_FLOAT_EQ(binomial_logit_lpmf(n[0], N, theta), + binomial_logit_glm_lpmf(n[0], N, x, alpha[0], beta)); + EXPECT_FLOAT_EQ(binomial_logit_lpmf(n, N[0], theta), + binomial_logit_glm_lpmf(n, N[0], x, alpha[0], beta)); + EXPECT_FLOAT_EQ(binomial_logit_lpmf(n[0], N[0], theta), + binomial_logit_glm_lpmf(n[0], N[0], x, alpha[0], beta)); + + theta = (alpha.array() + (x_row * beta)(0, 0)).matrix(); + + EXPECT_FLOAT_EQ(binomial_logit_lpmf(n, N, theta), + binomial_logit_glm_lpmf(n, N, x_row, alpha, beta)); + EXPECT_FLOAT_EQ(binomial_logit_lpmf(n[0], N, theta), + binomial_logit_glm_lpmf(n[0], N, x_row, alpha, beta)); + EXPECT_FLOAT_EQ(binomial_logit_lpmf(n, N[0], theta), + binomial_logit_glm_lpmf(n, N[0], x_row, alpha, beta)); + EXPECT_FLOAT_EQ(binomial_logit_lpmf(n[0], N[0], theta), + binomial_logit_glm_lpmf(n[0], N[0], x_row, alpha, beta)); +} + +TEST(ProbBinomialLogitGLM, throwsCorrectly) { + using stan::math::binomial_logit_glm_lpmf; + using stan::math::INFTY; + + std::vector n{1, 2}; + std::vector N{5, 4}; + Eigen::MatrixXd x = Eigen::MatrixXd::Random(2, 2); + Eigen::VectorXd alpha = Eigen::VectorXd::Random(2); + Eigen::VectorXd beta = Eigen::VectorXd::Random(2); + + std::vector N_mismatch_size{5, 4, 10}; + EXPECT_THROW(binomial_logit_glm_lpmf(n, N_mismatch_size, x, alpha, beta), + std::invalid_argument); + EXPECT_THROW(binomial_logit_glm_lpmf(500, 1, x, alpha, beta), + std::domain_error); + EXPECT_THROW(binomial_logit_glm_lpmf(-10, N, x, alpha, beta), + std::domain_error); + EXPECT_THROW(binomial_logit_glm_lpmf(n, -10, x, alpha, beta), + std::domain_error); + + Eigen::VectorXd alpha_inf = alpha; + alpha[0] = INFTY; + Eigen::VectorXd beta_inf = beta; + beta[0] = INFTY; + Eigen::MatrixXd x_inf = x; + x(0, 0) = INFTY; + + EXPECT_THROW(binomial_logit_glm_lpmf(n, N, x_inf, alpha, beta), + std::domain_error); + EXPECT_THROW(binomial_logit_glm_lpmf(n, N, x, alpha_inf, beta), + std::domain_error); + EXPECT_THROW(binomial_logit_glm_lpmf(n, N, x, alpha, beta_inf), + std::domain_error); +} From 337ae79178490151ec640b4224d350f186beabb1 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 20 Sep 2023 09:13:08 +0300 Subject: [PATCH 2/6] Cpplint --- test/unit/math/mix/prob/binomial_logit_glm_lpmf_test.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/test/unit/math/mix/prob/binomial_logit_glm_lpmf_test.cpp b/test/unit/math/mix/prob/binomial_logit_glm_lpmf_test.cpp index f71124b696a..2d7d0d08b3a 100644 --- a/test/unit/math/mix/prob/binomial_logit_glm_lpmf_test.cpp +++ b/test/unit/math/mix/prob/binomial_logit_glm_lpmf_test.cpp @@ -31,5 +31,4 @@ TEST(mathMixScalFun, binomial_logit_glm_lpmf) { stan::test::expect_ad(f(n_arr[0], N_arr), x_rowvec, alpha[0], beta); stan::test::expect_ad(f(n_arr, N_arr[0]), x_rowvec, alpha[0], beta); stan::test::expect_ad(f(n_arr, N_arr), x_rowvec, alpha[0], beta); - } From f3625c34d40091c90e0d16cb397a0ea50e4f5d76 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Wed, 20 Sep 2023 02:29:12 -0400 Subject: [PATCH 3/6] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/prob/binomial_logit_glm_lpmf.hpp | 8 ++++---- .../prim/prob/binomial_logit_glm_lpmf_test.cpp | 16 ++++++++-------- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/stan/math/prim/prob/binomial_logit_glm_lpmf.hpp b/stan/math/prim/prob/binomial_logit_glm_lpmf.hpp index c46a17d998c..8d8d5a75308 100644 --- a/stan/math/prim/prob/binomial_logit_glm_lpmf.hpp +++ b/stan/math/prim/prob/binomial_logit_glm_lpmf.hpp @@ -53,8 +53,8 @@ namespace math { template * = nullptr> return_type_t binomial_logit_glm_lpmf( - const T_n& n, const T_N& N, const T_x& x, - const T_alpha& alpha, const T_beta& beta) { + const T_n& n, const T_N& N, const T_x& x, const T_alpha& alpha, + const T_beta& beta) { constexpr int T_x_rows = T_x::RowsAtCompileTime; using T_xbeta_partials = partials_return_t; using T_partials_return = partials_return_t; @@ -165,8 +165,8 @@ return_type_t binomial_logit_glm_lpmf( template inline return_type_t binomial_logit_glm_lpmf( - const T_n& n, const T_N& N, const T_x& x, - const T_alpha& alpha, const T_beta& beta) { + const T_n& n, const T_N& N, const T_x& x, const T_alpha& alpha, + const T_beta& beta) { return binomial_logit_glm_lpmf(n, N, x, alpha, beta); } } // namespace math diff --git a/test/unit/math/prim/prob/binomial_logit_glm_lpmf_test.cpp b/test/unit/math/prim/prob/binomial_logit_glm_lpmf_test.cpp index f23e78b890c..fb3298861eb 100644 --- a/test/unit/math/prim/prob/binomial_logit_glm_lpmf_test.cpp +++ b/test/unit/math/prim/prob/binomial_logit_glm_lpmf_test.cpp @@ -2,8 +2,8 @@ #include TEST(ProbBinomialLogitGLM, matchesNonGLM) { - using stan::math::binomial_logit_lpmf; using stan::math::binomial_logit_glm_lpmf; + using stan::math::binomial_logit_lpmf; std::vector n{1, 2}; std::vector N{5, 4}; @@ -58,13 +58,13 @@ TEST(ProbBinomialLogitGLM, throwsCorrectly) { std::vector N_mismatch_size{5, 4, 10}; EXPECT_THROW(binomial_logit_glm_lpmf(n, N_mismatch_size, x, alpha, beta), - std::invalid_argument); + std::invalid_argument); EXPECT_THROW(binomial_logit_glm_lpmf(500, 1, x, alpha, beta), - std::domain_error); + std::domain_error); EXPECT_THROW(binomial_logit_glm_lpmf(-10, N, x, alpha, beta), - std::domain_error); + std::domain_error); EXPECT_THROW(binomial_logit_glm_lpmf(n, -10, x, alpha, beta), - std::domain_error); + std::domain_error); Eigen::VectorXd alpha_inf = alpha; alpha[0] = INFTY; @@ -74,9 +74,9 @@ TEST(ProbBinomialLogitGLM, throwsCorrectly) { x(0, 0) = INFTY; EXPECT_THROW(binomial_logit_glm_lpmf(n, N, x_inf, alpha, beta), - std::domain_error); + std::domain_error); EXPECT_THROW(binomial_logit_glm_lpmf(n, N, x, alpha_inf, beta), - std::domain_error); + std::domain_error); EXPECT_THROW(binomial_logit_glm_lpmf(n, N, x, alpha, beta_inf), - std::domain_error); + std::domain_error); } From 96c757eb7cc829191279f0971366b22555ec2a42 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 20 Sep 2023 09:29:28 +0300 Subject: [PATCH 4/6] Doc update --- stan/math/prim/prob/binomial_logit_glm_lpmf.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/stan/math/prim/prob/binomial_logit_glm_lpmf.hpp b/stan/math/prim/prob/binomial_logit_glm_lpmf.hpp index 8d8d5a75308..8556a5067e8 100644 --- a/stan/math/prim/prob/binomial_logit_glm_lpmf.hpp +++ b/stan/math/prim/prob/binomial_logit_glm_lpmf.hpp @@ -26,9 +26,9 @@ namespace math { * simplified gradients. * If containers are supplied, returns the log sum of the probabilities. * - * @tparam T_y type of binary vector of successes variables; + * @tparam T_n type of binary vector of successes variables; * this can also be a single binary value; - * @tparam T_y type of binary vector of population size variables; + * @tparam T_N type of binary vector of population size variables; * this can also be a single binary value; * @tparam T_x type of the matrix of independent variables (features) * @tparam T_alpha type of the intercept(s); @@ -42,7 +42,7 @@ namespace math { * broadcast - used for all instances. * @param x design matrix or row vector. If it is a row vector it will be * broadcast - used for all instances. - * @param alpha intercept (in log odds) + * @param alpha intercept * @param beta weight vector * @return log probability or log sum of probabilities * @throw std::domain_error if x, beta or alpha is infinite. From d0bb70ee813501039f2e8642a7df63c58972f7c9 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 20 Sep 2023 09:42:48 +0300 Subject: [PATCH 5/6] Missed header --- stan/math/prim/prob/binomial_logit_glm_lpmf.hpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/stan/math/prim/prob/binomial_logit_glm_lpmf.hpp b/stan/math/prim/prob/binomial_logit_glm_lpmf.hpp index 8556a5067e8..a62cacc8f6f 100644 --- a/stan/math/prim/prob/binomial_logit_glm_lpmf.hpp +++ b/stan/math/prim/prob/binomial_logit_glm_lpmf.hpp @@ -13,6 +13,7 @@ #include #include #include +#include namespace stan { namespace math { @@ -58,9 +59,6 @@ return_type_t binomial_logit_glm_lpmf( constexpr int T_x_rows = T_x::RowsAtCompileTime; using T_xbeta_partials = partials_return_t; using T_partials_return = partials_return_t; - using T_theta_tmp = - typename std::conditional_t>; using T_xbeta_tmp = typename std::conditional_t>; From b89147268908262ad7d9c308fdd1385d436306e8 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Thu, 5 Oct 2023 15:03:09 +0300 Subject: [PATCH 6/6] Update doc, ref handlings --- .../prim/prob/binomial_logit_glm_lpmf.hpp | 21 +++++++++---------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/stan/math/prim/prob/binomial_logit_glm_lpmf.hpp b/stan/math/prim/prob/binomial_logit_glm_lpmf.hpp index a62cacc8f6f..afd038cd23c 100644 --- a/stan/math/prim/prob/binomial_logit_glm_lpmf.hpp +++ b/stan/math/prim/prob/binomial_logit_glm_lpmf.hpp @@ -18,7 +18,7 @@ namespace stan { namespace math { -/** \ingroup multivar_dists +/** \ingroup prob_dists * Returns the log PMF of the Generalized Linear Model (GLM) * with Binomial distribution and logit link function. * The idea is that binomial_logit_glm_lpmf(n | N, x, alpha, beta) should @@ -93,15 +93,15 @@ return_type_t binomial_logit_glm_lpmf( check_consistent_size(function, "Weight vector", beta, N_attributes); check_consistent_size(function, "Vector of intercepts", alpha, N_instances); - const auto& n_val = as_value_column_array_or_scalar(n_ref); - const auto& N_val = as_value_column_array_or_scalar(N_ref); + auto&& n_val = as_value_column_array_or_scalar(n_ref); + auto&& N_val = as_value_column_array_or_scalar(N_ref); check_bounded(function, "Successes variable", n_val, 0, N_val); check_nonnegative(function, "Population size parameter", N_val); - const auto& alpha_val = as_value_column_array_or_scalar(alpha_ref); - const auto& beta_val = as_value_column_vector_or_scalar(beta_ref); - const auto& x_val = value_of(x_ref); + auto&& alpha_val = as_value_column_array_or_scalar(alpha_ref); + auto&& beta_val = as_value_column_vector_or_scalar(beta_ref); + auto&& x_val = value_of(x_ref); Eigen::Array theta(N_instances); if (T_x_rows == 1) { theta = forward_as((x_val * beta_val)(0, 0)) + alpha_val; @@ -109,11 +109,11 @@ return_type_t binomial_logit_glm_lpmf( theta = (x_val * beta_val).array() + alpha_val; } - const auto& log_inv_logit_theta = log_inv_logit(theta); - const auto& log1m_inv_logit_theta = log1m_inv_logit(theta); + constexpr bool gradients_calc = !is_constant_all::value; + auto&& log_inv_logit_theta = to_ref_if(log_inv_logit(theta)); T_partials_return logp = sum(n_val * log_inv_logit_theta - + (N_val - n_val) * log1m_inv_logit_theta); + + (N_val - n_val) * log1m_inv_logit(theta)); using std::isfinite; if (!isfinite(logp)) { @@ -128,8 +128,7 @@ return_type_t binomial_logit_glm_lpmf( } auto ops_partials = make_partials_propagator(x_ref, alpha_ref, beta_ref); - - if (!is_constant_all::value) { + if (gradients_calc) { Eigen::Matrix theta_derivative = n_val - N_val * exp(log_inv_logit_theta);