Fix weibull lcdf

stan/math/opencl/prim/weibull_cdf.hpp

@@ -67,8 +67,9 @@ return_type_t<T_y_cl, T_shape_cl, T_scale_cl> weibull_cdf(
   auto exp_n = exp(-pow_n);
   auto cdf_n = 1.0 - exp_n;
   auto cdf_expr = colwise_prod(cdf_n);
+  auto cdf_zero_expr = cdf_expr == 0;
 
-  auto rep_deriv = elt_multiply(exp_n, elt_divide(pow_n, cdf_n));
+  auto rep_deriv = elt_multiply(exp_n, pow_n);
   auto deriv_y_sigma = elt_multiply(rep_deriv, alpha_val);
   auto y_deriv_tmp = elt_divide(deriv_y_sigma, y_val);
   auto sigma_deriv_tmp = elt_divide(deriv_y_sigma, -sigma_val);
@@ -89,6 +90,10 @@ return_type_t<T_y_cl, T_shape_cl, T_scale_cl> weibull_cdf(
                     calc_if<!is_constant<T_shape_cl>::value>(alpha_deriv_tmp),
                     calc_if<!is_constant<T_scale_cl>::value>(sigma_deriv_tmp));
 
+  if (from_matrix_cl(cdf_zero_expr).any()) {
+    return 0;
+  }
+
   T_partials_return cdf = (from_matrix_cl(cdf_cl)).prod();
 
   auto alpha_deriv = alpha_deriv_cl * cdf;

stan/math/prim/prob/weibull_cdf.hpp

@@ -3,14 +3,15 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
 #include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
 #include <stan/math/prim/fun/exp.hpp>
 #include <stan/math/prim/fun/log.hpp>
-#include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/size.hpp>
 #include <stan/math/prim/fun/size_zero.hpp>
+#include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
 #include <stan/math/prim/functor/partials_propagator.hpp>
@@ -40,59 +41,69 @@ return_type_t<T_y, T_shape, T_scale> weibull_cdf(const T_y& y,
                                                  const T_shape& alpha,
                                                  const T_scale& sigma) {
   using T_partials_return = partials_return_t<T_y, T_shape, T_scale>;
-  using T_y_ref = ref_type_if_not_constant_t<T_y>;
-  using T_alpha_ref = ref_type_if_not_constant_t<T_shape>;
-  using T_sigma_ref = ref_type_if_not_constant_t<T_scale>;
+  using T_y_ref = ref_type_t<T_y>;
+  using T_alpha_ref = ref_type_t<T_shape>;
+  using T_sigma_ref = ref_type_t<T_scale>;
   using std::pow;
   static constexpr const char* function = "weibull_cdf";
 
   T_y_ref y_ref = y;
   T_alpha_ref alpha_ref = alpha;
   T_sigma_ref sigma_ref = sigma;
 
-  decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref));
-  decltype(auto) alpha_val = to_ref(as_value_column_array_or_scalar(alpha_ref));
-  decltype(auto) sigma_val = to_ref(as_value_column_array_or_scalar(sigma_ref));
-
-  check_nonnegative(function, "Random variable", y_val);
-  check_positive_finite(function, "Shape parameter", alpha_val);
-  check_positive_finite(function, "Scale parameter", sigma_val);
+  check_nonnegative(function, "Random variable", y_ref);
+  check_positive_finite(function, "Shape parameter", alpha_ref);
+  check_positive_finite(function, "Scale parameter", sigma_ref);
 
   if (size_zero(y, alpha, sigma)) {
     return 1.0;
   }
 
