stan-dev · t4c1 · Oct 5, 2020 · Sep 15, 2020 · Sep 15, 2020 · Sep 15, 2020
diff --git a/stan/math/prim/prob/discrete_range_cdf.hpp b/stan/math/prim/prob/discrete_range_cdf.hpp
@@ -55,17 +55,16 @@ double discrete_range_cdf(const T_y& y, const T_lower& lower,
     if (y_dbl < lower_vec[n]) {
       return 0;
     }
-    if (y_dbl > upper_vec[n]) {
-      return 1;
-    }
   }
 
   double cdf(1.0);
   for (size_t n = 0; n < N; n++) {
     const double y_dbl = y_vec[n];
-    const double lower_dbl = lower_vec[n];
-    const double upper_dbl = upper_vec[n];
-    cdf *= (y_dbl - lower_dbl + 1) / (upper_dbl - lower_dbl + 1);
+    if (y_dbl <= upper_vec[n]) {
+      const double lower_dbl = lower_vec[n];
+      const double upper_dbl = upper_vec[n];
+      cdf *= (y_dbl - lower_dbl + 1) / (upper_dbl - lower_dbl + 1);
+    }
   }
   return cdf;
 }

diff --git a/stan/math/prim/prob/discrete_range_lccdf.hpp b/stan/math/prim/prob/discrete_range_lccdf.hpp
@@ -52,20 +52,19 @@ double discrete_range_lccdf(const T_y& y, const T_lower& lower,
 
   for (size_t n = 0; n < N; ++n) {
     const int y_dbl = y_vec[n];
-    if (y_dbl < lower_vec[n]) {
-      return 0;
-    }
-    if (y_dbl > upper_vec[n]) {
+    if (y_dbl >= upper_vec[n]) {
       return LOG_ZERO;
     }
   }
 
   double ccdf(0.0);
   for (size_t n = 0; n < N; n++) {
     const int y_dbl = y_vec[n];
-    const int lower_dbl = lower_vec[n];
-    const int upper_dbl = upper_vec[n];
-    ccdf += log(upper_dbl - y_dbl) - log(upper_dbl - lower_dbl + 1);
+    if (y_dbl >= lower_vec[n]) {
+      const int lower_dbl = lower_vec[n];
+      const int upper_dbl = upper_vec[n];
+      ccdf += log(upper_dbl - y_dbl) - log(upper_dbl - lower_dbl + 1);
+    }
   }
   return ccdf;
 }

diff --git a/stan/math/prim/prob/discrete_range_lcdf.hpp b/stan/math/prim/prob/discrete_range_lcdf.hpp
@@ -55,17 +55,16 @@ double discrete_range_lcdf(const T_y& y, const T_lower& lower,
     if (y_dbl < lower_vec[n]) {
       return LOG_ZERO;
     }
-    if (y_dbl > upper_vec[n]) {
-      return 0;
-    }
   }
 
   double cdf(0.0);
   for (size_t n = 0; n < N; n++) {
     const int y_dbl = y_vec[n];
-    const int lower_dbl = lower_vec[n];
-    const int upper_dbl = upper_vec[n];
-    cdf += log(y_dbl - lower_dbl + 1) - log(upper_dbl - lower_dbl + 1);
+    if (y_dbl < upper_vec[n]) {
+      const int lower_dbl = lower_vec[n];
+      const int upper_dbl = upper_vec[n];
+      cdf += log(y_dbl - lower_dbl + 1) - log(upper_dbl - lower_dbl + 1);
+    }
   }
   return cdf;
 }

