Add probability tests to check vectorized/scalar lpmfs/lpdfs are the same

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;
 }

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;
 }

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;
 }

stan/math/prim/prob/frechet_lpdf.hpp

@@ -81,7 +81,7 @@ return_type_t<T_y, T_shape, T_scale> frechet_lpdf(const T_y& y,
                 T_partials_return, T_y>
       inv_y(size(y));
   for (size_t i = 0; i < stan::math::size(y); i++) {
-    inv_y[i] = 1.0 / value_of(y_vec[i]);
+    inv_y[i] = inv(value_of(y_vec[i]));
   }
 
   VectorBuilder<include_summand<propto, T_y, T_shape, T_scale>::value,
@@ -107,10 +107,9 @@ return_type_t<T_y, T_shape, T_scale> frechet_lpdf(const T_y& y,
     logp -= sigma_div_y_pow_alpha[n];
 
     if (!is_constant_all<T_y>::value) {
-      const T_partials_return inv_y_dbl = value_of(inv_y[n]);
       ops_partials.edge1_.partials_[n]
-          += -(alpha_dbl + 1.0) * inv_y_dbl
-             + alpha_dbl * sigma_div_y_pow_alpha[n] * inv_y_dbl;
+          += -(alpha_dbl + 1.0) * inv_y[n]
+             + alpha_dbl * sigma_div_y_pow_alpha[n] * inv_y[n];
     }
     if (!is_constant_all<T_shape>::value) {
       ops_partials.edge2_.partials_[n]

test/prob/binomial/binomial_cdf_test.hpp

@@ -52,10 +52,10 @@ class AgradCdfBinomial : public AgradCdfTest {
     using std::exp;
     using std::log;
 
-    stan::return_type_t<T_prob> cdf(1);
+    stan::return_type_t<T_prob> cdf(0);
 
     for (int i = 0; i <= n; i++) {
-      cdf *= binomial_coefficient<double>(N, i)
+      cdf += binomial_coefficient<double>(N, i)
              * exp(i * log(theta) + (N - i) * log(1 - theta));
     }
 

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 1.0;
+    }
+
+    if (y >= upper) {
       return stan::math::LOG_ZERO;
     }
 

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));
   }
 };

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);
   }
 };

test/prob/exponential/exponential_test.hpp

@@ -21,10 +21,10 @@ class AgradDistributionsExponential : public AgradDistributionTest {
     parameters.push_back(param);
     log_prob.push_back(-57.13902344686439249699);  // expected log_prob
 
-    param[0] = 1e-08;
+    param[0] = 1e-2;
     param[1] = 3.9;
     parameters.push_back(param);
-    log_prob.push_back(1.3609765141356007);
+    log_prob.push_back(1.3219765531356007);
   }
 
   void invalid_values(vector<size_t>& index, vector<double>& value) {

test/prob/frechet/frechet_test.hpp

@@ -75,7 +75,6 @@ class AgradDistributionsFrechet : public AgradDistributionTest {
       const T4&, const T5&) {
     using stan::math::include_summand;
     using stan::math::multiply_log;
-    using stan::math::value_of;
     using std::log;
     using std::pow;
 

test/prob/gamma/gamma_ccdf_log_test.hpp

@@ -32,11 +32,11 @@ class AgradCcdfLogGamma : public AgradCcdfLogTest {
     ccdf_log.push_back(
         std::log(1.0 - 0.6321205588285576659757));  // expected ccdf_log
 
-    param[0] = 16.0;  // y
-    param[1] = 3.0;   // alpha
-    param[2] = 3.0;   // beta
+    param[0] = 9.0;  // y
+    param[1] = 1.2;  // alpha
+    param[2] = 1.2;  // beta
     parameters.push_back(param);
-    ccdf_log.push_back(std::log(1.7116220633718619e-18));  // expected ccdf_log
+    ccdf_log.push_back(-10.221534317077268);  // expected ccdf_log
   }
 
   void invalid_values(vector<size_t>& index, vector<double>& value) {
@@ -90,10 +90,10 @@ class AgradCcdfLogGamma : public AgradCcdfLogTest {
   stan::return_type_t<T_y, T_shape, T_inv_scale> ccdf_log_function(
       const T_y& y, const T_shape& alpha, const T_inv_scale& beta, const T3&,
       const T4&, const T5&) {
-    using boost::math::gamma_q;
-    using stan::math::gamma_q;
-    using std::log;
+    using boost::math::gamma_p;
+    using stan::math::gamma_p;
+    using stan::math::log1m;
 
-    return log(gamma_q(alpha, beta * y));
+    return log1m(gamma_p(alpha, beta * y));
   }
 };

test/prob/gumbel/gumbel_test.hpp

@@ -88,6 +88,6 @@ class AgradDistributionGumbel : public AgradDistributionTest {
   stan::return_type_t<T_y, T_loc, T_scale> log_prob_function(
       const T_y& y, const T_loc& mu, const T_scale& beta, const T3&, const T4&,
       const T5&) {
-    return -log(beta) - (y - mu) / beta + exp((mu - y) / beta);
+    return -log(beta) - (y - mu) / beta - exp((mu - y) / beta);
   }
 };

test/prob/inv_chi_square/inv_chi_square_ccdf_log_test.hpp

@@ -56,6 +56,6 @@ class AgradCcdfLogInvChiSquare : public AgradCcdfLogTest {
                                                     const T5&) {
     using stan::math::gamma_p;
 
-    return log(1.0 - gamma_p(0.5 * nu, 0.5 / y));
+    return log(gamma_p(0.5 * nu, 0.5 / y));
   }
 };