-  auto ops_partials = make_partials_propagator(y_ref, alpha_ref, sigma_ref);
+  T_partials_return cdf(1.0);
+  operands_and_partials<T_y_ref, T_alpha_ref, T_sigma_ref> ops_partials(
+      y_ref, alpha_ref, sigma_ref);
+
+  scalar_seq_view<T_y_ref> y_vec(y_ref);
+  scalar_seq_view<T_alpha_ref> alpha_vec(alpha_ref);
+  scalar_seq_view<T_sigma_ref> sigma_vec(sigma_ref);
+
+  size_t max_size_seq_view = max_size(y, alpha, sigma);
+  size_t size_y = stan::math::size(y);
+
+  // Explicit return for extreme values
+  // The gradients are technically ill-defined, but treated as zero
+  for (size_t i = 0; i < size_y; i++) {
+    if (y_vec[i] == 0) {
+      return ops_partials.build(0);
+    }
+  }
+
+  for (size_t i = 0; i < max_size_seq_view; i++) {
+    const T_partials_return pow_n
+        = pow(y_vec.val(i) / sigma_vec.val(i), alpha_vec.val(i));
+    const T_partials_return exp_n = exp(-pow_n);
+    const T_partials_return cdf_n = 1 - exp_n;
+    const T_partials_return rep_deriv = exp_n * pow_n;
 
-  constexpr bool any_derivs = !is_constant_all<T_y, T_shape, T_scale>::value;
-  const auto& pow_n = to_ref_if<any_derivs>(pow(y_val / sigma_val, alpha_val));
-  const auto& exp_n = to_ref_if<any_derivs>(exp(-pow_n));
-  const auto& cdf_n = to_ref_if<any_derivs>(1 - exp_n);
+    cdf *= cdf_n;
 
-  T_partials_return cdf = prod(cdf_n);
+    if (!is_constant_all<T_y, T_scale, T_shape>::value) {
+      const T_partials_return deriv_y_sigma = rep_deriv * alpha_vec.val(i);
 
-  if (any_derivs) {
-    const auto& rep_deriv = to_ref_if<(!is_constant_all<T_y, T_scale>::value
-                                       && !is_constant_all<T_shape>::value)>(
-        exp_n * pow_n * cdf / cdf_n);
-    if (!is_constant_all<T_y, T_scale>::value) {
-      const auto& deriv_y_sigma = to_ref_if<(
-          !is_constant_all<T_y>::value && !is_constant_all<T_scale>::value)>(
-          rep_deriv * alpha_val);
       if (!is_constant_all<T_y>::value) {
-        partials<0>(ops_partials) = deriv_y_sigma / y_val;
+        ops_partials.edge1_.partials_[i] += deriv_y_sigma / y_vec.val(i);
       }
       if (!is_constant_all<T_scale>::value) {
-        partials<2>(ops_partials) = -deriv_y_sigma / sigma_val;
+        ops_partials.edge3_.partials_[i] += -deriv_y_sigma / sigma_vec.val(i);
       }
     }
     if (!is_constant_all<T_shape>::value) {
-      partials<1>(ops_partials) = rep_deriv * log(y_val / sigma_val);
+      ops_partials.edge2_.partials_[i]
+          += rep_deriv * log(y_vec.val(i) / sigma_vec.val(i));
     }
   }
   return ops_partials.build(cdf);
 }
-
 }  // namespace math
 }  // namespace stan
 #endif

stan/math/prim/prob/weibull_lcdf.hpp

@@ -3,6 +3,7 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
 #include <stan/math/prim/fun/as_value_column_array_or_scalar.hpp>
@@ -11,6 +12,7 @@
 #include <stan/math/prim/fun/log1m.hpp>
 #include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
 #include <stan/math/prim/functor/partials_propagator.hpp>