diff --git a/stan/math/prim/prob/skew_normal_lpdf.hpp b/stan/math/prim/prob/skew_normal_lpdf.hpp
@@ -78,13 +78,17 @@ return_type_t<T_y, T_loc, T_scale, T_shape> skew_normal_lpdf(
     if (include_summand<propto, T_y, T_loc, T_scale>::value) {
       logp -= y_minus_mu_over_sigma * y_minus_mu_over_sigma / 2.0;
     }
-    logp += log(erfc(-alpha_dbl * y_minus_mu_over_sigma / SQRT_TWO));
+
+    T_partials_return log_erfc_alpha_z
+        = log(erfc(-alpha_dbl * y_minus_mu_over_sigma / SQRT_TWO));
+
+    logp += log_erfc_alpha_z;
 
     T_partials_return deriv_logerf
         = TWO_OVER_SQRT_PI
-          * exp(-alpha_dbl * y_minus_mu_over_sigma / SQRT_TWO * alpha_dbl
-                * y_minus_mu_over_sigma / SQRT_TWO)
-          / (1 + erf(alpha_dbl * y_minus_mu_over_sigma / SQRT_TWO));
+          * exp(-square(alpha_dbl * y_minus_mu_over_sigma / SQRT_TWO)
+                - log_erfc_alpha_z);
+
     if (!is_constant_all<T_y>::value) {
       ops_partials.edge1_.partials_[n]
           += -y_minus_mu_over_sigma / sigma_dbl

diff --git a/test/prob/discrete_range/discrete_range_ccdf_log_test.hpp b/test/prob/discrete_range/discrete_range_ccdf_log_test.hpp
@@ -56,7 +56,11 @@ class AgradCcdfLogDiscreteRange : public AgradCcdfLogTest {
   stan::return_type_t<T_y, T_lower, T_upper> ccdf_log_function(
       const T_y& y, const T_lower& lower, const T_upper& upper, const T3&,
       const T4&, const T5&) {
-    if (y < lower || y > upper) {
+    if (y < lower) {
+      return 0.0;
+    }
+
+    if (y >= upper) {
       return stan::math::LOG_ZERO;
     }
 

diff --git a/test/prob/discrete_range/discrete_range_cdf_log_test.hpp b/test/prob/discrete_range/discrete_range_cdf_log_test.hpp
@@ -54,10 +54,14 @@ class AgradCdfLogDiscreteRange : public AgradCdfLogTest {
   double cdf_log_function(const T_y& y, const T_lower& lower,
                           const T_upper& upper, const T3&, const T4&,
                           const T5&) {
-    if (y < lower || y > upper) {
+    if (y < lower) {
       return stan::math::LOG_ZERO;
     }
 
+    if (y > upper) {
+      return 0.0;
+    }
+
     return log((y - lower + 1.0) / (upper - lower + 1.0));
   }
 };
diff --git a/test/prob/discrete_range/discrete_range_cdf_test.hpp b/test/prob/discrete_range/discrete_range_cdf_test.hpp
@@ -52,10 +52,14 @@ class AgradCdfDiscreteRange : public AgradCdfTest {
             typename T5>
   double cdf_function(const T_y& y, const T_lower& lower, const T_upper& upper,
                       const T3&, const T4&, const T5&) {
-    if (y < lower || y > upper) {
+    if (y < lower) {
       return 0;
     }
 
+    if (y >= upper) {
+      return 1;
+    }
+
     return (y - lower + 1.0) / (upper - lower + 1.0);
   }
 };
diff --git a/test/prob/gumbel/gumbel_cdf_test.hpp b/test/prob/gumbel/gumbel_cdf_test.hpp
@@ -34,13 +34,11 @@ class AgradCdfGumbel : public AgradCdfTest {
         0.0006179789893310934986195216040530260548886143651007);  // expected
                                                                   // cdf
 
-    param[0] = -3.5;  // y
+    param[0] = -1.5;  // y
     param[1] = 1.9;   // mu
     param[2] = 7.2;   // beta
     parameters.push_back(param);
-    cdf.push_back(
-        0.1203922620798295861862650786832089422663975274508450);  // expected
-                                                                  // cdf
+    cdf.push_back(0.20118031381610754216);  // expected cdf
   }
 
   void invalid_values(vector<size_t>& index, vector<double>& value) {

diff --git a/test/prob/logistic/logistic_test.hpp b/test/prob/logistic/logistic_test.hpp
@@ -19,9 +19,9 @@ class AgradDistributionsLogistic : public AgradDistributionTest {
 
     param[0] = -1.0;  // y
     param[1] = 0.2;   // mu
-    param[2] = 0.25;  // sigma
+    param[2] = 0.5;   // sigma
     parameters.push_back(param);
-    log_prob.push_back(-3.430097773556644469295);  // expected log_prob
+    log_prob.push_back(-1.8805251237479538862);  // expected log_prob
   }
 
   void invalid_values(vector<size_t>& index, vector<double>& value) {

diff --git a/test/prob/neg_binomial/neg_binomial_test.hpp b/test/prob/neg_binomial/neg_binomial_test.hpp
@@ -24,11 +24,11 @@ class AgradDistributionsNegBinomial : public AgradDistributionTest {
     parameters.push_back(param);
     log_prob.push_back(-142.6147368129045105434);  // expected log_prob
 
-    param[0] = 13;
-    param[1] = 1e11;
-    param[2] = 1e10;
+    param[0] = 10;
+    param[1] = 1e6;
+    param[2] = 1e5;
     parameters.push_back(param);
-    log_prob.push_back(-2.6185576442208003);  // expected log_prob
+    log_prob.push_back(-2.0785666431081630812);  // expected log_prob
   }
 
   void invalid_values(vector<size_t>& index, vector<double>& value) {
@@ -75,8 +75,6 @@ class AgradDistributionsNegBinomial : public AgradDistributionTest {
     using stan::math::multiply_log;
     using std::log;
 
-    if (alpha > 1e10)
-      return -lgamma(n + 1.0) + multiply_log(n, alpha / beta) - alpha / beta;
     if (n != 0)
       return binomial_coefficient_log<
                  typename stan::scalar_type<T_shape>::type>(n + alpha - 1.0, n)