From cfcff04b2d6f87ae0e8876ac55e3312aa2d47a4a Mon Sep 17 00:00:00 2001 From: Ben Date: Thu, 5 Mar 2020 18:30:22 -0500 Subject: [PATCH 01/16] Changed logic for select_var_param select the nth parameter instead of always picking the 0th parameter (Issue #1763) --- test/prob/utility.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/prob/utility.hpp b/test/prob/utility.hpp index 6d3a4fb42f2..ca4a136d7a3 100644 --- a/test/prob/utility.hpp +++ b/test/prob/utility.hpp @@ -596,7 +596,7 @@ typename scalar_type::type select_var_param( const vector>& parameters, const size_t n, const size_t p) { typename scalar_type::type param(0); if (p < parameters[0].size()) { - if (is_vector::value && !is_constant_all::value) + if (!is_vector::value && !is_constant_all::value) param = parameters[n][p]; else param = parameters[0][p]; From 4ffbd3c0c6fb5dcef4b557821f7511a9abcb1301 Mon Sep 17 00:00:00 2001 From: Ben Date: Sat, 2 May 2020 11:49:48 -0400 Subject: [PATCH 02/16] Check value_of_rec in results so that propto check works for all autodiff types (Issue #1763) --- test/prob/test_fixture_distr.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/prob/test_fixture_distr.hpp b/test/prob/test_fixture_distr.hpp index cfcd73e6b4e..d16c72565fa 100644 --- a/test/prob/test_fixture_distr.hpp +++ b/test/prob/test_fixture_distr.hpp @@ -12,6 +12,7 @@ using stan::is_constant_all; using stan::is_vector; using stan::math::fvar; using stan::math::value_of; +using stan::math::value_of_rec; using stan::math::var; using stan::scalar_type; using std::vector; @@ -251,8 +252,8 @@ class AgradDistributionTestFixture : public ::testing::Test { Scalar3, Scalar4, Scalar5>(p0, p1, p2, p3, p4, p5); - EXPECT_TRUE(reference_logprob_false - logprob_false - == reference_logprob_true - logprob_true) + EXPECT_TRUE(value_of_rec(reference_logprob_false - logprob_false) + == value_of_rec(reference_logprob_true - logprob_true)) << "Proportional test failed at index: " << n << std::endl << " reference params: " << parameters[0] << std::endl << " current params: " << parameters[n] << std::endl From 56988042310b7d28eaddbf2c4fde2799884a3339 Mon Sep 17 00:00:00 2001 From: Ben Date: Sat, 2 May 2020 12:51:39 -0400 Subject: [PATCH 03/16] Removed poisson approximation from negative binomial and updated the numerics a bit (Issue #1763) --- stan/math/prim/prob/neg_binomial_lpmf.hpp | 71 +++++++------------ test/prob/neg_binomial/neg_binomial_test.hpp | 6 +- .../math/prim/prob/neg_binomial_log_test.cpp | 3 + 3 files changed, 30 insertions(+), 50 deletions(-) diff --git a/stan/math/prim/prob/neg_binomial_lpmf.hpp b/stan/math/prim/prob/neg_binomial_lpmf.hpp index bef6448e730..e21e7646b29 100644 --- a/stan/math/prim/prob/neg_binomial_lpmf.hpp +++ b/stan/math/prim/prob/neg_binomial_lpmf.hpp @@ -64,69 +64,46 @@ return_type_t neg_binomial_lpmf(const T_n& n, } } - VectorBuilder log_beta(size_beta); + VectorBuilder log1p_inv_beta(size_beta); VectorBuilder log1p_beta(size_beta); - VectorBuilder log_beta_m_log1p_beta( - size_beta); for (size_t i = 0; i < size_beta; ++i) { const T_partials_return beta_dbl = value_of(beta_vec[i]); - log_beta[i] = log(beta_dbl); + log1p_inv_beta[i] = log1p(inv(beta_dbl)); log1p_beta[i] = log1p(beta_dbl); - log_beta_m_log1p_beta[i] = log_beta[i] - log1p_beta[i]; } - VectorBuilder lambda( - size_alpha_beta); - VectorBuilder - alpha_log_beta_over_1p_beta(size_alpha_beta); VectorBuilder::value, T_partials_return, T_shape, T_inv_scale> lambda_m_alpha_over_1p_beta(size_alpha_beta); - for (size_t i = 0; i < size_alpha_beta; ++i) { - const T_partials_return alpha_dbl = value_of(alpha_vec[i]); - const T_partials_return beta_dbl = value_of(beta_vec[i]); - lambda[i] = alpha_dbl / beta_dbl; - alpha_log_beta_over_1p_beta[i] = alpha_dbl * log_beta_m_log1p_beta[i]; - if (!is_constant_all::value) { - lambda_m_alpha_over_1p_beta[i] = lambda[i] - alpha_dbl / (1 + beta_dbl); + if (!is_constant_all::value) { + for (size_t i = 0; i < size_alpha_beta; ++i) { + const T_partials_return alpha_dbl = value_of(alpha_vec[i]); + const T_partials_return beta_dbl = value_of(beta_vec[i]); + lambda_m_alpha_over_1p_beta[i] = alpha_dbl / beta_dbl - alpha_dbl / (1 + beta_dbl); } } for (size_t i = 0; i < max_size_seq_view; i++) { - if (alpha_vec[i] > internal::neg_binomial_alpha_cutoff) { - // reduces numerically to Poisson - if (include_summand::value) { - logp -= lgamma(n_vec[i] + 1.0); - } - logp += multiply_log(n_vec[i], lambda[i]) - lambda[i]; + const T_partials_return alpha_dbl = value_of(alpha_vec[i]); + const T_partials_return beta_dbl = value_of(beta_vec[i]); - if (!is_constant_all::value) { - ops_partials.edge1_.partials_[i] - += n_vec[i] / value_of(alpha_vec[i]) - 1.0 / value_of(beta_vec[i]); - } - if (!is_constant_all::value) { - ops_partials.edge2_.partials_[i] - += (lambda[i] - n_vec[i]) / value_of(beta_vec[i]); + if (include_summand::value) { + if (n_vec[i] != 0) { + logp += binomial_coefficient_log( + n_vec[i] + alpha_dbl - 1.0, alpha_dbl - 1.0); } - } else { // standard density definition - if (include_summand::value) { - if (n_vec[i] != 0) { - logp += binomial_coefficient_log( - n_vec[i] + value_of(alpha_vec[i]) - 1.0, n_vec[i]); - } - } - logp += alpha_log_beta_over_1p_beta[i] - n_vec[i] * log1p_beta[i]; + } + logp += -alpha_dbl * log1p_inv_beta[i] - n_vec[i] * log1p_beta[i]; - if (!is_constant_all::value) { - ops_partials.edge1_.partials_[i] - += digamma(value_of(alpha_vec[i]) + n_vec[i]) - digamma_alpha[i] - + log_beta_m_log1p_beta[i]; - } - if (!is_constant_all::value) { - ops_partials.edge2_.partials_[i] - += lambda_m_alpha_over_1p_beta[i] - - n_vec[i] / (value_of(beta_vec[i]) + 1.0); - } + if (!is_constant_all::value) { + ops_partials.edge1_.partials_[i] + += digamma(alpha_dbl + n_vec[i]) - digamma_alpha[i] + - log1p_inv_beta[i]; + } + if (!is_constant_all::value) { + ops_partials.edge2_.partials_[i] + += lambda_m_alpha_over_1p_beta[i] + - n_vec[i] / (beta_dbl + 1.0); } } diff --git a/test/prob/neg_binomial/neg_binomial_test.hpp b/test/prob/neg_binomial/neg_binomial_test.hpp index 265ecaf129a..e7007131787 100644 --- a/test/prob/neg_binomial/neg_binomial_test.hpp +++ b/test/prob/neg_binomial/neg_binomial_test.hpp @@ -25,10 +25,10 @@ class AgradDistributionsNegBinomial : public AgradDistributionTest { log_prob.push_back(-142.6147368129045105434); // expected log_prob param[0] = 13; - param[1] = 1e11; // alpha > 1e10, causes redux to Poisson - param[2] = 1e10; // equiv to Poisson(1e11/1e10) = Poisson(10) + param[1] = 1e11; + param[2] = 1e10; parameters.push_back(param); - log_prob.push_back(-2.6185576442008289933); // log poisson(13|10) + log_prob.push_back(-2.6185576442208003); // expected log_prob } void invalid_values(vector& index, vector& value) { diff --git a/test/unit/math/prim/prob/neg_binomial_log_test.cpp b/test/unit/math/prim/prob/neg_binomial_log_test.cpp index 25754929b80..c5f6ec0f493 100644 --- a/test/unit/math/prim/prob/neg_binomial_log_test.cpp +++ b/test/unit/math/prim/prob/neg_binomial_log_test.cpp @@ -6,6 +6,9 @@ TEST(ProbNegBinomial, log_matches_lpmf) { double alpha = 1.1; double beta = 2.3; + EXPECT_FLOAT_EQ((stan::math::neg_binomial_lpmf(13, 1e11, 1e10)), + -2.6185576442208003); + EXPECT_FLOAT_EQ((stan::math::neg_binomial_lpmf(y, alpha, beta)), (stan::math::neg_binomial_log(y, alpha, beta))); EXPECT_FLOAT_EQ((stan::math::neg_binomial_lpmf(y, alpha, beta)), From fc3a47be7cc30018110df8bec4a99080e2e4e635 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sat, 2 May 2020 12:56:12 -0400 Subject: [PATCH 04/16] [Jenkins] auto-formatting by clang-format version 6.0.0 --- stan/math/prim/prob/neg_binomial_lpmf.hpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/stan/math/prim/prob/neg_binomial_lpmf.hpp b/stan/math/prim/prob/neg_binomial_lpmf.hpp index e21e7646b29..e5ddc0cb804 100644 --- a/stan/math/prim/prob/neg_binomial_lpmf.hpp +++ b/stan/math/prim/prob/neg_binomial_lpmf.hpp @@ -79,7 +79,8 @@ return_type_t neg_binomial_lpmf(const T_n& n, for (size_t i = 0; i < size_alpha_beta; ++i) { const T_partials_return alpha_dbl = value_of(alpha_vec[i]); const T_partials_return beta_dbl = value_of(beta_vec[i]); - lambda_m_alpha_over_1p_beta[i] = alpha_dbl / beta_dbl - alpha_dbl / (1 + beta_dbl); + lambda_m_alpha_over_1p_beta[i] + = alpha_dbl / beta_dbl - alpha_dbl / (1 + beta_dbl); } } @@ -89,21 +90,20 @@ return_type_t neg_binomial_lpmf(const T_n& n, if (include_summand::value) { if (n_vec[i] != 0) { - logp += binomial_coefficient_log( - n_vec[i] + alpha_dbl - 1.0, alpha_dbl - 1.0); + logp += binomial_coefficient_log(n_vec[i] + alpha_dbl - 1.0, + alpha_dbl - 1.0); } } logp += -alpha_dbl * log1p_inv_beta[i] - n_vec[i] * log1p_beta[i]; if (!is_constant_all::value) { - ops_partials.edge1_.partials_[i] - += digamma(alpha_dbl + n_vec[i]) - digamma_alpha[i] - - log1p_inv_beta[i]; + ops_partials.edge1_.partials_[i] += digamma(alpha_dbl + n_vec[i]) + - digamma_alpha[i] + - log1p_inv_beta[i]; } if (!is_constant_all::value) { ops_partials.edge2_.partials_[i] - += lambda_m_alpha_over_1p_beta[i] - - n_vec[i] / (beta_dbl + 1.0); + += lambda_m_alpha_over_1p_beta[i] - n_vec[i] / (beta_dbl + 1.0); } } From ff6046d1a6a22dbf0d7a26f60e9aa4c849718893 Mon Sep 17 00:00:00 2001 From: Ben Date: Fri, 17 Jul 2020 18:40:41 -0400 Subject: [PATCH 05/16] Removed vector check from select_var_param and added gradient checks to neg_binomial (Issue #1763) --- test/prob/utility.hpp | 2 +- test/unit/math/mix/prob/neg_binomial_test.cpp | 20 +++++++++++++++++++ 2 files changed, 21 insertions(+), 1 deletion(-) create mode 100644 test/unit/math/mix/prob/neg_binomial_test.cpp diff --git a/test/prob/utility.hpp b/test/prob/utility.hpp index 070da54cffa..53f9e966d0f 100644 --- a/test/prob/utility.hpp +++ b/test/prob/utility.hpp @@ -596,7 +596,7 @@ typename scalar_type::type select_var_param( const vector>& parameters, const size_t n, const size_t p) { typename scalar_type::type param(0); if (p < parameters[0].size()) { - if (!is_vector::value && !is_constant_all::value) + if (!is_constant_all::value) param = parameters[n][p]; else param = parameters[0][p]; diff --git a/test/unit/math/mix/prob/neg_binomial_test.cpp b/test/unit/math/mix/prob/neg_binomial_test.cpp new file mode 100644 index 00000000000..42df66a9343 --- /dev/null +++ b/test/unit/math/mix/prob/neg_binomial_test.cpp @@ -0,0 +1,20 @@ +#include +#include +#include +#include + +TEST(mathMixScalFun, neg_binomial_lpmf_derivatives) { + auto f1 = [](const auto& alpha, const auto& beta) { + return stan::math::neg_binomial_lpmf(0, alpha, beta); + }; + auto f2 = [](const auto& alpha, const auto& beta) { + return stan::math::neg_binomial_lpmf(6, alpha, beta); + }; + + stan::test::expect_ad(f1, 1.5, 4.1); + stan::test::expect_ad(f1, std::vector({ 1.2, 2.0 }), + std::vector({ 1.1, 1.2 })); + stan::test::expect_ad(f2, 1.5, 4.1); + stan::test::expect_ad(f2, std::vector({ 1.7, 2.0 }), + std::vector({ 1.1, 2.3 })); +} From 2ff7275b24f1f4cc41a04984cec81a529b322c3e Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Fri, 17 Jul 2020 22:42:09 +0000 Subject: [PATCH 06/16] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- test/unit/math/mix/prob/neg_binomial_test.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/unit/math/mix/prob/neg_binomial_test.cpp b/test/unit/math/mix/prob/neg_binomial_test.cpp index 42df66a9343..03326fcd25c 100644 --- a/test/unit/math/mix/prob/neg_binomial_test.cpp +++ b/test/unit/math/mix/prob/neg_binomial_test.cpp @@ -12,9 +12,9 @@ TEST(mathMixScalFun, neg_binomial_lpmf_derivatives) { }; stan::test::expect_ad(f1, 1.5, 4.1); - stan::test::expect_ad(f1, std::vector({ 1.2, 2.0 }), - std::vector({ 1.1, 1.2 })); + stan::test::expect_ad(f1, std::vector({1.2, 2.0}), + std::vector({1.1, 1.2})); stan::test::expect_ad(f2, 1.5, 4.1); - stan::test::expect_ad(f2, std::vector({ 1.7, 2.0 }), - std::vector({ 1.1, 2.3 })); + stan::test::expect_ad(f2, std::vector({1.7, 2.0}), + std::vector({1.1, 2.3})); } From 49de095257818f96b2583afc48f8e099963f960f Mon Sep 17 00:00:00 2001 From: Ben Date: Fri, 17 Jul 2020 19:42:25 -0400 Subject: [PATCH 07/16] Use long double in test to avoid precision loss (Issue #1763) --- test/unit/math/prim/prob/neg_binomial_test.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/unit/math/prim/prob/neg_binomial_test.cpp b/test/unit/math/prim/prob/neg_binomial_test.cpp index 34b0fcd705a..ac5b0f611cb 100644 --- a/test/unit/math/prim/prob/neg_binomial_test.cpp +++ b/test/unit/math/prim/prob/neg_binomial_test.cpp @@ -73,12 +73,12 @@ TEST(ProbDistributionsNegBinomial, error_check) { void expected_bin_sizes(double* expect, const int K, const int N, const double alpha, const double beta) { - double p = 0; + long double p = 0; for (int i = 0; i < K; i++) { expect[i] = N * std::exp(stan::math::neg_binomial_log(i, alpha, beta)); p += std::exp(stan::math::neg_binomial_log(i, alpha, beta)); } - expect[K - 1] = N * (1.0 - p); + expect[K - 1] = N * static_cast(1.0 - p); } TEST(ProbDistributionsNegBinomial, chiSquareGoodnessFitTest) { From 483a3c0945689cf183ec3a198e891d5e972b12cd Mon Sep 17 00:00:00 2001 From: Ben Date: Sun, 19 Jul 2020 15:36:19 -0400 Subject: [PATCH 08/16] Replaced EXPECT_EQ with EXPECT_NEAR for float comparison of normalizing constants (Issue #1763) --- test/prob/test_fixture_distr.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/prob/test_fixture_distr.hpp b/test/prob/test_fixture_distr.hpp index 90dca29a539..77f1d89a704 100644 --- a/test/prob/test_fixture_distr.hpp +++ b/test/prob/test_fixture_distr.hpp @@ -252,8 +252,8 @@ class AgradDistributionTestFixture : public ::testing::Test { Scalar3, Scalar4, Scalar5>(p0, p1, p2, p3, p4, p5); - EXPECT_TRUE(value_of_rec(reference_logprob_false - logprob_false) - == value_of_rec(reference_logprob_true - logprob_true)) + EXPECT_NEAR(value_of_rec(reference_logprob_false - logprob_false), + value_of_rec(reference_logprob_true - logprob_true), 1e-12) << "Proportional test failed at index: " << n << std::endl << " reference params: " << parameters[0] << std::endl << " current params: " << parameters[n] << std::endl From ed2b7f2e9ab3c36c5c0a3be5d7f0bdfa284a93bf Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 19 Jul 2020 19:42:47 +0000 Subject: [PATCH 09/16] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- test/prob/test_fixture_distr.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/prob/test_fixture_distr.hpp b/test/prob/test_fixture_distr.hpp index 77f1d89a704..e5fa3cae6e1 100644 --- a/test/prob/test_fixture_distr.hpp +++ b/test/prob/test_fixture_distr.hpp @@ -253,7 +253,7 @@ class AgradDistributionTestFixture : public ::testing::Test { p3, p4, p5); EXPECT_NEAR(value_of_rec(reference_logprob_false - logprob_false), - value_of_rec(reference_logprob_true - logprob_true), 1e-12) + value_of_rec(reference_logprob_true - logprob_true), 1e-12) << "Proportional test failed at index: " << n << std::endl << " reference params: " << parameters[0] << std::endl << " current params: " << parameters[n] << std::endl From 8213edf7be66e191d4d1021f142fccd7776eaf43 Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 20 Jul 2020 11:53:39 -0400 Subject: [PATCH 10/16] Added missing T_loc arguments to include_summand terms in pareto_type_2 (Issue #1763) --- stan/math/prim/prob/pareto_type_2_lpdf.hpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/stan/math/prim/prob/pareto_type_2_lpdf.hpp b/stan/math/prim/prob/pareto_type_2_lpdf.hpp index d3472bae72c..486321b4fd8 100644 --- a/stan/math/prim/prob/pareto_type_2_lpdf.hpp +++ b/stan/math/prim/prob/pareto_type_2_lpdf.hpp @@ -67,9 +67,9 @@ return_type_t pareto_type_2_lpdf( } } - VectorBuilder::value, T_partials_return, T_shape> + VectorBuilder::value, T_partials_return, T_shape> inv_alpha(size(alpha)); - if (!is_constant_all::value) { + if (include_summand::value) { for (size_t n = 0; n < stan::math::size(alpha); n++) { inv_alpha[n] = 1 / value_of(alpha_vec[n]); } @@ -85,9 +85,9 @@ return_type_t pareto_type_2_lpdf( const T_partials_return deriv_1_2 = inv_sum + alpha_div_sum; const T_partials_return log1p_scaled_diff - = include_summand::value - ? log1p((y_dbl - mu_dbl) / lambda_dbl) - : 0; + = include_summand::value + ? log1p((y_dbl - mu_dbl) / lambda_dbl) + : 0; if (include_summand::value) { logp += log_alpha[n]; @@ -95,7 +95,8 @@ return_type_t pareto_type_2_lpdf( if (include_summand::value) { logp -= log_lambda[n]; } - if (include_summand::value) { + + if (include_summand::value) { logp -= (alpha_dbl + 1.0) * log1p_scaled_diff; } From ea096c09f4d4668306bd9a1eb81a66129893357a Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 20 Jul 2020 15:58:59 +0000 Subject: [PATCH 11/16] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/prob/pareto_type_2_lpdf.hpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/stan/math/prim/prob/pareto_type_2_lpdf.hpp b/stan/math/prim/prob/pareto_type_2_lpdf.hpp index 486321b4fd8..7a08527ec14 100644 --- a/stan/math/prim/prob/pareto_type_2_lpdf.hpp +++ b/stan/math/prim/prob/pareto_type_2_lpdf.hpp @@ -67,7 +67,8 @@ return_type_t pareto_type_2_lpdf( } } - VectorBuilder::value, T_partials_return, T_shape> + VectorBuilder::value, T_partials_return, + T_shape> inv_alpha(size(alpha)); if (include_summand::value) { for (size_t n = 0; n < stan::math::size(alpha); n++) { @@ -85,9 +86,9 @@ return_type_t pareto_type_2_lpdf( const T_partials_return deriv_1_2 = inv_sum + alpha_div_sum; const T_partials_return log1p_scaled_diff - = include_summand::value - ? log1p((y_dbl - mu_dbl) / lambda_dbl) - : 0; + = include_summand::value + ? log1p((y_dbl - mu_dbl) / lambda_dbl) + : 0; if (include_summand::value) { logp += log_alpha[n]; @@ -95,7 +96,7 @@ return_type_t pareto_type_2_lpdf( if (include_summand::value) { logp -= log_lambda[n]; } - + if (include_summand::value) { logp -= (alpha_dbl + 1.0) * log1p_scaled_diff; } From ee2cf8dcf72a269cf19f2f0c8e1de15293d6ea88 Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 20 Jul 2020 12:22:32 -0400 Subject: [PATCH 12/16] Added n = 0, n = 1 tests to negative binomial pdf (Issue #1763) --- test/unit/math/prim/prob/neg_binomial_log_test.cpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/test/unit/math/prim/prob/neg_binomial_log_test.cpp b/test/unit/math/prim/prob/neg_binomial_log_test.cpp index c5f6ec0f493..6286d069e46 100644 --- a/test/unit/math/prim/prob/neg_binomial_log_test.cpp +++ b/test/unit/math/prim/prob/neg_binomial_log_test.cpp @@ -6,6 +6,13 @@ TEST(ProbNegBinomial, log_matches_lpmf) { double alpha = 1.1; double beta = 2.3; + long double alpha2 = 1e11; + long double beta2 = 1e10; + + EXPECT_FLOAT_EQ(stan::math::neg_binomial_lpmf(0, 1e11, 1e10), + -alpha2 * std::log1p(1.0 / beta2)); + EXPECT_FLOAT_EQ(stan::math::neg_binomial_lpmf(1, 1e11, 1e10), + std::log(alpha2 - 1.0) - alpha2 * std::log1p(1.0 / beta2) - std::log1p(beta2)); EXPECT_FLOAT_EQ((stan::math::neg_binomial_lpmf(13, 1e11, 1e10)), -2.6185576442208003); From 283a2b030f51031dc0f6b799eb20055acaac7283 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 20 Jul 2020 16:23:50 +0000 Subject: [PATCH 13/16] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- test/unit/math/prim/prob/neg_binomial_log_test.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/unit/math/prim/prob/neg_binomial_log_test.cpp b/test/unit/math/prim/prob/neg_binomial_log_test.cpp index 6286d069e46..5b5ae00c634 100644 --- a/test/unit/math/prim/prob/neg_binomial_log_test.cpp +++ b/test/unit/math/prim/prob/neg_binomial_log_test.cpp @@ -10,9 +10,10 @@ TEST(ProbNegBinomial, log_matches_lpmf) { long double beta2 = 1e10; EXPECT_FLOAT_EQ(stan::math::neg_binomial_lpmf(0, 1e11, 1e10), - -alpha2 * std::log1p(1.0 / beta2)); + -alpha2 * std::log1p(1.0 / beta2)); EXPECT_FLOAT_EQ(stan::math::neg_binomial_lpmf(1, 1e11, 1e10), - std::log(alpha2 - 1.0) - alpha2 * std::log1p(1.0 / beta2) - std::log1p(beta2)); + std::log(alpha2 - 1.0) - alpha2 * std::log1p(1.0 / beta2) + - std::log1p(beta2)); EXPECT_FLOAT_EQ((stan::math::neg_binomial_lpmf(13, 1e11, 1e10)), -2.6185576442208003); From 4e4ad0633ef75a5f059fe4d1df7c1cbf52855f33 Mon Sep 17 00:00:00 2001 From: Ben Date: Tue, 21 Jul 2020 11:25:28 -0400 Subject: [PATCH 14/16] Pulled alpha/beta values into their own containers (Issue #1763) --- stan/math/prim/prob/neg_binomial_lpmf.hpp | 27 ++++++++++------------- 1 file changed, 12 insertions(+), 15 deletions(-) diff --git a/stan/math/prim/prob/neg_binomial_lpmf.hpp b/stan/math/prim/prob/neg_binomial_lpmf.hpp index 7a2a0118b1a..df2682387f8 100644 --- a/stan/math/prim/prob/neg_binomial_lpmf.hpp +++ b/stan/math/prim/prob/neg_binomial_lpmf.hpp @@ -57,20 +57,22 @@ return_type_t neg_binomial_lpmf(const T_n& n, size_t size_alpha_beta = max_size(alpha, beta); size_t max_size_seq_view = max_size(n, alpha, beta); + scalar_seq_view alpha_dbl(value_of(alpha)); + scalar_seq_view beta_dbl(value_of(beta)); + VectorBuilder::value, T_partials_return, T_shape> digamma_alpha(size_alpha); if (!is_constant_all::value) { for (size_t i = 0; i < size_alpha; ++i) { - digamma_alpha[i] = digamma(value_of(alpha_vec[i])); + digamma_alpha[i] = digamma(alpha_dbl[i]); } } VectorBuilder log1p_inv_beta(size_beta); VectorBuilder log1p_beta(size_beta); for (size_t i = 0; i < size_beta; ++i) { - const T_partials_return beta_dbl = value_of(beta_vec[i]); - log1p_inv_beta[i] = log1p(inv(beta_dbl)); - log1p_beta[i] = log1p(beta_dbl); + log1p_inv_beta[i] = log1p(inv(beta_dbl[i])); + log1p_beta[i] = log1p(beta_dbl[i]); } VectorBuilder::value, T_partials_return, @@ -78,33 +80,28 @@ return_type_t neg_binomial_lpmf(const T_n& n, lambda_m_alpha_over_1p_beta(size_alpha_beta); if (!is_constant_all::value) { for (size_t i = 0; i < size_alpha_beta; ++i) { - const T_partials_return alpha_dbl = value_of(alpha_vec[i]); - const T_partials_return beta_dbl = value_of(beta_vec[i]); lambda_m_alpha_over_1p_beta[i] - = alpha_dbl / beta_dbl - alpha_dbl / (1 + beta_dbl); + = alpha_dbl[i] / beta_dbl[i] - alpha_dbl[i] / (1 + beta_dbl[i]); } } for (size_t i = 0; i < max_size_seq_view; i++) { - const T_partials_return alpha_dbl = value_of(alpha_vec[i]); - const T_partials_return beta_dbl = value_of(beta_vec[i]); - if (include_summand::value) { if (n_vec[i] != 0) { - logp += binomial_coefficient_log(n_vec[i] + alpha_dbl - 1.0, - alpha_dbl - 1.0); + logp += binomial_coefficient_log(n_vec[i] + alpha_dbl[i] - 1.0, + alpha_dbl[i] - 1.0); } } - logp += -alpha_dbl * log1p_inv_beta[i] - n_vec[i] * log1p_beta[i]; + logp -= alpha_dbl[i] * log1p_inv_beta[i] + n_vec[i] * log1p_beta[i]; if (!is_constant_all::value) { - ops_partials.edge1_.partials_[i] += digamma(alpha_dbl + n_vec[i]) + ops_partials.edge1_.partials_[i] += digamma(alpha_dbl[i] + n_vec[i]) - digamma_alpha[i] - log1p_inv_beta[i]; } if (!is_constant_all::value) { ops_partials.edge2_.partials_[i] - += lambda_m_alpha_over_1p_beta[i] - n_vec[i] / (beta_dbl + 1.0); + += lambda_m_alpha_over_1p_beta[i] - n_vec[i] / (beta_dbl[i] + 1.0); } } From c1822fcfe3720018076b3b69e6c6ad43dc22ac5d Mon Sep 17 00:00:00 2001 From: Ben Date: Tue, 21 Jul 2020 13:55:33 -0400 Subject: [PATCH 15/16] Switching back to logs of value_ofs in neg_binomial. Adding some changes I forgot previously (Issue #1763) --- stan/math/prim/prob/neg_binomial_lpmf.hpp | 27 ++++++++++--------- stan/math/prim/prob/pareto_type_2_lpdf.hpp | 7 +++-- test/unit/math/mix/prob/neg_binomial_test.cpp | 16 +++++------ 3 files changed, 25 insertions(+), 25 deletions(-) diff --git a/stan/math/prim/prob/neg_binomial_lpmf.hpp b/stan/math/prim/prob/neg_binomial_lpmf.hpp index df2682387f8..901bc9a6e64 100644 --- a/stan/math/prim/prob/neg_binomial_lpmf.hpp +++ b/stan/math/prim/prob/neg_binomial_lpmf.hpp @@ -57,22 +57,20 @@ return_type_t neg_binomial_lpmf(const T_n& n, size_t size_alpha_beta = max_size(alpha, beta); size_t max_size_seq_view = max_size(n, alpha, beta); - scalar_seq_view alpha_dbl(value_of(alpha)); - scalar_seq_view beta_dbl(value_of(beta)); - VectorBuilder::value, T_partials_return, T_shape> digamma_alpha(size_alpha); if (!is_constant_all::value) { for (size_t i = 0; i < size_alpha; ++i) { - digamma_alpha[i] = digamma(alpha_dbl[i]); + digamma_alpha[i] = digamma(value_of(alpha_vec[i])); } } VectorBuilder log1p_inv_beta(size_beta); VectorBuilder log1p_beta(size_beta); for (size_t i = 0; i < size_beta; ++i) { - log1p_inv_beta[i] = log1p(inv(beta_dbl[i])); - log1p_beta[i] = log1p(beta_dbl[i]); + const T_partials_return beta_dbl = value_of(beta_vec[i]); + log1p_inv_beta[i] = log1p(inv(beta_dbl)); + log1p_beta[i] = log1p(beta_dbl); } VectorBuilder::value, T_partials_return, @@ -80,28 +78,33 @@ return_type_t neg_binomial_lpmf(const T_n& n, lambda_m_alpha_over_1p_beta(size_alpha_beta); if (!is_constant_all::value) { for (size_t i = 0; i < size_alpha_beta; ++i) { + const T_partials_return alpha_dbl = value_of(alpha_vec[i]); + const T_partials_return beta_dbl = value_of(beta_vec[i]); lambda_m_alpha_over_1p_beta[i] - = alpha_dbl[i] / beta_dbl[i] - alpha_dbl[i] / (1 + beta_dbl[i]); + = alpha_dbl / beta_dbl - alpha_dbl / (1 + beta_dbl); } } for (size_t i = 0; i < max_size_seq_view; i++) { + const T_partials_return alpha_dbl = value_of(alpha_vec[i]); + const T_partials_return beta_dbl = value_of(beta_vec[i]); + if (include_summand::value) { if (n_vec[i] != 0) { - logp += binomial_coefficient_log(n_vec[i] + alpha_dbl[i] - 1.0, - alpha_dbl[i] - 1.0); + logp += binomial_coefficient_log(n_vec[i] + alpha_dbl - 1.0, + alpha_dbl - 1.0); } } - logp -= alpha_dbl[i] * log1p_inv_beta[i] + n_vec[i] * log1p_beta[i]; + logp -= alpha_dbl * log1p_inv_beta[i] + n_vec[i] * log1p_beta[i]; if (!is_constant_all::value) { - ops_partials.edge1_.partials_[i] += digamma(alpha_dbl[i] + n_vec[i]) + ops_partials.edge1_.partials_[i] += digamma(alpha_dbl + n_vec[i]) - digamma_alpha[i] - log1p_inv_beta[i]; } if (!is_constant_all::value) { ops_partials.edge2_.partials_[i] - += lambda_m_alpha_over_1p_beta[i] - n_vec[i] / (beta_dbl[i] + 1.0); + += lambda_m_alpha_over_1p_beta[i] - n_vec[i] / (beta_dbl + 1.0); } } diff --git a/stan/math/prim/prob/pareto_type_2_lpdf.hpp b/stan/math/prim/prob/pareto_type_2_lpdf.hpp index 7a08527ec14..eb0e1789f9a 100644 --- a/stan/math/prim/prob/pareto_type_2_lpdf.hpp +++ b/stan/math/prim/prob/pareto_type_2_lpdf.hpp @@ -67,12 +67,11 @@ return_type_t pareto_type_2_lpdf( } } - VectorBuilder::value, T_partials_return, - T_shape> + VectorBuilder::value, T_partials_return, T_shape> inv_alpha(size(alpha)); - if (include_summand::value) { + if (!is_constant_all::value) { for (size_t n = 0; n < stan::math::size(alpha); n++) { - inv_alpha[n] = 1 / value_of(alpha_vec[n]); + inv_alpha[n] = inv(value_of(alpha_vec[n])); } } diff --git a/test/unit/math/mix/prob/neg_binomial_test.cpp b/test/unit/math/mix/prob/neg_binomial_test.cpp index 03326fcd25c..88223a3b1e8 100644 --- a/test/unit/math/mix/prob/neg_binomial_test.cpp +++ b/test/unit/math/mix/prob/neg_binomial_test.cpp @@ -4,17 +4,15 @@ #include TEST(mathMixScalFun, neg_binomial_lpmf_derivatives) { - auto f1 = [](const auto& alpha, const auto& beta) { - return stan::math::neg_binomial_lpmf(0, alpha, beta); - }; - auto f2 = [](const auto& alpha, const auto& beta) { - return stan::math::neg_binomial_lpmf(6, alpha, beta); + auto f = [](const int y) { + return + [=](const auto& alpha, const auto& beta) { return stan::math::neg_binomial_lpmf(y, alpha, beta); }; }; - stan::test::expect_ad(f1, 1.5, 4.1); - stan::test::expect_ad(f1, std::vector({1.2, 2.0}), + stan::test::expect_ad(f(0), 1.5, 4.1); + stan::test::expect_ad(f(0), std::vector({1.2, 2.0}), std::vector({1.1, 1.2})); - stan::test::expect_ad(f2, 1.5, 4.1); - stan::test::expect_ad(f2, std::vector({1.7, 2.0}), + stan::test::expect_ad(f(6), 1.5, 4.1); + stan::test::expect_ad(f(6), std::vector({1.7, 2.0}), std::vector({1.1, 2.3})); } From 55cf763fede2f4bc24a6f92e9e972c3a0e0862a0 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Tue, 21 Jul 2020 19:58:55 +0000 Subject: [PATCH 16/16] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- test/unit/math/mix/prob/neg_binomial_test.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/unit/math/mix/prob/neg_binomial_test.cpp b/test/unit/math/mix/prob/neg_binomial_test.cpp index 88223a3b1e8..f106f51ab69 100644 --- a/test/unit/math/mix/prob/neg_binomial_test.cpp +++ b/test/unit/math/mix/prob/neg_binomial_test.cpp @@ -5,8 +5,9 @@ TEST(mathMixScalFun, neg_binomial_lpmf_derivatives) { auto f = [](const int y) { - return - [=](const auto& alpha, const auto& beta) { return stan::math::neg_binomial_lpmf(y, alpha, beta); }; + return [=](const auto& alpha, const auto& beta) { + return stan::math::neg_binomial_lpmf(y, alpha, beta); + }; }; stan::test::expect_ad(f(0), 1.5, 4.1);