@@ -40,59 +42,67 @@ return_type_t<T_y, T_shape, T_scale> weibull_lcdf(const T_y& y,
                                                   const T_shape& alpha,
                                                   const T_scale& sigma) {
   using T_partials_return = partials_return_t<T_y, T_shape, T_scale>;
-  using T_y_ref = ref_type_if_not_constant_t<T_y>;
-  using T_alpha_ref = ref_type_if_not_constant_t<T_shape>;
-  using T_sigma_ref = ref_type_if_not_constant_t<T_scale>;
+  using T_y_ref = ref_type_t<T_y>;
+  using T_alpha_ref = ref_type_t<T_shape>;
+  using T_sigma_ref = ref_type_t<T_scale>;
   using std::pow;
   static constexpr const char* function = "weibull_lcdf";
 
   T_y_ref y_ref = y;
   T_alpha_ref alpha_ref = alpha;
   T_sigma_ref sigma_ref = sigma;
 
-  decltype(auto) y_val = to_ref(as_value_column_array_or_scalar(y_ref));
-  decltype(auto) alpha_val = to_ref(as_value_column_array_or_scalar(alpha_ref));
-  decltype(auto) sigma_val = to_ref(as_value_column_array_or_scalar(sigma_ref));
-
-  check_nonnegative(function, "Random variable", y_val);
-  check_positive_finite(function, "Shape parameter", alpha_val);
-  check_positive_finite(function, "Scale parameter", sigma_val);
+  check_nonnegative(function, "Random variable", y_ref);
+  check_positive_finite(function, "Shape parameter", alpha_ref);
+  check_positive_finite(function, "Scale parameter", sigma_ref);
 
   if (size_zero(y, alpha, sigma)) {
     return 0.0;
   }
 
-  auto ops_partials = make_partials_propagator(y_ref, alpha_ref, sigma_ref);
+  T_partials_return cdf_log(0);
+  operands_and_partials<T_y_ref, T_alpha_ref, T_sigma_ref> ops_partials(
+      y_ref, alpha_ref, sigma_ref);
+
+  scalar_seq_view<T_y_ref> y_vec(y_ref);
+  scalar_seq_view<T_alpha_ref> alpha_vec(alpha_ref);
+  scalar_seq_view<T_sigma_ref> sigma_vec(sigma_ref);
+
+  size_t max_size_seq_view = max_size(y, alpha, sigma);
+  size_t size_y = stan::math::size(y);
+
+  for (size_t i = 0; i < size_y; i++) {
+    if (y_vec[i] == 0) {
+      return LOG_ZERO;
+    }
+  }
+
+  for (size_t i = 0; i < max_size_seq_view; i++) {
+    const T_partials_return pow_n
+        = pow(y_vec.val(i) / sigma_vec.val(i), alpha_vec.val(i));
+    const T_partials_return exp_n = exp(-pow_n);
+    const T_partials_return log1m_exp_n = log1m_exp(-pow_n);
+    const T_partials_return rep_deriv = exp_n * pow_n / exp(log1m_exp_n);
 
-  constexpr bool any_derivs = !is_constant_all<T_y, T_shape, T_scale>::value;
-  const auto& pow_n = to_ref_if<any_derivs>(pow(y_val / sigma_val, alpha_val));
-  const auto& exp_n = to_ref_if<any_derivs>(exp(-pow_n));
+    cdf_log += log1m_exp_n;
 
-  // TODO(Andrew) Further simplify derivatives and log1m_exp below
-  T_partials_return cdf_log = sum(log1m(exp_n));
+    if (!is_constant_all<T_y, T_scale, T_shape>::value) {
+      const T_partials_return deriv_y_sigma = rep_deriv * alpha_vec.val(i);
 
-  if (!is_constant_all<T_y, T_scale, T_shape>::value) {
-    const auto& rep_deriv = to_ref_if<(!is_constant_all<T_y, T_scale>::value
-                                       && !is_constant_all<T_shape>::value)>(
-        pow_n / (1.0 / exp_n - 1.0));
-    if (!is_constant_all<T_y, T_scale>::value) {
-      const auto& deriv_y_sigma = to_ref_if<(
-          !is_constant_all<T_y>::value && !is_constant_all<T_scale>::value)>(
-          rep_deriv * alpha_val);
       if (!is_constant_all<T_y>::value) {
-        partials<0>(ops_partials) = deriv_y_sigma / y_val;
+        ops_partials.edge1_.partials_[i] += deriv_y_sigma / y_vec.val(i);
       }
       if (!is_constant_all<T_scale>::value) {
-        partials<2>(ops_partials) = -deriv_y_sigma / sigma_val;
+        ops_partials.edge3_.partials_[i] += -deriv_y_sigma / sigma_vec.val(i);
       }
     }
     if (!is_constant_all<T_shape>::value) {
-      partials<1>(ops_partials) = rep_deriv * log(y_val / sigma_val);
+      ops_partials.edge2_.partials_[i]
+          += rep_deriv * log(y_vec.val(i) / sigma_vec.val(i));
     }
   }
   return ops_partials.build(cdf_log);
 }
-
 }  // namespace math
 }  // namespace stan
 #endif

test/unit/math/mix/prob/weibull_cdf_test.cpp

@@ -0,0 +1,17 @@
+#include <stan/math/mix.hpp>
+#include <test/unit/math/test_ad.hpp>
+
+TEST(mathMixScalFun, weibull_cdf) {
+  auto f = [](const double alpha, const double sigma) {
+    return
+        [=](const auto& y) { return stan::math::weibull_cdf(y, alpha, sigma); };
+  };
+
+  stan::test::expect_ad(f(0.1, 1.0), 0.1);
+  stan::test::expect_ad(f(1.0, 1.0), 0.001);
+  stan::test::expect_ad(f(0.1, 1.0), 1.25);
+  stan::test::expect_ad(f(1.0, 1.0), 0.15);
+  stan::test::expect_ad(f(0.5, 0.5), 0.25);
+  stan::test::expect_ad(f(1.0, 0.5), 2.0);
+  stan::test::expect_ad(f(2.5, 3.0), 1.5);
+}

test/unit/math/mix/prob/weibull_lcdf_test.cpp

@@ -0,0 +1,18 @@
+#include <stan/math/mix.hpp>
+#include <test/unit/math/test_ad.hpp>
+
+TEST(mathMixScalFun, weibull_lcdf) {
+  auto f = [](const double alpha, const double sigma) {
+    return [=](const auto& y) {
+      return stan::math::weibull_lcdf(y, alpha, sigma);
+    };
+  };
+
+  stan::test::expect_ad(f(0.1, 1.0), 0.1);
+  stan::test::expect_ad(f(1.0, 1.0), 0.001);
+  stan::test::expect_ad(f(0.1, 1.0), 1.25);
+  stan::test::expect_ad(f(1.0, 1.0), 0.15);
+  stan::test::expect_ad(f(0.5, 0.5), 0.25);
+  stan::test::expect_ad(f(1.0, 0.5), 2.0);
+  stan::test::expect_ad(f(2.5, 3.0), 1.5);
+}