stan-dev · bbbales2 · Jul 22, 2020 · Mar 5, 2020 · May 2, 2020 · May 2, 2020
diff --git a/stan/math/prim/prob/neg_binomial_lpmf.hpp b/stan/math/prim/prob/neg_binomial_lpmf.hpp
@@ -65,69 +65,46 @@ return_type_t<T_shape, T_inv_scale> neg_binomial_lpmf(const T_n& n,
     }
   }
 
-  VectorBuilder<true, T_partials_return, T_inv_scale> log_beta(size_beta);
+  VectorBuilder<true, T_partials_return, T_inv_scale> log1p_inv_beta(size_beta);
   VectorBuilder<true, T_partials_return, T_inv_scale> log1p_beta(size_beta);
-  VectorBuilder<true, T_partials_return, T_inv_scale> 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<true, T_partials_return, T_shape, T_inv_scale> lambda(
-      size_alpha_beta);
-  VectorBuilder<true, T_partials_return, T_shape, T_inv_scale>
-      alpha_log_beta_over_1p_beta(size_alpha_beta);
   VectorBuilder<!is_constant_all<T_inv_scale>::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<T_inv_scale>::value) {
-      lambda_m_alpha_over_1p_beta[i] = lambda[i] - alpha_dbl / (1 + beta_dbl);
+  if (!is_constant_all<T_inv_scale>::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<propto>::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<T_shape>::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<T_inv_scale>::value) {
-        ops_partials.edge2_.partials_[i]
-            += (lambda[i] - n_vec[i]) / value_of(beta_vec[i]);
+    if (include_summand<propto, T_shape>::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<propto, T_shape>::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<T_shape>::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<T_inv_scale>::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<T_shape>::value) {
+      ops_partials.edge1_.partials_[i] += digamma(alpha_dbl + n_vec[i])
+                                          - digamma_alpha[i]
+                                          - log1p_inv_beta[i];
+    }
+    if (!is_constant_all<T_inv_scale>::value) {
+      ops_partials.edge2_.partials_[i]
+          += 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
@@ -71,7 +71,7 @@ return_type_t<T_y, T_loc, T_scale, T_shape> pareto_type_2_lpdf(
       inv_alpha(size(alpha));
   if (!is_constant_all<T_shape>::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]));
     }
   }
 
@@ -85,7 +85,7 @@ return_type_t<T_y, T_loc, T_scale, T_shape> 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<propto, T_y, T_scale, T_shape>::value
+        = include_summand<propto, T_y, T_loc, T_scale, T_shape>::value
               ? log1p((y_dbl - mu_dbl) / lambda_dbl)
               : 0;
 
@@ -95,7 +95,8 @@ return_type_t<T_y, T_loc, T_scale, T_shape> pareto_type_2_lpdf(
     if (include_summand<propto, T_scale>::value) {
       logp -= log_lambda[n];
     }
-    if (include_summand<propto, T_y, T_scale, T_shape>::value) {
+
+    if (include_summand<propto, T_y, T_loc, T_scale, T_shape>::value) {
       logp -= (alpha_dbl + 1.0) * log1p_scaled_diff;
     }
 

diff --git 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<size_t>& index, vector<double>& value) {

diff --git a/test/prob/test_fixture_distr.hpp b/test/prob/test_fixture_distr.hpp
@@ -13,6 +13,7 @@ using stan::is_vector;
 using stan::scalar_type;
 using stan::math::fvar;
 using stan::math::value_of;
+using stan::math::value_of_rec;
 using stan::math::var;
 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_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

diff --git a/test/prob/utility.hpp b/test/prob/utility.hpp
@@ -596,7 +596,7 @@ typename scalar_type<T>::type select_var_param(
     const vector<vector<double>>& parameters, const size_t n, const size_t p) {
   typename scalar_type<T>::type param(0);
   if (p < parameters[0].size()) {
-    if (is_vector<T>::value && !is_constant_all<T>::value)
+    if (!is_constant_all<T>::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
@@ -0,0 +1,19 @@
+#include <stan/math/rev.hpp>
+#include <test/unit/math/test_ad.hpp>
+#include <gtest/gtest.h>
+#include <boost/math/differentiation/finite_difference.hpp>
+
+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);
+    };
+  };
+
+  stan::test::expect_ad(f(0), 1.5, 4.1);
+  stan::test::expect_ad(f(0), std::vector<double>({1.2, 2.0}),
+                        std::vector<double>({1.1, 1.2}));
+  stan::test::expect_ad(f(6), 1.5, 4.1);
+  stan::test::expect_ad(f(6), std::vector<double>({1.7, 2.0}),
+                        std::vector<double>({1.1, 2.3}));
+}
diff --git 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,17 @@ 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);
+
   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<true>(y, alpha, beta)),

diff --git 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<double>(1.0 - p);
 }
 
 TEST(ProbDistributionsNegBinomial, chiSquareGoodnessFitTest) {