From 1ec777a52c64e02d44fbefcb912aa645b5ee316c Mon Sep 17 00:00:00 2001 From: Marco Colombo Date: Fri, 31 Jan 2020 15:50:12 +0000 Subject: [PATCH 1/9] Loop over the correct number of elements when constructing VectorBuilders --- stan/math/prim/prob/neg_binomial_2_cdf.hpp | 18 +++--- .../prim/prob/neg_binomial_2_log_lpmf.hpp | 56 +++++++++--------- stan/math/prim/prob/neg_binomial_2_lpmf.hpp | 57 +++++++++---------- stan/math/prim/prob/neg_binomial_cdf.hpp | 16 +++--- stan/math/prim/prob/neg_binomial_lccdf.hpp | 37 +++++++----- stan/math/prim/prob/neg_binomial_lcdf.hpp | 37 +++++++----- 6 files changed, 116 insertions(+), 105 deletions(-) diff --git a/stan/math/prim/prob/neg_binomial_2_cdf.hpp b/stan/math/prim/prob/neg_binomial_2_cdf.hpp index a080ebe55cb..5cdaa662d90 100644 --- a/stan/math/prim/prob/neg_binomial_2_cdf.hpp +++ b/stan/math/prim/prob/neg_binomial_2_cdf.hpp @@ -34,6 +34,8 @@ return_type_t neg_binomial_2_cdf( scalar_seq_view n_vec(n); scalar_seq_view mu_vec(mu); scalar_seq_view phi_vec(phi); + size_t size_phi = size(phi); + size_t size_n_phi = max_size(n, phi); size_t max_size_seq_view = max_size(n, mu, phi); operands_and_partials ops_partials(mu, phi); @@ -48,18 +50,18 @@ return_type_t neg_binomial_2_cdf( VectorBuilder::value, T_partials_return, T_precision> - digamma_phi_vec(size(phi)); - - VectorBuilder::value, T_partials_return, + digamma_phi_vec(size_phi); + VectorBuilder::value, T_partials_return, T_n, T_precision> - digamma_sum_vec(size(phi)); + digamma_sum_vec(size_n_phi); if (!is_constant_all::value) { - for (size_t i = 0; i < size(phi); i++) { + for (size_t i = 0; i < size_phi; i++) { + digamma_phi_vec[i] = digamma(value_of(phi_vec[i])); + } + for (size_t i = 0; i < size_n_phi; i++) { const T_partials_return n_dbl = value_of(n_vec[i]); const T_partials_return phi_dbl = value_of(phi_vec[i]); - - digamma_phi_vec[i] = digamma(phi_dbl); digamma_sum_vec[i] = digamma(n_dbl + phi_dbl + 1); } } @@ -104,7 +106,7 @@ return_type_t neg_binomial_2_cdf( } if (!is_constant_all::value) { - for (size_t i = 0; i < size(phi); ++i) { + for (size_t i = 0; i < size_phi; ++i) { ops_partials.edge2_.partials_[i] *= P; } } diff --git a/stan/math/prim/prob/neg_binomial_2_log_lpmf.hpp b/stan/math/prim/prob/neg_binomial_2_log_lpmf.hpp index 30162f83565..06a376efaa5 100644 --- a/stan/math/prim/prob/neg_binomial_2_log_lpmf.hpp +++ b/stan/math/prim/prob/neg_binomial_2_log_lpmf.hpp @@ -48,58 +48,54 @@ return_type_t neg_binomial_2_log_lpmf( scalar_seq_view n_vec(n); scalar_seq_view eta_vec(eta); scalar_seq_view phi_vec(phi); + size_t size_eta = size(eta); + size_t size_phi = size(phi); + size_t size_eta_phi = max_size(eta, phi); + size_t size_n_phi = max_size(n, phi); size_t max_size_seq_view = max_size(n, eta, phi); operands_and_partials ops_partials(eta, phi); - size_t len_ep = max_size(eta, phi); - size_t len_np = max_size(n, phi); - - VectorBuilder eta__(size(eta)); - for (size_t i = 0, max_size_seq_view = size(eta); i < max_size_seq_view; - ++i) { - eta__[i] = value_of(eta_vec[i]); - } - - VectorBuilder phi__(size(phi)); - for (size_t i = 0, max_size_seq_view = size(phi); i < max_size_seq_view; - ++i) { - phi__[i] = value_of(phi_vec[i]); + VectorBuilder eta_val(size_eta); + for (size_t i = 0; i < size_eta; ++i) { + eta_val[i] = value_of(eta_vec[i]); } - VectorBuilder log_phi(size(phi)); - for (size_t i = 0, max_size_seq_view = size(phi); i < max_size_seq_view; - ++i) { - log_phi[i] = log(phi__[i]); + VectorBuilder phi_val(size_phi); + VectorBuilder log_phi(size_phi); + for (size_t i = 0; i < size_phi; ++i) { + phi_val[i] = value_of(phi_vec[i]); + log_phi[i] = log(phi_val[i]); } VectorBuilder - logsumexp_eta_logphi(len_ep); - for (size_t i = 0; i < len_ep; ++i) { - logsumexp_eta_logphi[i] = log_sum_exp(eta__[i], log_phi[i]); + logsumexp_eta_logphi(size_eta_phi); + for (size_t i = 0; i < size_eta_phi; ++i) { + logsumexp_eta_logphi[i] = log_sum_exp(eta_val[i], log_phi[i]); } - VectorBuilder n_plus_phi(len_np); - for (size_t i = 0; i < len_np; ++i) { - n_plus_phi[i] = n_vec[i] + phi__[i]; + VectorBuilder n_plus_phi( + size_n_phi); + for (size_t i = 0; i < size_n_phi; ++i) { + n_plus_phi[i] = n_vec[i] + phi_val[i]; } for (size_t i = 0; i < max_size_seq_view; i++) { - if (phi__[i] > 1e5) { + if (phi_val[i] > 1e5) { // TODO(martinmodrak) This is wrong (doesn't pass propto information), // and inaccurate for n = 0, but shouldn't break most models. // Also the 1e5 cutoff is way too low. // Will be adressed better once PR #1497 is merged - logp += poisson_log_lpmf(n_vec[i], eta__[i]); + logp += poisson_log_lpmf(n_vec[i], eta_val[i]); } else { if (include_summand::value) { logp -= lgamma(n_vec[i] + 1.0); } if (include_summand::value) { - logp += multiply_log(phi__[i], phi__[i]) - lgamma(phi__[i]); + logp += multiply_log(phi_val[i], phi_val[i]) - lgamma(phi_val[i]); } if (include_summand::value) { - logp += n_vec[i] * eta__[i]; + logp += n_vec[i] * eta_val[i]; } if (include_summand::value) { logp += lgamma(n_plus_phi[i]); @@ -109,12 +105,12 @@ return_type_t neg_binomial_2_log_lpmf( if (!is_constant_all::value) { ops_partials.edge1_.partials_[i] - += n_vec[i] - n_plus_phi[i] / (phi__[i] / exp(eta__[i]) + 1.0); + += n_vec[i] - n_plus_phi[i] / (phi_val[i] / exp(eta_val[i]) + 1.0); } if (!is_constant_all::value) { ops_partials.edge2_.partials_[i] - += 1.0 - n_plus_phi[i] / (exp(eta__[i]) + phi__[i]) + log_phi[i] - - logsumexp_eta_logphi[i] - digamma(phi__[i]) + += 1.0 - n_plus_phi[i] / (exp(eta_val[i]) + phi_val[i]) + log_phi[i] + - logsumexp_eta_logphi[i] - digamma(phi_val[i]) + digamma(n_plus_phi[i]); } } diff --git a/stan/math/prim/prob/neg_binomial_2_lpmf.hpp b/stan/math/prim/prob/neg_binomial_2_lpmf.hpp index a2efdaf25b7..c1fa5365570 100644 --- a/stan/math/prim/prob/neg_binomial_2_lpmf.hpp +++ b/stan/math/prim/prob/neg_binomial_2_lpmf.hpp @@ -42,59 +42,55 @@ return_type_t neg_binomial_2_lpmf( scalar_seq_view n_vec(n); scalar_seq_view mu_vec(mu); scalar_seq_view phi_vec(phi); + size_t size_mu = size(mu); + size_t size_phi = size(phi); + size_t size_mu_phi = max_size(mu, phi); + size_t size_n_phi = max_size(n, phi); size_t max_size_seq_view = max_size(n, mu, phi); operands_and_partials ops_partials(mu, phi); - size_t len_ep = max_size(mu, phi); - size_t len_np = max_size(n, phi); - - VectorBuilder mu__(size(mu)); - - for (size_t i = 0, max_size_seq_view = size(mu); i < max_size_seq_view; ++i) { - mu__[i] = value_of(mu_vec[i]); - } - - VectorBuilder phi__(size(phi)); - for (size_t i = 0, max_size_seq_view = size(phi); i < max_size_seq_view; - ++i) { - phi__[i] = value_of(phi_vec[i]); + VectorBuilder mu_val(size_mu); + for (size_t i = 0; i < size_mu; ++i) { + mu_val[i] = value_of(mu_vec[i]); } - VectorBuilder log_phi(size(phi)); - for (size_t i = 0, max_size_seq_view = size(phi); i < max_size_seq_view; - ++i) { - log_phi[i] = log(phi__[i]); + VectorBuilder phi_val(size_phi); + VectorBuilder log_phi(size_phi); + for (size_t i = 0; i < size_phi; ++i) { + phi_val[i] = value_of(phi_vec[i]); + log_phi[i] = log(phi_val[i]); } VectorBuilder - log_mu_plus_phi(len_ep); - for (size_t i = 0; i < len_ep; ++i) { - log_mu_plus_phi[i] = log(mu__[i] + phi__[i]); + log_mu_plus_phi(size_mu_phi); + for (size_t i = 0; i < size_mu_phi; ++i) { + log_mu_plus_phi[i] = log(mu_val[i] + phi_val[i]); } - VectorBuilder n_plus_phi(len_np); - for (size_t i = 0; i < len_np; ++i) { - n_plus_phi[i] = n_vec[i] + phi__[i]; + VectorBuilder n_plus_phi( + size_n_phi); + for (size_t i = 0; i < size_n_phi; ++i) { + n_plus_phi[i] = n_vec[i] + phi_val[i]; } for (size_t i = 0; i < max_size_seq_view; i++) { // if phi is large we probably overflow, defer to Poisson: - if (phi__[i] > 1e5) { + if (phi_val[i] > 1e5) { // TODO(martinmodrak) This is wrong (doesn't pass propto information), // and inaccurate for n = 0, but shouldn't break most models. // Also the 1e5 cutoff is too small. // Will be adressed better in PR #1497 - logp += poisson_lpmf(n_vec[i], mu__[i]); + logp += poisson_lpmf(n_vec[i], mu_val[i]); } else { if (include_summand::value) { logp -= lgamma(n_vec[i] + 1.0); } if (include_summand::value) { - logp += multiply_log(phi__[i], phi__[i]) - lgamma(phi__[i]); + logp += multiply_log(phi_val[i], phi_val[i]) - lgamma(phi_val[i]); } if (include_summand::value) { - logp += multiply_log(n_vec[i], mu__[i]); + logp += multiply_log(n_vec[i], mu_val[i]); } if (include_summand::value) { logp += lgamma(n_plus_phi[i]); @@ -104,12 +100,13 @@ return_type_t neg_binomial_2_lpmf( if (!is_constant_all::value) { ops_partials.edge1_.partials_[i] - += n_vec[i] / mu__[i] - (n_vec[i] + phi__[i]) / (mu__[i] + phi__[i]); + += n_vec[i] / mu_val[i] - n_plus_phi[i] / (mu_val[i] + phi_val[i]); } if (!is_constant_all::value) { ops_partials.edge2_.partials_[i] - += 1.0 - n_plus_phi[i] / (mu__[i] + phi__[i]) + log_phi[i] - - log_mu_plus_phi[i] - digamma(phi__[i]) + digamma(n_plus_phi[i]); + += 1.0 - n_plus_phi[i] / (mu_val[i] + phi_val[i]) + log_phi[i] + - log_mu_plus_phi[i] - digamma(phi_val[i]) + + digamma(n_plus_phi[i]); } } return ops_partials.build(logp); diff --git a/stan/math/prim/prob/neg_binomial_cdf.hpp b/stan/math/prim/prob/neg_binomial_cdf.hpp index 0dca1aae2bf..0af303da6af 100644 --- a/stan/math/prim/prob/neg_binomial_cdf.hpp +++ b/stan/math/prim/prob/neg_binomial_cdf.hpp @@ -36,6 +36,8 @@ return_type_t neg_binomial_cdf(const T_n& n, scalar_seq_view n_vec(n); scalar_seq_view alpha_vec(alpha); scalar_seq_view beta_vec(beta); + size_t size_alpha = size(alpha); + size_t size_n_alpha = max_size(n, alpha); size_t max_size_seq_view = max_size(n, alpha, beta); operands_and_partials ops_partials(alpha, beta); @@ -49,17 +51,17 @@ return_type_t neg_binomial_cdf(const T_n& n, } VectorBuilder::value, T_partials_return, T_shape> - digamma_alpha_vec(size(alpha)); - + digamma_alpha_vec(size_alpha); VectorBuilder::value, T_partials_return, T_shape> - digamma_sum_vec(size(alpha)); + digamma_sum_vec(size_n_alpha); if (!is_constant_all::value) { - for (size_t i = 0; i < size(alpha); i++) { + for (size_t i = 0; i < size_alpha; i++) { + digamma_alpha_vec[i] = digamma(value_of(alpha_vec[i])); + } + for (size_t i = 0; i < size_n_alpha; i++) { const T_partials_return n_dbl = value_of(n_vec[i]); const T_partials_return alpha_dbl = value_of(alpha_vec[i]); - - digamma_alpha_vec[i] = digamma(alpha_dbl); digamma_sum_vec[i] = digamma(n_dbl + alpha_dbl + 1); } } @@ -96,7 +98,7 @@ return_type_t neg_binomial_cdf(const T_n& n, } if (!is_constant_all::value) { - for (size_t i = 0; i < size(alpha); ++i) { + for (size_t i = 0; i < size_alpha; ++i) { ops_partials.edge1_.partials_[i] *= P; } } diff --git a/stan/math/prim/prob/neg_binomial_lccdf.hpp b/stan/math/prim/prob/neg_binomial_lccdf.hpp index 15323fe2a78..243662ed968 100644 --- a/stan/math/prim/prob/neg_binomial_lccdf.hpp +++ b/stan/math/prim/prob/neg_binomial_lccdf.hpp @@ -33,39 +33,46 @@ return_type_t neg_binomial_lccdf( check_consistent_sizes(function, "Failures variable", n, "Shape parameter", alpha, "Inverse scale parameter", beta); - scalar_seq_view n_vec(n); - scalar_seq_view alpha_vec(alpha); - scalar_seq_view beta_vec(beta); - size_t max_size_seq_view = max_size(n, alpha, beta); - using std::exp; using std::log; using std::pow; operands_and_partials ops_partials(alpha, beta); + scalar_seq_view n_vec(n); + scalar_seq_view alpha_vec(alpha); + scalar_seq_view beta_vec(beta); + size_t size_n = size(n); + size_t size_alpha = size(alpha); + size_t size_n_alpha = max_size(n, alpha); + size_t max_size_seq_view = max_size(n, alpha, beta); + // Explicit return for extreme values // The gradients are technically ill-defined, but treated as zero - for (size_t i = 0; i < size(n); i++) { + for (size_t i = 0; i < size_n; i++) { if (value_of(n_vec[i]) < 0) { return ops_partials.build(0.0); } } + VectorBuilder::value, T_partials_return, T_n> + digammaN_vec(size_n); VectorBuilder::value, T_partials_return, T_shape> - digammaN_vec(size(alpha)); - VectorBuilder::value, T_partials_return, T_shape> - digammaAlpha_vec(size(alpha)); - VectorBuilder::value, T_partials_return, T_shape> - digammaSum_vec(size(alpha)); + digammaAlpha_vec(size_alpha); + VectorBuilder::value, T_partials_return, T_n, + T_shape> + digammaSum_vec(size_n_alpha); if (!is_constant_all::value) { - for (size_t i = 0; i < size(alpha); i++) { + for (size_t i = 0; i < size_n; i++) { + digammaN_vec[i] = digamma(value_of(n_vec[i]) + 1); + } + for (size_t i = 0; i < size_alpha; i++) { + digammaAlpha_vec[i] = digamma(value_of(alpha_vec[i])); + } + for (size_t i = 0; i < size_n_alpha; i++) { const T_partials_return n_dbl = value_of(n_vec[i]); const T_partials_return alpha_dbl = value_of(alpha_vec[i]); - - digammaN_vec[i] = digamma(n_dbl + 1); - digammaAlpha_vec[i] = digamma(alpha_dbl); digammaSum_vec[i] = digamma(n_dbl + alpha_dbl + 1); } } diff --git a/stan/math/prim/prob/neg_binomial_lcdf.hpp b/stan/math/prim/prob/neg_binomial_lcdf.hpp index 30a10433f18..35a7affb303 100644 --- a/stan/math/prim/prob/neg_binomial_lcdf.hpp +++ b/stan/math/prim/prob/neg_binomial_lcdf.hpp @@ -34,39 +34,46 @@ return_type_t neg_binomial_lcdf(const T_n& n, check_consistent_sizes(function, "Failures variable", n, "Shape parameter", alpha, "Inverse scale parameter", beta); - scalar_seq_view n_vec(n); - scalar_seq_view alpha_vec(alpha); - scalar_seq_view beta_vec(beta); - size_t max_size_seq_view = max_size(n, alpha, beta); - using std::exp; using std::log; using std::pow; operands_and_partials ops_partials(alpha, beta); + scalar_seq_view n_vec(n); + scalar_seq_view alpha_vec(alpha); + scalar_seq_view beta_vec(beta); + size_t size_n = size(n); + size_t size_alpha = size(alpha); + size_t size_n_alpha = max_size(n, alpha); + size_t max_size_seq_view = max_size(n, alpha, beta); + // Explicit return for extreme values // The gradients are technically ill-defined, but treated as zero - for (size_t i = 0; i < size(n); i++) { + for (size_t i = 0; i < size_n; i++) { if (value_of(n_vec[i]) < 0) { return ops_partials.build(negative_infinity()); } } + VectorBuilder::value, T_partials_return, T_n> + digammaN_vec(size_n); VectorBuilder::value, T_partials_return, T_shape> - digammaN_vec(size(alpha)); - VectorBuilder::value, T_partials_return, T_shape> - digammaAlpha_vec(size(alpha)); - VectorBuilder::value, T_partials_return, T_shape> - digammaSum_vec(size(alpha)); + digammaAlpha_vec(size_alpha); + VectorBuilder::value, T_partials_return, T_n, + T_shape> + digammaSum_vec(size_n_alpha); if (!is_constant_all::value) { - for (size_t i = 0; i < size(alpha); i++) { + for (size_t i = 0; i < size_n; i++) { + digammaN_vec[i] = digamma(value_of(n_vec[i]) + 1); + } + for (size_t i = 0; i < size_alpha; i++) { + digammaAlpha_vec[i] = digamma(value_of(alpha_vec[i])); + } + for (size_t i = 0; i < size_n_alpha; i++) { const T_partials_return n_dbl = value_of(n_vec[i]); const T_partials_return alpha_dbl = value_of(alpha_vec[i]); - - digammaN_vec[i] = digamma(n_dbl + 1); - digammaAlpha_vec[i] = digamma(alpha_dbl); digammaSum_vec[i] = digamma(n_dbl + alpha_dbl + 1); } } From dbaac119e3943acb4ff802cd20b105da6614ed5e Mon Sep 17 00:00:00 2001 From: Marco Colombo Date: Fri, 31 Jan 2020 18:21:33 +0000 Subject: [PATCH 2/9] Avoid repeated computations of the same quantities --- stan/math/prim/prob/neg_binomial_2_cdf.hpp | 27 ++++--- stan/math/prim/prob/neg_binomial_2_lccdf.hpp | 3 +- stan/math/prim/prob/neg_binomial_2_lcdf.hpp | 17 ++-- .../prim/prob/neg_binomial_2_log_lpmf.hpp | 13 +++- stan/math/prim/prob/neg_binomial_2_lpmf.hpp | 21 ++--- stan/math/prim/prob/neg_binomial_cdf.hpp | 8 +- stan/math/prim/prob/neg_binomial_lccdf.hpp | 15 ++-- stan/math/prim/prob/neg_binomial_lcdf.hpp | 15 ++-- stan/math/prim/prob/neg_binomial_lpmf.hpp | 77 ++++++++----------- 9 files changed, 104 insertions(+), 92 deletions(-) diff --git a/stan/math/prim/prob/neg_binomial_2_cdf.hpp b/stan/math/prim/prob/neg_binomial_2_cdf.hpp index 5cdaa662d90..dad7a41555c 100644 --- a/stan/math/prim/prob/neg_binomial_2_cdf.hpp +++ b/stan/math/prim/prob/neg_binomial_2_cdf.hpp @@ -3,11 +3,13 @@ #include #include -#include #include #include #include #include +#include +#include +#include #include #include @@ -73,29 +75,30 @@ return_type_t neg_binomial_2_cdf( return ops_partials.build(1.0); } - const T_partials_return n_dbl = value_of(n_vec[i]); + const T_partials_return n_dbl_p1 = value_of(n_vec[i]) + 1; const T_partials_return mu_dbl = value_of(mu_vec[i]); const T_partials_return phi_dbl = value_of(phi_vec[i]); - - const T_partials_return p_dbl = phi_dbl / (mu_dbl + phi_dbl); - const T_partials_return d_dbl - = 1.0 / ((mu_dbl + phi_dbl) * (mu_dbl + phi_dbl)); - - const T_partials_return P_i = inc_beta(phi_dbl, n_dbl + 1.0, p_dbl); + const T_partials_return inv_mu_plus_phi = inv(mu_dbl + phi_dbl); + const T_partials_return p_dbl = phi_dbl * inv_mu_plus_phi; + const T_partials_return d_dbl = square(inv_mu_plus_phi); + const T_partials_return P_i = inc_beta(phi_dbl, n_dbl_p1, p_dbl); + const T_partials_return inc_beta_ddz_i + = is_constant_all::value + ? 0 + : inc_beta_ddz(phi_dbl, n_dbl_p1, p_dbl) * d_dbl / P_i; P *= P_i; if (!is_constant_all::value) { - ops_partials.edge1_.partials_[i] - += -inc_beta_ddz(phi_dbl, n_dbl + 1.0, p_dbl) * phi_dbl * d_dbl / P_i; + ops_partials.edge1_.partials_[i] -= inc_beta_ddz_i * phi_dbl; } if (!is_constant_all::value) { ops_partials.edge2_.partials_[i] - += inc_beta_dda(phi_dbl, n_dbl + 1, p_dbl, digamma_phi_vec[i], + += inc_beta_dda(phi_dbl, n_dbl_p1, p_dbl, digamma_phi_vec[i], digamma_sum_vec[i]) / P_i - + inc_beta_ddz(phi_dbl, n_dbl + 1.0, p_dbl) * mu_dbl * d_dbl / P_i; + + inc_beta_ddz_i * mu_dbl; } } diff --git a/stan/math/prim/prob/neg_binomial_2_lccdf.hpp b/stan/math/prim/prob/neg_binomial_2_lccdf.hpp index 36a34b06f2b..2a2c383b199 100644 --- a/stan/math/prim/prob/neg_binomial_2_lccdf.hpp +++ b/stan/math/prim/prob/neg_binomial_2_lccdf.hpp @@ -27,7 +27,6 @@ return_type_t neg_binomial_2_lccdf( scalar_seq_view mu_vec(mu); scalar_seq_view phi_vec(phi); - size_t size_beta = max_size(mu, phi); VectorBuilder, T_location, @@ -37,7 +36,7 @@ return_type_t neg_binomial_2_lccdf( beta_vec[i] = phi_vec[i] / mu_vec[i]; } - return neg_binomial_ccdf_log(n, phi, beta_vec.data()); + return neg_binomial_lccdf(n, phi, beta_vec.data()); } } // namespace math diff --git a/stan/math/prim/prob/neg_binomial_2_lcdf.hpp b/stan/math/prim/prob/neg_binomial_2_lcdf.hpp index 966e1102f9a..f4f36ff9282 100644 --- a/stan/math/prim/prob/neg_binomial_2_lcdf.hpp +++ b/stan/math/prim/prob/neg_binomial_2_lcdf.hpp @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -29,8 +30,15 @@ return_type_t neg_binomial_2_lcdf( scalar_seq_view n_vec(n); scalar_seq_view mu_vec(mu); scalar_seq_view phi_vec(phi); - + size_t size_n = size(n); size_t size_phi_mu = max_size(mu, phi); + + for (size_t i = 0; i < size_n; i++) { + if (n_vec[i] < 0) { + return LOG_ZERO; + } + } + VectorBuilder, T_location, T_precision> phi_mu(size_phi_mu); @@ -38,14 +46,9 @@ return_type_t neg_binomial_2_lcdf( phi_mu[i] = phi_vec[i] / (phi_vec[i] + mu_vec[i]); } - size_t size_n = size(n); VectorBuilder, T_n> np1(size_n); for (size_t i = 0; i < size_n; i++) { - if (n_vec[i] < 0) { - return log(0.0); - } else { - np1[i] = n_vec[i] + 1.0; - } + np1[i] = n_vec[i] + 1.0; } return beta_cdf_log(phi_mu.data(), phi, np1.data()); diff --git a/stan/math/prim/prob/neg_binomial_2_log_lpmf.hpp b/stan/math/prim/prob/neg_binomial_2_log_lpmf.hpp index 06a376efaa5..3989e1da857 100644 --- a/stan/math/prim/prob/neg_binomial_2_log_lpmf.hpp +++ b/stan/math/prim/prob/neg_binomial_2_log_lpmf.hpp @@ -61,6 +61,15 @@ return_type_t neg_binomial_2_log_lpmf( eta_val[i] = value_of(eta_vec[i]); } + VectorBuilder::value, + T_partials_return, T_log_location> + exp_eta(size_eta); + if (!is_constant_all::value) { + for (size_t i = 0; i < size_eta; ++i) { + exp_eta[i] = exp(eta_val[i]); + } + } + VectorBuilder phi_val(size_phi); VectorBuilder log_phi(size_phi); for (size_t i = 0; i < size_phi; ++i) { @@ -105,11 +114,11 @@ return_type_t neg_binomial_2_log_lpmf( if (!is_constant_all::value) { ops_partials.edge1_.partials_[i] - += n_vec[i] - n_plus_phi[i] / (phi_val[i] / exp(eta_val[i]) + 1.0); + += n_vec[i] - n_plus_phi[i] / (phi_val[i] / exp_eta[i] + 1.0); } if (!is_constant_all::value) { ops_partials.edge2_.partials_[i] - += 1.0 - n_plus_phi[i] / (exp(eta_val[i]) + phi_val[i]) + log_phi[i] + += 1.0 - n_plus_phi[i] / (exp_eta[i] + phi_val[i]) + log_phi[i] - logsumexp_eta_logphi[i] - digamma(phi_val[i]) + digamma(n_plus_phi[i]); } diff --git a/stan/math/prim/prob/neg_binomial_2_lpmf.hpp b/stan/math/prim/prob/neg_binomial_2_lpmf.hpp index c1fa5365570..1fe639fa771 100644 --- a/stan/math/prim/prob/neg_binomial_2_lpmf.hpp +++ b/stan/math/prim/prob/neg_binomial_2_lpmf.hpp @@ -3,10 +3,10 @@ #include #include -#include -#include #include #include +#include +#include #include #include #include @@ -62,10 +62,13 @@ return_type_t neg_binomial_2_lpmf( log_phi[i] = log(phi_val[i]); } + VectorBuilder mu_plus_phi( + size_mu_phi); VectorBuilder log_mu_plus_phi(size_mu_phi); for (size_t i = 0; i < size_mu_phi; ++i) { - log_mu_plus_phi[i] = log(mu_val[i] + phi_val[i]); + mu_plus_phi[i] = mu_val[i] + phi_val[i]; + log_mu_plus_phi[i] = log(mu_plus_phi[i]); } VectorBuilder n_plus_phi( @@ -95,18 +98,18 @@ return_type_t neg_binomial_2_lpmf( if (include_summand::value) { logp += lgamma(n_plus_phi[i]); } - logp -= (n_plus_phi[i]) * log_mu_plus_phi[i]; + logp -= n_plus_phi[i] * log_mu_plus_phi[i]; } if (!is_constant_all::value) { ops_partials.edge1_.partials_[i] - += n_vec[i] / mu_val[i] - n_plus_phi[i] / (mu_val[i] + phi_val[i]); + += n_vec[i] / mu_val[i] - n_plus_phi[i] / mu_plus_phi[i]; } if (!is_constant_all::value) { - ops_partials.edge2_.partials_[i] - += 1.0 - n_plus_phi[i] / (mu_val[i] + phi_val[i]) + log_phi[i] - - log_mu_plus_phi[i] - digamma(phi_val[i]) - + digamma(n_plus_phi[i]); + ops_partials.edge2_.partials_[i] += 1.0 - n_plus_phi[i] / mu_plus_phi[i] + + log_phi[i] - log_mu_plus_phi[i] + - digamma(phi_val[i]) + + digamma(n_plus_phi[i]); } } return ops_partials.build(logp); diff --git a/stan/math/prim/prob/neg_binomial_cdf.hpp b/stan/math/prim/prob/neg_binomial_cdf.hpp index 0af303da6af..0a5d77d5c02 100644 --- a/stan/math/prim/prob/neg_binomial_cdf.hpp +++ b/stan/math/prim/prob/neg_binomial_cdf.hpp @@ -7,7 +7,9 @@ #include #include #include +#include #include +#include #include #include #include @@ -76,9 +78,9 @@ return_type_t neg_binomial_cdf(const T_n& n, const T_partials_return n_dbl = value_of(n_vec[i]); const T_partials_return alpha_dbl = value_of(alpha_vec[i]); const T_partials_return beta_dbl = value_of(beta_vec[i]); - - const T_partials_return p_dbl = beta_dbl / (1.0 + beta_dbl); - const T_partials_return d_dbl = 1.0 / ((1.0 + beta_dbl) * (1.0 + beta_dbl)); + const T_partials_return inv_beta_p1 = inv(beta_dbl + 1); + const T_partials_return p_dbl = beta_dbl * inv_beta_p1; + const T_partials_return d_dbl = square(inv_beta_p1); const T_partials_return P_i = inc_beta(alpha_dbl, n_dbl + 1.0, p_dbl); diff --git a/stan/math/prim/prob/neg_binomial_lccdf.hpp b/stan/math/prim/prob/neg_binomial_lccdf.hpp index 243662ed968..5152f630f82 100644 --- a/stan/math/prim/prob/neg_binomial_lccdf.hpp +++ b/stan/math/prim/prob/neg_binomial_lccdf.hpp @@ -3,13 +3,15 @@ #include #include -#include +#include #include -#include #include -#include #include #include +#include +#include +#include +#include #include #include @@ -87,8 +89,9 @@ return_type_t neg_binomial_lccdf( const T_partials_return n_dbl = value_of(n_vec[i]); const T_partials_return alpha_dbl = value_of(alpha_vec[i]); const T_partials_return beta_dbl = value_of(beta_vec[i]); - const T_partials_return p_dbl = beta_dbl / (1.0 + beta_dbl); - const T_partials_return d_dbl = 1.0 / ((1.0 + beta_dbl) * (1.0 + beta_dbl)); + const T_partials_return inv_beta_p1 = inv(beta_dbl + 1); + const T_partials_return p_dbl = beta_dbl * inv_beta_p1; + const T_partials_return d_dbl = square(inv_beta_p1); const T_partials_return Pi = 1.0 - inc_beta(alpha_dbl, n_dbl + 1.0, p_dbl); const T_partials_return beta_func = stan::math::beta(n_dbl + 1, alpha_dbl); @@ -106,7 +109,7 @@ return_type_t neg_binomial_lccdf( if (!is_constant_all::value) { ops_partials.edge2_.partials_[i] -= d_dbl * pow(1 - p_dbl, n_dbl) * pow(p_dbl, alpha_dbl - 1) - / beta_func / Pi; + / (beta_func * Pi); } } diff --git a/stan/math/prim/prob/neg_binomial_lcdf.hpp b/stan/math/prim/prob/neg_binomial_lcdf.hpp index 35a7affb303..dc7839ae9a2 100644 --- a/stan/math/prim/prob/neg_binomial_lcdf.hpp +++ b/stan/math/prim/prob/neg_binomial_lcdf.hpp @@ -3,13 +3,15 @@ #include #include -#include +#include #include -#include #include -#include #include #include +#include +#include +#include +#include #include #include @@ -88,8 +90,9 @@ return_type_t neg_binomial_lcdf(const T_n& n, const T_partials_return n_dbl = value_of(n_vec[i]); const T_partials_return alpha_dbl = value_of(alpha_vec[i]); const T_partials_return beta_dbl = value_of(beta_vec[i]); - const T_partials_return p_dbl = beta_dbl / (1.0 + beta_dbl); - const T_partials_return d_dbl = 1.0 / ((1.0 + beta_dbl) * (1.0 + beta_dbl)); + const T_partials_return inv_beta_p1 = inv(beta_dbl + 1); + const T_partials_return p_dbl = beta_dbl * inv_beta_p1; + const T_partials_return d_dbl = square(inv_beta_p1); const T_partials_return Pi = inc_beta(alpha_dbl, n_dbl + 1.0, p_dbl); const T_partials_return beta_func = stan::math::beta(n_dbl + 1, alpha_dbl); @@ -107,7 +110,7 @@ return_type_t neg_binomial_lcdf(const T_n& n, if (!is_constant_all::value) { ops_partials.edge2_.partials_[i] += d_dbl * pow(1 - p_dbl, n_dbl) * pow(p_dbl, alpha_dbl - 1) - / beta_func / Pi; + / (beta_func * Pi); } } diff --git a/stan/math/prim/prob/neg_binomial_lpmf.hpp b/stan/math/prim/prob/neg_binomial_lpmf.hpp index e4fa3fabcce..81667bb195b 100644 --- a/stan/math/prim/prob/neg_binomial_lpmf.hpp +++ b/stan/math/prim/prob/neg_binomial_lpmf.hpp @@ -3,13 +3,13 @@ #include #include -#include -#include -#include #include -#include +#include #include #include +#include +#include +#include #include namespace stan { @@ -44,60 +44,46 @@ return_type_t neg_binomial_lpmf(const T_n& n, scalar_seq_view n_vec(n); scalar_seq_view alpha_vec(alpha); scalar_seq_view beta_vec(beta); + size_t size_alpha = size(alpha); + size_t size_beta = size(beta); + size_t size_alpha_beta = max_size(alpha, beta); size_t max_size_seq_view = max_size(n, alpha, beta); operands_and_partials ops_partials(alpha, beta); - size_t len_ab = max_size(alpha, beta); - VectorBuilder lambda(len_ab); - for (size_t i = 0; i < len_ab; ++i) { - lambda[i] = value_of(alpha_vec[i]) / value_of(beta_vec[i]); - } - - VectorBuilder log1p_beta(size(beta)); - for (size_t i = 0; i < size(beta); ++i) { - log1p_beta[i] = log1p(value_of(beta_vec[i])); - } - - VectorBuilder log_beta_m_log1p_beta( - size(beta)); - for (size_t i = 0; i < size(beta); ++i) { - log_beta_m_log1p_beta[i] = log(value_of(beta_vec[i])) - log1p_beta[i]; - } - - VectorBuilder - alpha_times_log_beta_over_1p_beta(len_ab); - for (size_t i = 0; i < len_ab; ++i) { - alpha_times_log_beta_over_1p_beta[i] - = value_of(alpha_vec[i]) - * log(value_of(beta_vec[i]) / (1.0 + value_of(beta_vec[i]))); - } - VectorBuilder::value, T_partials_return, T_shape> - digamma_alpha(size(alpha)); + digamma_alpha(size_alpha); if (!is_constant_all::value) { - for (size_t i = 0; i < size(alpha); ++i) { + for (size_t i = 0; i < size_alpha; ++i) { digamma_alpha[i] = digamma(value_of(alpha_vec[i])); } } - VectorBuilder::value, T_partials_return, - T_inv_scale> - log_beta(size(beta)); - if (!is_constant_all::value) { - for (size_t i = 0; i < size(beta); ++i) { - log_beta[i] = log(value_of(beta_vec[i])); - } + VectorBuilder log_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_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(len_ab); - if (!is_constant_all::value) { - for (size_t i = 0; i < len_ab; ++i) { - lambda_m_alpha_over_1p_beta[i] - = lambda[i] - - (value_of(alpha_vec[i]) / (1.0 + value_of(beta_vec[i]))); + 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.0 + beta_dbl); } } @@ -123,7 +109,7 @@ return_type_t neg_binomial_lpmf(const T_n& n, n_vec[i] + value_of(alpha_vec[i]) - 1.0, n_vec[i]); } } - logp += alpha_times_log_beta_over_1p_beta[i] - n_vec[i] * log1p_beta[i]; + logp += alpha_log_beta_over_1p_beta[i] - n_vec[i] * log1p_beta[i]; if (!is_constant_all::value) { ops_partials.edge1_.partials_[i] @@ -137,6 +123,7 @@ return_type_t neg_binomial_lpmf(const T_n& n, } } } + return ops_partials.build(logp); } From dd0aee41b58f1effbb516f7dda3ad64463a00982 Mon Sep 17 00:00:00 2001 From: Marco Colombo Date: Tue, 4 Feb 2020 10:32:04 +0000 Subject: [PATCH 3/9] Add missing type to VectorBuilder call --- stan/math/prim/prob/neg_binomial_cdf.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/stan/math/prim/prob/neg_binomial_cdf.hpp b/stan/math/prim/prob/neg_binomial_cdf.hpp index 0a5d77d5c02..70d564c53ab 100644 --- a/stan/math/prim/prob/neg_binomial_cdf.hpp +++ b/stan/math/prim/prob/neg_binomial_cdf.hpp @@ -54,7 +54,8 @@ return_type_t neg_binomial_cdf(const T_n& n, VectorBuilder::value, T_partials_return, T_shape> digamma_alpha_vec(size_alpha); - VectorBuilder::value, T_partials_return, T_shape> + VectorBuilder::value, T_partials_return, T_n, + T_shape> digamma_sum_vec(size_n_alpha); if (!is_constant_all::value) { From 8ed2e371c0392808b6e12752ba9051f4e4239691 Mon Sep 17 00:00:00 2001 From: Marco Colombo Date: Tue, 4 Feb 2020 10:32:50 +0000 Subject: [PATCH 4/9] Use 1 instead of 1.0 --- stan/math/prim/prob/neg_binomial_2_lcdf.hpp | 2 +- stan/math/prim/prob/neg_binomial_2_log_lpmf.hpp | 2 +- stan/math/prim/prob/neg_binomial_lpmf.hpp | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/stan/math/prim/prob/neg_binomial_2_lcdf.hpp b/stan/math/prim/prob/neg_binomial_2_lcdf.hpp index f4f36ff9282..4de35ae9eff 100644 --- a/stan/math/prim/prob/neg_binomial_2_lcdf.hpp +++ b/stan/math/prim/prob/neg_binomial_2_lcdf.hpp @@ -48,7 +48,7 @@ return_type_t neg_binomial_2_lcdf( VectorBuilder, T_n> np1(size_n); for (size_t i = 0; i < size_n; i++) { - np1[i] = n_vec[i] + 1.0; + np1[i] = n_vec[i] + 1; } return beta_cdf_log(phi_mu.data(), phi, np1.data()); diff --git a/stan/math/prim/prob/neg_binomial_2_log_lpmf.hpp b/stan/math/prim/prob/neg_binomial_2_log_lpmf.hpp index 3989e1da857..9479973cbd1 100644 --- a/stan/math/prim/prob/neg_binomial_2_log_lpmf.hpp +++ b/stan/math/prim/prob/neg_binomial_2_log_lpmf.hpp @@ -114,7 +114,7 @@ return_type_t neg_binomial_2_log_lpmf( if (!is_constant_all::value) { ops_partials.edge1_.partials_[i] - += n_vec[i] - n_plus_phi[i] / (phi_val[i] / exp_eta[i] + 1.0); + += n_vec[i] - n_plus_phi[i] / (phi_val[i] / exp_eta[i] + 1); } if (!is_constant_all::value) { ops_partials.edge2_.partials_[i] diff --git a/stan/math/prim/prob/neg_binomial_lpmf.hpp b/stan/math/prim/prob/neg_binomial_lpmf.hpp index 81667bb195b..60f5e5c5be6 100644 --- a/stan/math/prim/prob/neg_binomial_lpmf.hpp +++ b/stan/math/prim/prob/neg_binomial_lpmf.hpp @@ -83,7 +83,7 @@ return_type_t neg_binomial_lpmf(const T_n& n, 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.0 + beta_dbl); + lambda_m_alpha_over_1p_beta[i] = lambda[i] - alpha_dbl / (1 + beta_dbl); } } From 50702ef7c7244b9fa0386a21032a16945840df54 Mon Sep 17 00:00:00 2001 From: Marco Colombo Date: Tue, 4 Feb 2020 11:00:59 +0000 Subject: [PATCH 5/9] Move mix test to the mix directory --- .../math/mix/prob/neg_binomial_2_log_test.cpp | 18 ++++++++++++++++++ .../math/rev/prob/neg_binomial_2_log_test.cpp | 19 ++----------------- 2 files changed, 20 insertions(+), 17 deletions(-) create mode 100644 test/unit/math/mix/prob/neg_binomial_2_log_test.cpp diff --git a/test/unit/math/mix/prob/neg_binomial_2_log_test.cpp b/test/unit/math/mix/prob/neg_binomial_2_log_test.cpp new file mode 100644 index 00000000000..a04071c0ccb --- /dev/null +++ b/test/unit/math/mix/prob/neg_binomial_2_log_test.cpp @@ -0,0 +1,18 @@ +#include +#include +#include +#include + +TEST(mathMixScalFun, neg_binomial_2_log_lpmf_derivatives) { + auto f1 = [](const auto& eta, const auto& phi) { + return stan::math::neg_binomial_2_log_lpmf(0, eta, phi); + }; + auto f2 = [](const auto& eta, const auto& phi) { + return stan::math::neg_binomial_2_log_lpmf(6, eta, phi); + }; + + stan::test::expect_ad(f1, -1.5, 4.1); + stan::test::expect_ad(f1, 2.0, 1.1); + stan::test::expect_ad(f2, -1.5, 4.1); + stan::test::expect_ad(f2, 2.0, 1.1); +} diff --git a/test/unit/math/rev/prob/neg_binomial_2_log_test.cpp b/test/unit/math/rev/prob/neg_binomial_2_log_test.cpp index ef48ba12592..ea2261fe995 100644 --- a/test/unit/math/rev/prob/neg_binomial_2_log_test.cpp +++ b/test/unit/math/rev/prob/neg_binomial_2_log_test.cpp @@ -1,11 +1,10 @@ #include -#include #include #include #include -#include #include #include +#include namespace neg_binomial_2_log_test_internal { struct TestValue { @@ -455,7 +454,7 @@ TEST(ProbDistributionsNegBinomial2Log, derivativesPrecomputed) { } } -TEST(ProbDistributionsNegBinomial2, derivativesComplexStep) { +TEST(ProbDistributionsNegBinomial2Log, derivativesComplexStep) { using boost::math::differentiation::complex_step_derivative; using stan::math::is_nan; using stan::math::neg_binomial_2_log_lpmf; @@ -548,17 +547,3 @@ TEST(ProbDistributionsNegBinomial2, derivativesComplexStep) { } } } - -TEST(mathMixScalFun, neg_binomial_2_log_lpmf_derivatives) { - auto f1 = [](const auto& eta, const auto& phi) { - return stan::math::neg_binomial_2_log_lpmf(0, eta, phi); - }; - auto f2 = [](const auto& eta, const auto& phi) { - return stan::math::neg_binomial_2_log_lpmf(6, eta, phi); - }; - - stan::test::expect_ad(f1, -1.5, 4.1); - stan::test::expect_ad(f1, 2.0, 1.1); - stan::test::expect_ad(f2, -1.5, 4.1); - stan::test::expect_ad(f2, 2.0, 1.1); -} From 5ca140052cad80000a172c823801ade6f30bb542 Mon Sep 17 00:00:00 2001 From: Marco Colombo Date: Tue, 4 Feb 2020 11:18:30 +0000 Subject: [PATCH 6/9] Add tests for derivatives when inputs have different size --- .../rev/prob/neg_binomial_2_ccdf_log_test.cpp | 32 +++++++++++++++ .../rev/prob/neg_binomial_2_cdf_log_test.cpp | 32 +++++++++++++++ .../math/rev/prob/neg_binomial_2_cdf_test.cpp | 32 +++++++++++++++ .../math/rev/prob/neg_binomial_2_log_test.cpp | 26 +++++++++++++ .../math/rev/prob/neg_binomial_2_test.cpp | 39 ++++++++++++++++--- .../rev/prob/neg_binomial_ccdf_log_test.cpp | 32 +++++++++++++++ .../rev/prob/neg_binomial_cdf_log_test.cpp | 32 +++++++++++++++ .../math/rev/prob/neg_binomial_cdf_test.cpp | 32 +++++++++++++++ test/unit/math/rev/prob/neg_binomial_test.cpp | 29 ++++++++++++++ 9 files changed, 280 insertions(+), 6 deletions(-) create mode 100644 test/unit/math/rev/prob/neg_binomial_2_ccdf_log_test.cpp create mode 100644 test/unit/math/rev/prob/neg_binomial_2_cdf_log_test.cpp create mode 100644 test/unit/math/rev/prob/neg_binomial_2_cdf_test.cpp create mode 100644 test/unit/math/rev/prob/neg_binomial_ccdf_log_test.cpp create mode 100644 test/unit/math/rev/prob/neg_binomial_cdf_log_test.cpp create mode 100644 test/unit/math/rev/prob/neg_binomial_cdf_test.cpp create mode 100644 test/unit/math/rev/prob/neg_binomial_test.cpp diff --git a/test/unit/math/rev/prob/neg_binomial_2_ccdf_log_test.cpp b/test/unit/math/rev/prob/neg_binomial_2_ccdf_log_test.cpp new file mode 100644 index 00000000000..e9c48c8992a --- /dev/null +++ b/test/unit/math/rev/prob/neg_binomial_2_ccdf_log_test.cpp @@ -0,0 +1,32 @@ +#include +#include +#include + +TEST(ProbDistributionsNegBinomial2, derivatives_lccdf) { + using stan::math::neg_binomial_2_lccdf; + using stan::math::var; + + std::vector N{1, 2, 3}; + double alpha_dbl = 8; + double beta_dbl = 1.5; + + var alpha(alpha_dbl); + var beta(beta_dbl); + + var val = neg_binomial_2_lccdf(N, alpha, beta); + std::vector x{alpha, beta}; + std::vector gradients; + val.grad(x, gradients); + + double epsilon = 1e-6; + + double grad_diff1 = (neg_binomial_2_lccdf(N, alpha_dbl + epsilon, beta_dbl) + - neg_binomial_2_lccdf(N, alpha_dbl - epsilon, beta_dbl)) + / (2 * epsilon); + EXPECT_FLOAT_EQ(grad_diff1, gradients[0]); + + double grad_diff2 = (neg_binomial_2_lccdf(N, alpha_dbl, beta_dbl + epsilon) + - neg_binomial_2_lccdf(N, alpha_dbl, beta_dbl - epsilon)) + / (2 * epsilon); + EXPECT_FLOAT_EQ(grad_diff2, gradients[1]); +} diff --git a/test/unit/math/rev/prob/neg_binomial_2_cdf_log_test.cpp b/test/unit/math/rev/prob/neg_binomial_2_cdf_log_test.cpp new file mode 100644 index 00000000000..a29476379ee --- /dev/null +++ b/test/unit/math/rev/prob/neg_binomial_2_cdf_log_test.cpp @@ -0,0 +1,32 @@ +#include +#include +#include + +TEST(ProbDistributionsNegBinomial2, derivatives_lcdf) { + using stan::math::neg_binomial_2_lcdf; + using stan::math::var; + + std::vector N{1, 2, 3}; + double alpha_dbl = 8; + double beta_dbl = 1.5; + + var alpha(alpha_dbl); + var beta(beta_dbl); + + var val = neg_binomial_2_lcdf(N, alpha, beta); + std::vector x{alpha, beta}; + std::vector gradients; + val.grad(x, gradients); + + double epsilon = 1e-6; + + double grad_diff1 = (neg_binomial_2_lcdf(N, alpha_dbl + epsilon, beta_dbl) + - neg_binomial_2_lcdf(N, alpha_dbl - epsilon, beta_dbl)) + / (2 * epsilon); + EXPECT_FLOAT_EQ(grad_diff1, gradients[0]); + + double grad_diff2 = (neg_binomial_2_lcdf(N, alpha_dbl, beta_dbl + epsilon) + - neg_binomial_2_lcdf(N, alpha_dbl, beta_dbl - epsilon)) + / (2 * epsilon); + EXPECT_FLOAT_EQ(grad_diff2, gradients[1]); +} diff --git a/test/unit/math/rev/prob/neg_binomial_2_cdf_test.cpp b/test/unit/math/rev/prob/neg_binomial_2_cdf_test.cpp new file mode 100644 index 00000000000..0d670b84674 --- /dev/null +++ b/test/unit/math/rev/prob/neg_binomial_2_cdf_test.cpp @@ -0,0 +1,32 @@ +#include +#include +#include + +TEST(ProbDistributionsNegBinomial2, derivatives_cdf) { + using stan::math::neg_binomial_2_cdf; + using stan::math::var; + + std::vector N{1, 2, 3}; + double alpha_dbl = 8; + double beta_dbl = 1.5; + + var alpha(alpha_dbl); + var beta(beta_dbl); + + var val = neg_binomial_2_cdf(N, alpha, beta); + std::vector x{alpha, beta}; + std::vector gradients; + val.grad(x, gradients); + + double epsilon = 1e-6; + + double grad_diff1 = (neg_binomial_2_cdf(N, alpha_dbl + epsilon, beta_dbl) + - neg_binomial_2_cdf(N, alpha_dbl - epsilon, beta_dbl)) + / (2 * epsilon); + EXPECT_FLOAT_EQ(grad_diff1, gradients[0]); + + double grad_diff2 = (neg_binomial_2_cdf(N, alpha_dbl, beta_dbl + epsilon) + - neg_binomial_2_cdf(N, alpha_dbl, beta_dbl - epsilon)) + / (2 * epsilon); + EXPECT_FLOAT_EQ(grad_diff2, gradients[1]); +} diff --git a/test/unit/math/rev/prob/neg_binomial_2_log_test.cpp b/test/unit/math/rev/prob/neg_binomial_2_log_test.cpp index ea2261fe995..c5aa187fb01 100644 --- a/test/unit/math/rev/prob/neg_binomial_2_log_test.cpp +++ b/test/unit/math/rev/prob/neg_binomial_2_log_test.cpp @@ -547,3 +547,29 @@ TEST(ProbDistributionsNegBinomial2Log, derivativesComplexStep) { } } } + +TEST(ProbDistributionsNegBinomial2Log, derivatives_diff_sizes) { + using stan::math::neg_binomial_2_log_lpmf; + using stan::math::var; + + int N = 100; + double mu_dbl = 1.5; + std::vector phi_dbl{2, 4, 6, 8}; + + var mu(mu_dbl); + std::vector phi; + for (double i : phi_dbl) { + phi.push_back(var(i)); + } + var val = neg_binomial_2_log_lpmf(N, mu, phi); + + std::vector x{mu}; + std::vector gradients; + val.grad(x, gradients); + + double eps = 1e-6; + double grad_diff = (neg_binomial_2_log_lpmf(N, mu_dbl + eps, phi_dbl) + - neg_binomial_2_log_lpmf(N, mu_dbl - eps, phi_dbl)) + / (2 * eps); + EXPECT_FLOAT_EQ(grad_diff, gradients[0]); +} diff --git a/test/unit/math/rev/prob/neg_binomial_2_test.cpp b/test/unit/math/rev/prob/neg_binomial_2_test.cpp index 7f5e806acda..aeb66efc2af 100644 --- a/test/unit/math/rev/prob/neg_binomial_2_test.cpp +++ b/test/unit/math/rev/prob/neg_binomial_2_test.cpp @@ -4,7 +4,7 @@ TEST(ProbDistributionsNegBinomial, derivatives) { using stan::math::is_nan; - using stan::math::neg_binomial_2_log; + using stan::math::neg_binomial_2_lpmf; using stan::math::var; int N = 100; @@ -14,7 +14,7 @@ TEST(ProbDistributionsNegBinomial, derivatives) { for (int k = 0; k < 20; ++k) { var mu(mu_dbl); var phi(phi_dbl); - var val = neg_binomial_2_log(N, mu, phi); + var val = neg_binomial_2_lpmf(N, mu, phi); std::vector x; x.push_back(mu); @@ -30,10 +30,11 @@ TEST(ProbDistributionsNegBinomial, derivatives) { std::vector finite_diffs; double eps = 1e-10; double inv2e = 0.5 / eps; - double dmu = neg_binomial_2_log(N, mu_dbl + eps, phi_dbl) - - neg_binomial_2_log(N, mu_dbl - eps, phi_dbl); - double dphi = neg_binomial_2_log(N, mu_dbl, phi_dbl + eps) - - neg_binomial_2_log(N, mu_dbl, phi_dbl - eps); + + double dmu = neg_binomial_2_lpmf(N, mu_dbl + eps, phi_dbl) + - neg_binomial_2_lpmf(N, mu_dbl - eps, phi_dbl); + double dphi = neg_binomial_2_lpmf(N, mu_dbl, phi_dbl + eps) + - neg_binomial_2_lpmf(N, mu_dbl, phi_dbl - eps); finite_diffs.push_back(dmu * inv2e); finite_diffs.push_back(dphi * inv2e); @@ -44,3 +45,29 @@ TEST(ProbDistributionsNegBinomial, derivatives) { phi_dbl *= 10; } } + +TEST(ProbDistributionsNegBinomial, derivatives_diff_sizes) { + using stan::math::neg_binomial_2_lpmf; + using stan::math::var; + + int N = 100; + double mu_dbl = 1.5; + std::vector phi_dbl{2, 4, 6, 8}; + + var mu(mu_dbl); + std::vector phi; + for (double i : phi_dbl) { + phi.push_back(var(i)); + } + var val = neg_binomial_2_lpmf(N, mu, phi); + + std::vector x{mu}; + std::vector gradients; + val.grad(x, gradients); + + double eps = 1e-6; + double grad_diff = (neg_binomial_2_lpmf(N, mu_dbl + eps, phi_dbl) + - neg_binomial_2_lpmf(N, mu_dbl - eps, phi_dbl)) + / (2 * eps); + EXPECT_FLOAT_EQ(grad_diff, gradients[0]); +} diff --git a/test/unit/math/rev/prob/neg_binomial_ccdf_log_test.cpp b/test/unit/math/rev/prob/neg_binomial_ccdf_log_test.cpp new file mode 100644 index 00000000000..d062f935a7d --- /dev/null +++ b/test/unit/math/rev/prob/neg_binomial_ccdf_log_test.cpp @@ -0,0 +1,32 @@ +#include +#include +#include + +TEST(ProbDistributionsNegBinomial, derivatives_lccdf) { + using stan::math::neg_binomial_lccdf; + using stan::math::var; + + std::vector N{1, 2, 3}; + double alpha_dbl = 8; + double beta_dbl = 1.5; + + var alpha(alpha_dbl); + var beta(beta_dbl); + + var val = neg_binomial_lccdf(N, alpha, beta); + std::vector x{alpha, beta}; + std::vector gradients; + val.grad(x, gradients); + + double epsilon = 1e-6; + + double grad_diff1 = (neg_binomial_lccdf(N, alpha_dbl + epsilon, beta_dbl) + - neg_binomial_lccdf(N, alpha_dbl - epsilon, beta_dbl)) + / (2 * epsilon); + EXPECT_FLOAT_EQ(grad_diff1, gradients[0]); + + double grad_diff2 = (neg_binomial_lccdf(N, alpha_dbl, beta_dbl + epsilon) + - neg_binomial_lccdf(N, alpha_dbl, beta_dbl - epsilon)) + / (2 * epsilon); + EXPECT_FLOAT_EQ(grad_diff2, gradients[1]); +} diff --git a/test/unit/math/rev/prob/neg_binomial_cdf_log_test.cpp b/test/unit/math/rev/prob/neg_binomial_cdf_log_test.cpp new file mode 100644 index 00000000000..6c5ddcfef1a --- /dev/null +++ b/test/unit/math/rev/prob/neg_binomial_cdf_log_test.cpp @@ -0,0 +1,32 @@ +#include +#include +#include + +TEST(ProbDistributionsNegBinomial, derivatives_lcdf) { + using stan::math::neg_binomial_lcdf; + using stan::math::var; + + std::vector N{1, 2, 3}; + double alpha_dbl = 8; + double beta_dbl = 1.5; + + var alpha(alpha_dbl); + var beta(beta_dbl); + + var val = neg_binomial_lcdf(N, alpha, beta); + std::vector x{alpha, beta}; + std::vector gradients; + val.grad(x, gradients); + + double epsilon = 1e-6; + + double grad_diff1 = (neg_binomial_lcdf(N, alpha_dbl + epsilon, beta_dbl) + - neg_binomial_lcdf(N, alpha_dbl - epsilon, beta_dbl)) + / (2 * epsilon); + EXPECT_FLOAT_EQ(grad_diff1, gradients[0]); + + double grad_diff2 = (neg_binomial_lcdf(N, alpha_dbl, beta_dbl + epsilon) + - neg_binomial_lcdf(N, alpha_dbl, beta_dbl - epsilon)) + / (2 * epsilon); + EXPECT_FLOAT_EQ(grad_diff2, gradients[1]); +} diff --git a/test/unit/math/rev/prob/neg_binomial_cdf_test.cpp b/test/unit/math/rev/prob/neg_binomial_cdf_test.cpp new file mode 100644 index 00000000000..a70c0283f33 --- /dev/null +++ b/test/unit/math/rev/prob/neg_binomial_cdf_test.cpp @@ -0,0 +1,32 @@ +#include +#include +#include + +TEST(ProbDistributionsNegBinomial, derivatives_cdf) { + using stan::math::neg_binomial_cdf; + using stan::math::var; + + std::vector N{1, 2, 3}; + double alpha_dbl = 8; + double beta_dbl = 1.5; + + var alpha(alpha_dbl); + var beta(beta_dbl); + + var val = neg_binomial_cdf(N, alpha, beta); + std::vector x{alpha, beta}; + std::vector gradients; + val.grad(x, gradients); + + double epsilon = 1e-6; + + double grad_diff1 = (neg_binomial_cdf(N, alpha_dbl + epsilon, beta_dbl) + - neg_binomial_cdf(N, alpha_dbl - epsilon, beta_dbl)) + / (2 * epsilon); + EXPECT_FLOAT_EQ(grad_diff1, gradients[0]); + + double grad_diff2 = (neg_binomial_cdf(N, alpha_dbl, beta_dbl + epsilon) + - neg_binomial_cdf(N, alpha_dbl, beta_dbl - epsilon)) + / (2 * epsilon); + EXPECT_FLOAT_EQ(grad_diff2, gradients[1]); +} diff --git a/test/unit/math/rev/prob/neg_binomial_test.cpp b/test/unit/math/rev/prob/neg_binomial_test.cpp new file mode 100644 index 00000000000..7dafe1ea725 --- /dev/null +++ b/test/unit/math/rev/prob/neg_binomial_test.cpp @@ -0,0 +1,29 @@ +#include +#include +#include + +TEST(ProbDistributionsNegBinomial, derivatives_diff_sizes) { + using stan::math::neg_binomial_lpmf; + using stan::math::var; + + int N = 100; + double mu_dbl = 1.5; + std::vector phi_dbl{2, 4, 6, 8}; + + var mu(mu_dbl); + std::vector phi; + for (double i : phi_dbl) { + phi.push_back(var(i)); + } + var val = neg_binomial_lpmf(N, mu, phi); + + std::vector x{mu}; + std::vector gradients; + val.grad(x, gradients); + + double eps = 1e-6; + double grad_diff = (neg_binomial_lpmf(N, mu_dbl + eps, phi_dbl) + - neg_binomial_lpmf(N, mu_dbl - eps, phi_dbl)) + / (2 * eps); + EXPECT_FLOAT_EQ(grad_diff, gradients[0]); +} From 7585fb246a1b04032de92e961318612d2a855e68 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Tue, 4 Feb 2020 11:34:22 +0000 Subject: [PATCH 7/9] [Jenkins] auto-formatting by clang-format version 5.0.0-3~16.04.1 (tags/RELEASE_500/final) --- test/unit/math/rev/prob/neg_binomial_2_log_test.cpp | 2 +- test/unit/math/rev/prob/neg_binomial_2_test.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/test/unit/math/rev/prob/neg_binomial_2_log_test.cpp b/test/unit/math/rev/prob/neg_binomial_2_log_test.cpp index c5aa187fb01..fede0d0850f 100644 --- a/test/unit/math/rev/prob/neg_binomial_2_log_test.cpp +++ b/test/unit/math/rev/prob/neg_binomial_2_log_test.cpp @@ -570,6 +570,6 @@ TEST(ProbDistributionsNegBinomial2Log, derivatives_diff_sizes) { double eps = 1e-6; double grad_diff = (neg_binomial_2_log_lpmf(N, mu_dbl + eps, phi_dbl) - neg_binomial_2_log_lpmf(N, mu_dbl - eps, phi_dbl)) - / (2 * eps); + / (2 * eps); EXPECT_FLOAT_EQ(grad_diff, gradients[0]); } diff --git a/test/unit/math/rev/prob/neg_binomial_2_test.cpp b/test/unit/math/rev/prob/neg_binomial_2_test.cpp index aeb66efc2af..4a216d7b331 100644 --- a/test/unit/math/rev/prob/neg_binomial_2_test.cpp +++ b/test/unit/math/rev/prob/neg_binomial_2_test.cpp @@ -68,6 +68,6 @@ TEST(ProbDistributionsNegBinomial, derivatives_diff_sizes) { double eps = 1e-6; double grad_diff = (neg_binomial_2_lpmf(N, mu_dbl + eps, phi_dbl) - neg_binomial_2_lpmf(N, mu_dbl - eps, phi_dbl)) - / (2 * eps); + / (2 * eps); EXPECT_FLOAT_EQ(grad_diff, gradients[0]); } From 3e384704cb699f646a5deac92d31c28db90a4888 Mon Sep 17 00:00:00 2001 From: Marco Colombo Date: Tue, 4 Feb 2020 12:41:31 +0000 Subject: [PATCH 8/9] Make test names unique --- test/unit/math/rev/prob/neg_binomial_2_test.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/unit/math/rev/prob/neg_binomial_2_test.cpp b/test/unit/math/rev/prob/neg_binomial_2_test.cpp index 4a216d7b331..3dbc8cedd91 100644 --- a/test/unit/math/rev/prob/neg_binomial_2_test.cpp +++ b/test/unit/math/rev/prob/neg_binomial_2_test.cpp @@ -2,7 +2,7 @@ #include #include -TEST(ProbDistributionsNegBinomial, derivatives) { +TEST(ProbDistributionsNegBinomial2, derivatives) { using stan::math::is_nan; using stan::math::neg_binomial_2_lpmf; using stan::math::var; @@ -46,7 +46,7 @@ TEST(ProbDistributionsNegBinomial, derivatives) { } } -TEST(ProbDistributionsNegBinomial, derivatives_diff_sizes) { +TEST(ProbDistributionsNegBinomial2, derivatives_diff_sizes) { using stan::math::neg_binomial_2_lpmf; using stan::math::var; From ec1031de41703d3ca29c7dae08a4f2e168f4787a Mon Sep 17 00:00:00 2001 From: Marco Colombo Date: Tue, 4 Feb 2020 16:27:08 +0000 Subject: [PATCH 9/9] Fix distribution test failure --- stan/math/prim/prob/neg_binomial_2_lcdf.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/prob/neg_binomial_2_lcdf.hpp b/stan/math/prim/prob/neg_binomial_2_lcdf.hpp index 4de35ae9eff..f4f36ff9282 100644 --- a/stan/math/prim/prob/neg_binomial_2_lcdf.hpp +++ b/stan/math/prim/prob/neg_binomial_2_lcdf.hpp @@ -48,7 +48,7 @@ return_type_t neg_binomial_2_lcdf( VectorBuilder, T_n> np1(size_n); for (size_t i = 0; i < size_n; i++) { - np1[i] = n_vec[i] + 1; + np1[i] = n_vec[i] + 1.0; } return beta_cdf_log(phi_mu.data(), phi, np1.data());