From d2b5f26cb4cfdf01700b52638760208d14f662c3 Mon Sep 17 00:00:00 2001 From: Ling Zhi <1336265834@qq.com> Date: Tue, 24 Sep 2024 17:42:14 +0800 Subject: [PATCH 1/4] Feature/issue-3107 add beta negative binomial lpmf with test --- stan/math/prim/prob.hpp | 1 + .../math/prim/prob/beta_neg_binomial_lpmf.hpp | 220 ++++++++++++++++++ .../beta_neg_binomial_test.hpp | 98 ++++++++ 3 files changed, 319 insertions(+) create mode 100644 stan/math/prim/prob/beta_neg_binomial_lpmf.hpp create mode 100644 test/prob/beta_neg_binomial/beta_neg_binomial_test.hpp diff --git a/stan/math/prim/prob.hpp b/stan/math/prim/prob.hpp index f1b4efa5b6f..0760c087a1c 100644 --- a/stan/math/prim/prob.hpp +++ b/stan/math/prim/prob.hpp @@ -25,6 +25,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/prob/beta_neg_binomial_lpmf.hpp b/stan/math/prim/prob/beta_neg_binomial_lpmf.hpp new file mode 100644 index 00000000000..f55dc64d88b --- /dev/null +++ b/stan/math/prim/prob/beta_neg_binomial_lpmf.hpp @@ -0,0 +1,220 @@ +#ifndef STAN_MATH_PRIM_PROB_BETA_NEG_BINOMIAL_LPMF_HPP +#define STAN_MATH_PRIM_PROB_BETA_NEG_BINOMIAL_LPMF_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** \ingroup prob_dists + * Returns the log PMF of the Beta Negative Binomial distribution with given + * number of successes, prior success, and prior failure parameters. + * Given containers of matching sizes, returns the log sum of probabilities. + * + * @tparam T_n type of failure parameter + * @tparam T_r type of number of successes parameter + * @tparam T_alpha type of prior success parameter + * @tparam T_beta type of prior failure parameter + * + * @param n failure parameter + * @param r Number of successes parameter + * @param alpha prior success parameter + * @param beta prior failure parameter + * @return log probability or log sum of probabilities + * @throw std::domain_error if r, alpha, or beta fails to be positive + * @throw std::invalid_argument if container sizes mismatch + */ +template * = nullptr> +inline return_type_t beta_neg_binomial_lpmf( + const T_n& n, const T_r& r, const T_alpha& alpha, const T_beta& beta) { + using T_partials_return = partials_return_t; + using T_n_ref = ref_type_t; + using T_r_ref = ref_type_t; + using T_alpha_ref = ref_type_t; + using T_beta_ref = ref_type_t; + static constexpr const char* function = "beta_neg_binomial_lpmf"; + check_consistent_sizes( + function, "Failures variable", n, "Number of successes parameter", r, + "Prior success parameter", alpha, "Prior failure parameter", beta); + if (size_zero(n, r, alpha, beta)) { + return 0.0; + } + + T_n_ref n_ref = n; + T_r_ref r_ref = r; + T_alpha_ref alpha_ref = alpha; + T_beta_ref beta_ref = beta; + check_nonnegative(function, "Failures variable", n_ref); + check_positive_finite(function, "Number of successes parameter", r_ref); + check_positive_finite(function, "Prior success parameter", alpha_ref); + check_positive_finite(function, "Prior failure parameter", beta_ref); + + if (!include_summand::value) { + return 0.0; + } + + T_partials_return logp(0.0); + auto ops_partials = make_partials_propagator(r_ref, alpha_ref, beta_ref); + + scalar_seq_view n_vec(n); + scalar_seq_view r_vec(r_ref); + scalar_seq_view alpha_vec(alpha_ref); + scalar_seq_view beta_vec(beta_ref); + size_t size_n = stan::math::size(n); + size_t size_r = stan::math::size(r); + size_t size_alpha = stan::math::size(alpha); + size_t size_beta = stan::math::size(beta); + size_t size_n_r = max_size(n, r); + size_t size_r_alpha = max_size(r, alpha); + size_t size_n_beta = max_size(n, beta); + size_t size_alpha_beta = max_size(alpha, beta); + size_t max_size_seq_view = max_size(n, r, alpha, beta); + + VectorBuilder::value, T_partials_return, T_n> + normalizing_constant(size_n); + for (size_t i = 0; i < size_n; i++) + if (include_summand::value) + normalizing_constant[i] = -lgamma(n_vec[i] + 1); + + VectorBuilder lbeta_denominator( + size_r_alpha); + for (size_t i = 0; i < size_r_alpha; i++) { + lbeta_denominator[i] = lbeta(r_vec.val(i), alpha_vec.val(i)); + } + + VectorBuilder lgamma_denominator(size_beta); + for (size_t i = 0; i < size_beta; i++) { + lgamma_denominator[i] = lgamma(beta_vec.val(i)); + } + + VectorBuilder lgamma_numerator( + size_n_beta); + for (size_t i = 0; i < size_n_beta; i++) { + lgamma_numerator[i] = lgamma(n_vec[i] + beta_vec.val(i)); + } + + VectorBuilder lbeta_diff( + max_size_seq_view); + for (size_t i = 0; i < max_size_seq_view; i++) { + lbeta_diff[i] + = lbeta(n_vec[i] + r_vec.val(i), alpha_vec.val(i) + beta_vec.val(i)) + + lgamma_numerator[i] - lbeta_denominator[i] - lgamma_denominator[i]; + } + + // derivative w.r.t. r, alpha and beta + + VectorBuilder::value, + T_partials_return, T_n, T_r, T_alpha, T_beta> + digamma_n_r_alpha_beta(max_size_seq_view); + if (!is_constant_all::value) { + for (size_t i = 0; i < max_size_seq_view; i++) { + digamma_n_r_alpha_beta[i] = digamma(n_vec[i] + r_vec.val(i) + + alpha_vec.val(i) + beta_vec.val(i)); + } + } + + VectorBuilder::value, T_partials_return, + T_alpha, T_beta> + digamma_alpha_beta(size_alpha_beta); + if (!is_constant_all::value) { + for (size_t i = 0; i < size_alpha_beta; i++) { + digamma_alpha_beta[i] = digamma(alpha_vec.val(i) + beta_vec.val(i)); + } + } + + VectorBuilder::value, T_partials_return, T_n, T_r> + digamma_n_r(size_n_r); + if (!is_constant_all::value) { + for (size_t i = 0; i < size_n_r; i++) { + digamma_n_r[i] = digamma(n_vec[i] + r_vec.val(i)); + } + } + + VectorBuilder::value, T_partials_return, T_r, + T_alpha> + digamma_r_alpha(size_r_alpha); + if (!is_constant_all::value) { + for (size_t i = 0; i < size_r_alpha; i++) { + digamma_r_alpha[i] = digamma(r_vec.val(i) + alpha_vec.val(i)); + } + } + + VectorBuilder::value, T_partials_return, T_n, + T_beta> + digamma_n_beta(size_n_beta); + if (!is_constant_all::value) { + for (size_t i = 0; i < size_n_beta; i++) { + digamma_n_beta[i] = digamma(n_vec[i] + beta_vec.val(i)); + } + } + + VectorBuilder::value, T_partials_return, T_r> digamma_r( + size_r); + if (!is_constant_all::value) { + for (size_t i = 0; i < size_r; i++) { + digamma_r[i] = digamma(r_vec.val(i)); + } + } + + VectorBuilder::value, T_partials_return, T_alpha> + digamma_alpha(size_alpha); + if (!is_constant_all::value) { + for (size_t i = 0; i < size_alpha; i++) { + digamma_alpha[i] = digamma(alpha_vec.val(i)); + } + } + + VectorBuilder::value, T_partials_return, T_beta> + digamma_beta(size_beta); + if (!is_constant_all::value) { + for (size_t i = 0; i < size_beta; i++) { + digamma_beta[i] = digamma(beta_vec.val(i)); + } + } + + for (size_t i = 0; i < max_size_seq_view; i++) { + if (include_summand::value) + logp += normalizing_constant[i]; + logp += lbeta_diff[i]; + + if (!is_constant_all::value) { + partials<0>(ops_partials)[i] += digamma_n_r[i] - digamma_n_r_alpha_beta[i] + - (digamma_r[i] - digamma_r_alpha[i]); + } + if (!is_constant_all::value) { + partials<1>(ops_partials)[i] += digamma_alpha_beta[i] + - digamma_n_r_alpha_beta[i] + - (digamma_alpha[i] - digamma_r_alpha[i]); + } + if (!is_constant_all::value) { + partials<2>(ops_partials)[i] += digamma_alpha_beta[i] + - digamma_n_r_alpha_beta[i] + + digamma_n_beta[i] - digamma_beta[i]; + } + } + return ops_partials.build(logp); +} + +template +inline return_type_t beta_neg_binomial_lpmf( + const T_n& n, const T_r& r, const T_alpha& alpha, const T_beta& beta) { + return beta_neg_binomial_lpmf(n, r, alpha, beta); +} + +} // namespace math +} // namespace stan +#endif diff --git a/test/prob/beta_neg_binomial/beta_neg_binomial_test.hpp b/test/prob/beta_neg_binomial/beta_neg_binomial_test.hpp new file mode 100644 index 00000000000..5df5c58cb49 --- /dev/null +++ b/test/prob/beta_neg_binomial/beta_neg_binomial_test.hpp @@ -0,0 +1,98 @@ +// Arguments: Ints, Doubles, Doubles, Doubles +#include +#include +#include + +using stan::math::var; +using std::numeric_limits; +using std::vector; + +class AgradDistributionsBetaNegBinomial : public AgradDistributionTest { + public: + void valid_values(vector >& parameters, + vector& log_prob) { + vector param(4); + + param[0] = 5; // n + param[1] = 20.0; // r + param[2] = 10.0; // alpha + param[3] = 25.0; // beta + parameters.push_back(param); + log_prob.push_back(-10.3681267949788); // expected log_prob + + param[0] = 10; // n + param[1] = 5.5; // r + param[2] = 2.5; // alpha + param[3] = 0.5; // beta + parameters.push_back(param); + log_prob.push_back(-5.166741878823932); // expected log_prob + } + + void invalid_values(vector& index, vector& value) { + // n + index.push_back(0U); + value.push_back(-1); + + // r + index.push_back(1U); + value.push_back(0.0); + + index.push_back(1U); + value.push_back(-1.0); + + index.push_back(1U); + value.push_back(std::numeric_limits::infinity()); + + // alpha + index.push_back(2U); + value.push_back(0.0); + + index.push_back(2U); + value.push_back(-1.0); + + index.push_back(2U); + value.push_back(std::numeric_limits::infinity()); + + // beta + index.push_back(3U); + value.push_back(0.0); + + index.push_back(3U); + value.push_back(-1.0); + + index.push_back(3U); + value.push_back(std::numeric_limits::infinity()); + } + + template + stan::return_type_t log_prob(const T_n& n, + const T_r& r, + const T_size1& alpha, + const T_size2& beta, + const T4&, const T5&) { + return stan::math::beta_neg_binomial_lpmf(n, r, alpha, beta); + } + + template + stan::return_type_t log_prob(const T_n& n, + const T_r& r, + const T_size1& alpha, + const T_size2& beta, + const T4&, const T5&) { + return stan::math::beta_neg_binomial_lpmf(n, r, alpha, beta); + } + + template + stan::return_type_t log_prob_function( + const T_n& n, const T_r& r, const T_size1& alpha, const T_size2& beta, + const T4&, const T5&) { + using stan::math::lbeta; + using stan::math::lgamma; + + return lbeta(n + r, alpha + beta) - lbeta(r, alpha) + lgamma(n + beta) + - lgamma(n + 1) - lgamma(beta); + } +}; From 11ba2890de0a264b64a7865ca7562e3ae93aad21 Mon Sep 17 00:00:00 2001 From: Ling Zhi <1336265834@qq.com> Date: Thu, 26 Sep 2024 16:31:09 +0800 Subject: [PATCH 2/4] move loops into one --- .../math/prim/prob/beta_neg_binomial_lpmf.hpp | 142 ++++-------------- 1 file changed, 33 insertions(+), 109 deletions(-) diff --git a/stan/math/prim/prob/beta_neg_binomial_lpmf.hpp b/stan/math/prim/prob/beta_neg_binomial_lpmf.hpp index f55dc64d88b..01f550cf203 100644 --- a/stan/math/prim/prob/beta_neg_binomial_lpmf.hpp +++ b/stan/math/prim/prob/beta_neg_binomial_lpmf.hpp @@ -84,126 +84,50 @@ inline return_type_t beta_neg_binomial_lpmf( size_t size_alpha_beta = max_size(alpha, beta); size_t max_size_seq_view = max_size(n, r, alpha, beta); - VectorBuilder::value, T_partials_return, T_n> - normalizing_constant(size_n); - for (size_t i = 0; i < size_n; i++) - if (include_summand::value) - normalizing_constant[i] = -lgamma(n_vec[i] + 1); - - VectorBuilder lbeta_denominator( - size_r_alpha); - for (size_t i = 0; i < size_r_alpha; i++) { - lbeta_denominator[i] = lbeta(r_vec.val(i), alpha_vec.val(i)); - } - - VectorBuilder lgamma_denominator(size_beta); - for (size_t i = 0; i < size_beta; i++) { - lgamma_denominator[i] = lgamma(beta_vec.val(i)); - } - - VectorBuilder lgamma_numerator( - size_n_beta); - for (size_t i = 0; i < size_n_beta; i++) { - lgamma_numerator[i] = lgamma(n_vec[i] + beta_vec.val(i)); - } - - VectorBuilder lbeta_diff( - max_size_seq_view); for (size_t i = 0; i < max_size_seq_view; i++) { - lbeta_diff[i] - = lbeta(n_vec[i] + r_vec.val(i), alpha_vec.val(i) + beta_vec.val(i)) - + lgamma_numerator[i] - lbeta_denominator[i] - lgamma_denominator[i]; - } - - // derivative w.r.t. r, alpha and beta - - VectorBuilder::value, - T_partials_return, T_n, T_r, T_alpha, T_beta> - digamma_n_r_alpha_beta(max_size_seq_view); - if (!is_constant_all::value) { - for (size_t i = 0; i < max_size_seq_view; i++) { - digamma_n_r_alpha_beta[i] = digamma(n_vec[i] + r_vec.val(i) - + alpha_vec.val(i) + beta_vec.val(i)); + const T_partials_return lbeta_denominator + = lbeta(r_vec.val(i), alpha_vec.val(i)); + const T_partials_return lgamma_numerator + = lgamma(n_vec[i] + beta_vec.val(i)); + const T_partials_return lgamma_denominator = lgamma(beta_vec.val(i)); + const T_partials_return lbeta_numerator + = lbeta(n_vec[i] + r_vec.val(i), alpha_vec.val(i) + beta_vec.val(i)); + if (include_summand::value) { + logp -= lgamma(n_vec[i] + 1); } - } + logp += lbeta_numerator + lgamma_numerator - lbeta_denominator + - lgamma_denominator; - VectorBuilder::value, T_partials_return, - T_alpha, T_beta> - digamma_alpha_beta(size_alpha_beta); - if (!is_constant_all::value) { - for (size_t i = 0; i < size_alpha_beta; i++) { - digamma_alpha_beta[i] = digamma(alpha_vec.val(i) + beta_vec.val(i)); - } - } + T_partials_return digamma_n_r_alpha_beta + = is_constant_all::value + ? 0 + : digamma(n_vec[i] + r_vec.val(i) + alpha_vec.val(i) + + beta_vec.val(i)); - VectorBuilder::value, T_partials_return, T_n, T_r> - digamma_n_r(size_n_r); - if (!is_constant_all::value) { - for (size_t i = 0; i < size_n_r; i++) { - digamma_n_r[i] = digamma(n_vec[i] + r_vec.val(i)); - } - } + T_partials_return digamma_r_alpha + = is_constant_all::value + ? 0 + : digamma(r_vec.val(i) + alpha_vec.val(i)); - VectorBuilder::value, T_partials_return, T_r, - T_alpha> - digamma_r_alpha(size_r_alpha); - if (!is_constant_all::value) { - for (size_t i = 0; i < size_r_alpha; i++) { - digamma_r_alpha[i] = digamma(r_vec.val(i) + alpha_vec.val(i)); - } - } - - VectorBuilder::value, T_partials_return, T_n, - T_beta> - digamma_n_beta(size_n_beta); - if (!is_constant_all::value) { - for (size_t i = 0; i < size_n_beta; i++) { - digamma_n_beta[i] = digamma(n_vec[i] + beta_vec.val(i)); - } - } - - VectorBuilder::value, T_partials_return, T_r> digamma_r( - size_r); - if (!is_constant_all::value) { - for (size_t i = 0; i < size_r; i++) { - digamma_r[i] = digamma(r_vec.val(i)); - } - } - - VectorBuilder::value, T_partials_return, T_alpha> - digamma_alpha(size_alpha); - if (!is_constant_all::value) { - for (size_t i = 0; i < size_alpha; i++) { - digamma_alpha[i] = digamma(alpha_vec.val(i)); - } - } - - VectorBuilder::value, T_partials_return, T_beta> - digamma_beta(size_beta); - if (!is_constant_all::value) { - for (size_t i = 0; i < size_beta; i++) { - digamma_beta[i] = digamma(beta_vec.val(i)); - } - } - - for (size_t i = 0; i < max_size_seq_view; i++) { - if (include_summand::value) - logp += normalizing_constant[i]; - logp += lbeta_diff[i]; + T_partials_return digamma_alpha_beta + = is_constant_all::value + ? 0 + : digamma(alpha_vec.val(i) + beta_vec.val(i)); if (!is_constant_all::value) { - partials<0>(ops_partials)[i] += digamma_n_r[i] - digamma_n_r_alpha_beta[i] - - (digamma_r[i] - digamma_r_alpha[i]); + partials<0>(ops_partials)[i] + += digamma(n_vec[i] + r_vec.val(i)) - digamma_n_r_alpha_beta + - (digamma(r_vec.val(i)) - digamma_r_alpha); } if (!is_constant_all::value) { - partials<1>(ops_partials)[i] += digamma_alpha_beta[i] - - digamma_n_r_alpha_beta[i] - - (digamma_alpha[i] - digamma_r_alpha[i]); + partials<1>(ops_partials)[i] + += digamma_alpha_beta - digamma_n_r_alpha_beta + - (digamma(alpha_vec.val(i)) - digamma_r_alpha); } if (!is_constant_all::value) { - partials<2>(ops_partials)[i] += digamma_alpha_beta[i] - - digamma_n_r_alpha_beta[i] - + digamma_n_beta[i] - digamma_beta[i]; + partials<2>(ops_partials)[i] + += digamma_alpha_beta - digamma_n_r_alpha_beta + + digamma(n_vec[i] + beta_vec.val(i)) - digamma(beta_vec.val(i)); } } return ops_partials.build(logp); From 0731eb3dbbe36b5011563fccfbe1a8aa7231424d Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Mon, 30 Sep 2024 15:50:32 -0400 Subject: [PATCH 3/4] uses constexpr for whether to do extra calculations --- .../math/prim/prob/beta_neg_binomial_lpmf.hpp | 87 +++++++++---------- 1 file changed, 39 insertions(+), 48 deletions(-) diff --git a/stan/math/prim/prob/beta_neg_binomial_lpmf.hpp b/stan/math/prim/prob/beta_neg_binomial_lpmf.hpp index 01f550cf203..942d95a528b 100644 --- a/stan/math/prim/prob/beta_neg_binomial_lpmf.hpp +++ b/stan/math/prim/prob/beta_neg_binomial_lpmf.hpp @@ -63,71 +63,62 @@ inline return_type_t beta_neg_binomial_lpmf( check_positive_finite(function, "Prior success parameter", alpha_ref); check_positive_finite(function, "Prior failure parameter", beta_ref); - if (!include_summand::value) { + if constexpr (!include_summand::value) { return 0.0; } - T_partials_return logp(0.0); auto ops_partials = make_partials_propagator(r_ref, alpha_ref, beta_ref); scalar_seq_view n_vec(n); scalar_seq_view r_vec(r_ref); scalar_seq_view alpha_vec(alpha_ref); scalar_seq_view beta_vec(beta_ref); - size_t size_n = stan::math::size(n); - size_t size_r = stan::math::size(r); - size_t size_alpha = stan::math::size(alpha); - size_t size_beta = stan::math::size(beta); - size_t size_n_r = max_size(n, r); - size_t size_r_alpha = max_size(r, alpha); - size_t size_n_beta = max_size(n, beta); - size_t size_alpha_beta = max_size(alpha, beta); - size_t max_size_seq_view = max_size(n, r, alpha, beta); - + const size_t max_size_seq_view = max_size(n, r, alpha, beta); + T_partials_return logp(0.0); for (size_t i = 0; i < max_size_seq_view; i++) { - const T_partials_return lbeta_denominator + if constexpr (include_summand::value) { + logp -= lgamma(n_vec[i] + 1); + } + T_partials_return lbeta_denominator = lbeta(r_vec.val(i), alpha_vec.val(i)); - const T_partials_return lgamma_numerator + T_partials_return lgamma_numerator = lgamma(n_vec[i] + beta_vec.val(i)); - const T_partials_return lgamma_denominator = lgamma(beta_vec.val(i)); - const T_partials_return lbeta_numerator + T_partials_return lgamma_denominator = lgamma(beta_vec.val(i)); + T_partials_return lbeta_numerator = lbeta(n_vec[i] + r_vec.val(i), alpha_vec.val(i) + beta_vec.val(i)); - if (include_summand::value) { - logp -= lgamma(n_vec[i] + 1); - } logp += lbeta_numerator + lgamma_numerator - lbeta_denominator - lgamma_denominator; + if (!is_constant_all::value) { + T_partials_return digamma_n_r_alpha_beta + = digamma(n_vec[i] + r_vec.val(i) + alpha_vec.val(i) + + beta_vec.val(i)); - T_partials_return digamma_n_r_alpha_beta - = is_constant_all::value - ? 0 - : digamma(n_vec[i] + r_vec.val(i) + alpha_vec.val(i) - + beta_vec.val(i)); - - T_partials_return digamma_r_alpha - = is_constant_all::value - ? 0 - : digamma(r_vec.val(i) + alpha_vec.val(i)); - - T_partials_return digamma_alpha_beta - = is_constant_all::value - ? 0 - : digamma(alpha_vec.val(i) + beta_vec.val(i)); + if constexpr (!is_constant::value || !is_constant::value) { + T_partials_return digamma_r_alpha = digamma(r_vec.val(i) + alpha_vec.val(i)); + if constexpr (!is_constant_all::value) { + partials<0>(ops_partials)[i] + += digamma(n_vec[i] + r_vec.val(i)) - digamma_n_r_alpha_beta + - (digamma(r_vec.val(i)) - digamma_r_alpha); + } + if constexpr (!is_constant_all::value) { + partials<1>(ops_partials)[i] + += -digamma_n_r_alpha_beta + - (digamma(alpha_vec.val(i)) - digamma_r_alpha); + } + } + if constexpr (!is_constant::value || !is_constant::value) { + T_partials_return digamma_alpha_beta + = digamma(alpha_vec.val(i) + beta_vec.val(i)); + if constexpr (!is_constant_all::value) { + partials<2>(ops_partials)[i] + += digamma_alpha_beta - digamma_n_r_alpha_beta + + digamma(n_vec[i] + beta_vec.val(i)) - digamma(beta_vec.val(i)); + } + if constexpr (!is_constant_all::value) { + partials<1>(ops_partials)[i] += digamma_alpha_beta; + } + } - if (!is_constant_all::value) { - partials<0>(ops_partials)[i] - += digamma(n_vec[i] + r_vec.val(i)) - digamma_n_r_alpha_beta - - (digamma(r_vec.val(i)) - digamma_r_alpha); - } - if (!is_constant_all::value) { - partials<1>(ops_partials)[i] - += digamma_alpha_beta - digamma_n_r_alpha_beta - - (digamma(alpha_vec.val(i)) - digamma_r_alpha); - } - if (!is_constant_all::value) { - partials<2>(ops_partials)[i] - += digamma_alpha_beta - digamma_n_r_alpha_beta - + digamma(n_vec[i] + beta_vec.val(i)) - digamma(beta_vec.val(i)); } } return ops_partials.build(logp); From 25650dce7ed2d27ac57c9ffd6ac6fdd59a94802d Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 30 Sep 2024 15:54:14 -0400 Subject: [PATCH 4/4] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- .../math/prim/prob/beta_neg_binomial_lpmf.hpp | 29 +++++++++---------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/stan/math/prim/prob/beta_neg_binomial_lpmf.hpp b/stan/math/prim/prob/beta_neg_binomial_lpmf.hpp index 942d95a528b..55ae09828d4 100644 --- a/stan/math/prim/prob/beta_neg_binomial_lpmf.hpp +++ b/stan/math/prim/prob/beta_neg_binomial_lpmf.hpp @@ -79,46 +79,45 @@ inline return_type_t beta_neg_binomial_lpmf( if constexpr (include_summand::value) { logp -= lgamma(n_vec[i] + 1); } - T_partials_return lbeta_denominator - = lbeta(r_vec.val(i), alpha_vec.val(i)); - T_partials_return lgamma_numerator - = lgamma(n_vec[i] + beta_vec.val(i)); + T_partials_return lbeta_denominator = lbeta(r_vec.val(i), alpha_vec.val(i)); + T_partials_return lgamma_numerator = lgamma(n_vec[i] + beta_vec.val(i)); T_partials_return lgamma_denominator = lgamma(beta_vec.val(i)); T_partials_return lbeta_numerator = lbeta(n_vec[i] + r_vec.val(i), alpha_vec.val(i) + beta_vec.val(i)); logp += lbeta_numerator + lgamma_numerator - lbeta_denominator - lgamma_denominator; if (!is_constant_all::value) { - T_partials_return digamma_n_r_alpha_beta - = digamma(n_vec[i] + r_vec.val(i) + alpha_vec.val(i) - + beta_vec.val(i)); + T_partials_return digamma_n_r_alpha_beta = digamma( + n_vec[i] + r_vec.val(i) + alpha_vec.val(i) + beta_vec.val(i)); if constexpr (!is_constant::value || !is_constant::value) { - T_partials_return digamma_r_alpha = digamma(r_vec.val(i) + alpha_vec.val(i)); + T_partials_return digamma_r_alpha + = digamma(r_vec.val(i) + alpha_vec.val(i)); if constexpr (!is_constant_all::value) { partials<0>(ops_partials)[i] += digamma(n_vec[i] + r_vec.val(i)) - digamma_n_r_alpha_beta - - (digamma(r_vec.val(i)) - digamma_r_alpha); + - (digamma(r_vec.val(i)) - digamma_r_alpha); } if constexpr (!is_constant_all::value) { partials<1>(ops_partials)[i] += -digamma_n_r_alpha_beta - - (digamma(alpha_vec.val(i)) - digamma_r_alpha); + - (digamma(alpha_vec.val(i)) - digamma_r_alpha); } } - if constexpr (!is_constant::value || !is_constant::value) { + if constexpr (!is_constant::value + || !is_constant::value) { T_partials_return digamma_alpha_beta = digamma(alpha_vec.val(i) + beta_vec.val(i)); if constexpr (!is_constant_all::value) { - partials<2>(ops_partials)[i] - += digamma_alpha_beta - digamma_n_r_alpha_beta - + digamma(n_vec[i] + beta_vec.val(i)) - digamma(beta_vec.val(i)); + partials<2>(ops_partials)[i] += digamma_alpha_beta + - digamma_n_r_alpha_beta + + digamma(n_vec[i] + beta_vec.val(i)) + - digamma(beta_vec.val(i)); } if constexpr (!is_constant_all::value) { partials<1>(ops_partials)[i] += digamma_alpha_beta; } } - } } return ops_partials.build(logp);