From ae6ca1fec59fae41734bb964de7969b0d9635670 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sun, 24 Mar 2024 14:02:20 +0100 Subject: [PATCH 1/7] Numerically stable weibull_cdf --- stan/math/prim/prob/weibull_cdf.hpp | 30 +++++++++----------- test/unit/math/mix/prob/weibull_cdf_test.cpp | 10 +++++++ 2 files changed, 24 insertions(+), 16 deletions(-) create mode 100644 test/unit/math/mix/prob/weibull_cdf_test.cpp diff --git a/stan/math/prim/prob/weibull_cdf.hpp b/stan/math/prim/prob/weibull_cdf.hpp index 11bd26ec3b1..aa12f36405e 100644 --- a/stan/math/prim/prob/weibull_cdf.hpp +++ b/stan/math/prim/prob/weibull_cdf.hpp @@ -62,35 +62,33 @@ return_type_t weibull_cdf(const T_y& y, return 1.0; } - auto ops_partials = make_partials_propagator(y_ref, alpha_ref, sigma_ref); - constexpr bool any_derivs = !is_constant_all::value; - const auto& pow_n = to_ref_if(pow(y_val / sigma_val, alpha_val)); - const auto& exp_n = to_ref_if(exp(-pow_n)); - const auto& cdf_n = to_ref_if(1 - exp_n); + const auto& log_y = to_ref_if(log(y_val)); + const auto& log_sigma = to_ref_if(log(sigma_val)); + const auto& log_pow_n = to_ref_if( + alpha_val * log_y - alpha_val * log_sigma); + const auto& pow_n = to_ref_if(exp(log_pow_n)); + const auto& log_cdf_n = to_ref_if(log1m_exp(-pow_n)); - T_partials_return cdf = prod(cdf_n); + T_partials_return log_cdf = sum(log_cdf_n); + auto ops_partials = make_partials_propagator(y_ref, alpha_ref, sigma_ref); if (any_derivs) { - const auto& rep_deriv = to_ref_if<(!is_constant_all::value - && !is_constant_all::value)>( - exp_n * pow_n * cdf / cdf_n); + const auto& log_rep_deriv = to_ref(log_pow_n + log_cdf - log_cdf_n - pow_n); if (!is_constant_all::value) { - const auto& deriv_y_sigma = to_ref_if<( - !is_constant_all::value && !is_constant_all::value)>( - rep_deriv * alpha_val); + const auto& log_deriv_y_sigma = to_ref(log_rep_deriv + log(alpha_val)); if (!is_constant_all::value) { - partials<0>(ops_partials) = deriv_y_sigma / y_val; + partials<0>(ops_partials) = exp(log_deriv_y_sigma - log_y); } if (!is_constant_all::value) { - partials<2>(ops_partials) = -deriv_y_sigma / sigma_val; + partials<2>(ops_partials) = -exp(log_deriv_y_sigma - log_sigma); } } if (!is_constant_all::value) { - partials<1>(ops_partials) = rep_deriv * log(y_val / sigma_val); + partials<1>(ops_partials) = exp(log_rep_deriv) * (log_y - log_sigma); } } - return ops_partials.build(cdf); + return ops_partials.build(exp(log_cdf)); } } // namespace math diff --git a/test/unit/math/mix/prob/weibull_cdf_test.cpp b/test/unit/math/mix/prob/weibull_cdf_test.cpp new file mode 100644 index 00000000000..98eb80729a7 --- /dev/null +++ b/test/unit/math/mix/prob/weibull_cdf_test.cpp @@ -0,0 +1,10 @@ +#include +#include + +TEST(mathMixScalFun, weibull_cdf) { + auto f = [](const auto& y, const auto& alpha, const auto& sigma) { + return stan::math::weibull_cdf(y, alpha, sigma); + }; + + stan::test::expect_ad(f, 10, 10, 10); +} From 54257d9f9982071546b3215e79d2feb6372b26d7 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sun, 24 Mar 2024 14:51:16 +0100 Subject: [PATCH 2/7] Add y=0 handling, mix tests --- stan/math/prim/prob/weibull_cdf.hpp | 6 +++++- test/unit/math/mix/prob/weibull_cdf_test.cpp | 15 +++++++++++++-- 2 files changed, 18 insertions(+), 3 deletions(-) diff --git a/stan/math/prim/prob/weibull_cdf.hpp b/stan/math/prim/prob/weibull_cdf.hpp index aa12f36405e..856eae836cc 100644 --- a/stan/math/prim/prob/weibull_cdf.hpp +++ b/stan/math/prim/prob/weibull_cdf.hpp @@ -62,6 +62,11 @@ return_type_t weibull_cdf(const T_y& y, return 1.0; } + auto ops_partials = make_partials_propagator(y_ref, alpha_ref, sigma_ref); + if (any(value_of_rec(y_val) == 0)) { + return ops_partials.build(0.0); + } + constexpr bool any_derivs = !is_constant_all::value; const auto& log_y = to_ref_if(log(y_val)); const auto& log_sigma = to_ref_if(log(sigma_val)); @@ -72,7 +77,6 @@ return_type_t weibull_cdf(const T_y& y, T_partials_return log_cdf = sum(log_cdf_n); - auto ops_partials = make_partials_propagator(y_ref, alpha_ref, sigma_ref); if (any_derivs) { const auto& log_rep_deriv = to_ref(log_pow_n + log_cdf - log_cdf_n - pow_n); if (!is_constant_all::value) { diff --git a/test/unit/math/mix/prob/weibull_cdf_test.cpp b/test/unit/math/mix/prob/weibull_cdf_test.cpp index 98eb80729a7..0494c83532d 100644 --- a/test/unit/math/mix/prob/weibull_cdf_test.cpp +++ b/test/unit/math/mix/prob/weibull_cdf_test.cpp @@ -2,9 +2,20 @@ #include TEST(mathMixScalFun, weibull_cdf) { + // Inputs are tested on the log (i.e., unconstrained) scale so that the + // finite-diffs don't result in invalid inputs. auto f = [](const auto& y, const auto& alpha, const auto& sigma) { - return stan::math::weibull_cdf(y, alpha, sigma); + using stan::math::lb_constrain; + using stan::math::positive_constrain; + + return stan::math::weibull_cdf(lb_constrain(y, 0.0), + positive_constrain(alpha), + positive_constrain(sigma)); }; - stan::test::expect_ad(f, 10, 10, 10); + using stan::math::log; + + stan::test::expect_ad(f, log(1.2), log(4), log(20)); + stan::test::expect_ad(f, 0.0, 0.0, 0.0); + stan::test::expect_ad(f, stan::math::NEGATIVE_INFTY, 0.0, 0.0); // y = 0.0 } From e87b32ee7afa453b1c4cd811aa1d25bde4ad964f Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sun, 24 Mar 2024 15:17:40 +0100 Subject: [PATCH 3/7] Update lcdf and vectorise mix tests --- stan/math/prim/prob/weibull_cdf.hpp | 1 - stan/math/prim/prob/weibull_lcdf.hpp | 31 ++++++++-------- test/unit/math/mix/prob/weibull_cdf_test.cpp | 22 +++++++++-- test/unit/math/mix/prob/weibull_lcdf_test.cpp | 37 +++++++++++++++++++ 4 files changed, 71 insertions(+), 20 deletions(-) create mode 100644 test/unit/math/mix/prob/weibull_lcdf_test.cpp diff --git a/stan/math/prim/prob/weibull_cdf.hpp b/stan/math/prim/prob/weibull_cdf.hpp index 856eae836cc..9c9c0833475 100644 --- a/stan/math/prim/prob/weibull_cdf.hpp +++ b/stan/math/prim/prob/weibull_cdf.hpp @@ -43,7 +43,6 @@ return_type_t weibull_cdf(const T_y& y, using T_y_ref = ref_type_if_not_constant_t; using T_alpha_ref = ref_type_if_not_constant_t; using T_sigma_ref = ref_type_if_not_constant_t; - using std::pow; static constexpr const char* function = "weibull_cdf"; T_y_ref y_ref = y; diff --git a/stan/math/prim/prob/weibull_lcdf.hpp b/stan/math/prim/prob/weibull_lcdf.hpp index 98d927f602e..5dac327e365 100644 --- a/stan/math/prim/prob/weibull_lcdf.hpp +++ b/stan/math/prim/prob/weibull_lcdf.hpp @@ -43,7 +43,6 @@ return_type_t weibull_lcdf(const T_y& y, using T_y_ref = ref_type_if_not_constant_t; using T_alpha_ref = ref_type_if_not_constant_t; using T_sigma_ref = ref_type_if_not_constant_t; - using std::pow; static constexpr const char* function = "weibull_lcdf"; T_y_ref y_ref = y; @@ -63,34 +62,34 @@ return_type_t weibull_lcdf(const T_y& y, } auto ops_partials = make_partials_propagator(y_ref, alpha_ref, sigma_ref); + if (any(value_of_rec(y_val) == 0)) { + return ops_partials.build(stan::math::NEGATIVE_INFTY); + } constexpr bool any_derivs = !is_constant_all::value; - const auto& pow_n = to_ref_if(pow(y_val / sigma_val, alpha_val)); - const auto& exp_n = to_ref_if(exp(-pow_n)); + const auto& log_y = to_ref_if(log(y_val)); + const auto& log_sigma = to_ref_if(log(sigma_val)); + const auto& log_pow_n = to_ref_if( + alpha_val * log_y - alpha_val * log_sigma); + const auto& pow_n = to_ref_if(exp(log_pow_n)); - // TODO(Andrew) Further simplify derivatives and log1m_exp below - T_partials_return cdf_log = sum(log1m(exp_n)); + if (any_derivs) { + const auto& log_rep_deriv = to_ref(log_pow_n - log_diff_exp(pow_n, 0.0)); - if (!is_constant_all::value) { - const auto& rep_deriv = to_ref_if<(!is_constant_all::value - && !is_constant_all::value)>( - pow_n / (1.0 / exp_n - 1.0)); if (!is_constant_all::value) { - const auto& deriv_y_sigma = to_ref_if<( - !is_constant_all::value && !is_constant_all::value)>( - rep_deriv * alpha_val); + const auto& log_deriv_y_sigma = to_ref(log_rep_deriv + log(alpha_val)); if (!is_constant_all::value) { - partials<0>(ops_partials) = deriv_y_sigma / y_val; + partials<0>(ops_partials) = exp(log_deriv_y_sigma - log_y); } if (!is_constant_all::value) { - partials<2>(ops_partials) = -deriv_y_sigma / sigma_val; + partials<2>(ops_partials) = -exp(log_deriv_y_sigma - log_sigma); } } if (!is_constant_all::value) { - partials<1>(ops_partials) = rep_deriv * log(y_val / sigma_val); + partials<1>(ops_partials) = exp(log_rep_deriv) * (log_y - log_sigma); } } - return ops_partials.build(cdf_log); + return ops_partials.build(sum(log1m_exp(-pow_n))); } } // namespace math diff --git a/test/unit/math/mix/prob/weibull_cdf_test.cpp b/test/unit/math/mix/prob/weibull_cdf_test.cpp index 0494c83532d..988bff4637e 100644 --- a/test/unit/math/mix/prob/weibull_cdf_test.cpp +++ b/test/unit/math/mix/prob/weibull_cdf_test.cpp @@ -15,7 +15,23 @@ TEST(mathMixScalFun, weibull_cdf) { using stan::math::log; - stan::test::expect_ad(f, log(1.2), log(4), log(20)); - stan::test::expect_ad(f, 0.0, 0.0, 0.0); - stan::test::expect_ad(f, stan::math::NEGATIVE_INFTY, 0.0, 0.0); // y = 0.0 + Eigen::VectorXd y(3); + y << stan::math::NEGATIVE_INFTY, 1.2, 0.0; + + Eigen::VectorXd alpha(3); + alpha << 2.0, 3.0, 4.0; + + Eigen::VectorXd sigma(3); + sigma << 5.0, 6.0, 7.0; + + stan::test::expect_ad(f, y, alpha, sigma); + stan::test::expect_ad(f, y[0], alpha, sigma); + stan::test::expect_ad(f, y, alpha[0], sigma); + stan::test::expect_ad(f, y, alpha, sigma[0]); + + stan::test::expect_ad(f, y[0], alpha[0], sigma); + stan::test::expect_ad(f, y[0], alpha, sigma[0]); + stan::test::expect_ad(f, y, alpha[0], sigma[0]); + + stan::test::expect_ad(f, y[0], alpha[0], sigma[0]); } diff --git a/test/unit/math/mix/prob/weibull_lcdf_test.cpp b/test/unit/math/mix/prob/weibull_lcdf_test.cpp new file mode 100644 index 00000000000..b159ab08f03 --- /dev/null +++ b/test/unit/math/mix/prob/weibull_lcdf_test.cpp @@ -0,0 +1,37 @@ +#include +#include + +TEST(mathMixScalFun, weibull_lcdf) { + // Inputs are tested on the log (i.e., unconstrained) scale so that the + // finite-diffs don't result in invalid inputs. + auto f = [](const auto& y, const auto& alpha, const auto& sigma) { + using stan::math::lb_constrain; + using stan::math::positive_constrain; + + return stan::math::weibull_lcdf(lb_constrain(y, 0.0), + positive_constrain(alpha), + positive_constrain(sigma)); + }; + + using stan::math::log; + + Eigen::VectorXd y(3); + y << stan::math::NEGATIVE_INFTY, 1.2, 0.0; + + Eigen::VectorXd alpha(3); + alpha << 2.0, 3.0, 4.0; + + Eigen::VectorXd sigma(3); + sigma << 5.0, 6.0, 7.0; + + stan::test::expect_ad(f, y, alpha, sigma); + stan::test::expect_ad(f, y[0], alpha, sigma); + stan::test::expect_ad(f, y, alpha[0], sigma); + stan::test::expect_ad(f, y, alpha, sigma[0]); + + stan::test::expect_ad(f, y[0], alpha[0], sigma); + stan::test::expect_ad(f, y[0], alpha, sigma[0]); + stan::test::expect_ad(f, y, alpha[0], sigma[0]); + + stan::test::expect_ad(f, y[0], alpha[0], sigma[0]); +} From 5cef0e83cb58d52283c775064db4aadb65a17ca9 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sun, 24 Mar 2024 15:20:22 +0100 Subject: [PATCH 4/7] Unused typedef --- stan/math/prim/prob/weibull_lcdf.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/stan/math/prim/prob/weibull_lcdf.hpp b/stan/math/prim/prob/weibull_lcdf.hpp index 5dac327e365..c6509d29d7b 100644 --- a/stan/math/prim/prob/weibull_lcdf.hpp +++ b/stan/math/prim/prob/weibull_lcdf.hpp @@ -39,7 +39,6 @@ template weibull_lcdf(const T_y& y, const T_shape& alpha, const T_scale& sigma) { - using T_partials_return = partials_return_t; using T_y_ref = ref_type_if_not_constant_t; using T_alpha_ref = ref_type_if_not_constant_t; using T_sigma_ref = ref_type_if_not_constant_t; From 1c0f6be03688a3561bacf5039018ce3eb6542b5c Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sun, 24 Mar 2024 15:25:37 +0100 Subject: [PATCH 5/7] Update comments --- test/unit/math/mix/prob/weibull_cdf_test.cpp | 2 +- test/unit/math/mix/prob/weibull_lcdf_test.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/test/unit/math/mix/prob/weibull_cdf_test.cpp b/test/unit/math/mix/prob/weibull_cdf_test.cpp index 988bff4637e..f3612ea2046 100644 --- a/test/unit/math/mix/prob/weibull_cdf_test.cpp +++ b/test/unit/math/mix/prob/weibull_cdf_test.cpp @@ -16,7 +16,7 @@ TEST(mathMixScalFun, weibull_cdf) { using stan::math::log; Eigen::VectorXd y(3); - y << stan::math::NEGATIVE_INFTY, 1.2, 0.0; + y << stan::math::NEGATIVE_INFTY, 1.2, 0.0; // lb_constrain(y[0], 0.0) = 0.0 Eigen::VectorXd alpha(3); alpha << 2.0, 3.0, 4.0; diff --git a/test/unit/math/mix/prob/weibull_lcdf_test.cpp b/test/unit/math/mix/prob/weibull_lcdf_test.cpp index b159ab08f03..214f14cb84b 100644 --- a/test/unit/math/mix/prob/weibull_lcdf_test.cpp +++ b/test/unit/math/mix/prob/weibull_lcdf_test.cpp @@ -16,7 +16,7 @@ TEST(mathMixScalFun, weibull_lcdf) { using stan::math::log; Eigen::VectorXd y(3); - y << stan::math::NEGATIVE_INFTY, 1.2, 0.0; + y << stan::math::NEGATIVE_INFTY, 1.2, 0.0; // lb_constrain(y[0], 0.0) = 0.0 Eigen::VectorXd alpha(3); alpha << 2.0, 3.0, 4.0; From ffdc348be4135cf2553f824e205cc4d67f4890bf Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 24 Mar 2024 10:31:38 -0400 Subject: [PATCH 6/7] [Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1 --- stan/math/prim/prob/weibull_cdf.hpp | 4 ++-- stan/math/prim/prob/weibull_lcdf.hpp | 4 ++-- test/unit/math/mix/prob/weibull_cdf_test.cpp | 6 +++--- test/unit/math/mix/prob/weibull_lcdf_test.cpp | 2 +- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/stan/math/prim/prob/weibull_cdf.hpp b/stan/math/prim/prob/weibull_cdf.hpp index 9c9c0833475..d6b87b23588 100644 --- a/stan/math/prim/prob/weibull_cdf.hpp +++ b/stan/math/prim/prob/weibull_cdf.hpp @@ -69,8 +69,8 @@ return_type_t weibull_cdf(const T_y& y, constexpr bool any_derivs = !is_constant_all::value; const auto& log_y = to_ref_if(log(y_val)); const auto& log_sigma = to_ref_if(log(sigma_val)); - const auto& log_pow_n = to_ref_if( - alpha_val * log_y - alpha_val * log_sigma); + const auto& log_pow_n + = to_ref_if(alpha_val * log_y - alpha_val * log_sigma); const auto& pow_n = to_ref_if(exp(log_pow_n)); const auto& log_cdf_n = to_ref_if(log1m_exp(-pow_n)); diff --git a/stan/math/prim/prob/weibull_lcdf.hpp b/stan/math/prim/prob/weibull_lcdf.hpp index c6509d29d7b..8152b759c60 100644 --- a/stan/math/prim/prob/weibull_lcdf.hpp +++ b/stan/math/prim/prob/weibull_lcdf.hpp @@ -68,8 +68,8 @@ return_type_t weibull_lcdf(const T_y& y, constexpr bool any_derivs = !is_constant_all::value; const auto& log_y = to_ref_if(log(y_val)); const auto& log_sigma = to_ref_if(log(sigma_val)); - const auto& log_pow_n = to_ref_if( - alpha_val * log_y - alpha_val * log_sigma); + const auto& log_pow_n + = to_ref_if(alpha_val * log_y - alpha_val * log_sigma); const auto& pow_n = to_ref_if(exp(log_pow_n)); if (any_derivs) { diff --git a/test/unit/math/mix/prob/weibull_cdf_test.cpp b/test/unit/math/mix/prob/weibull_cdf_test.cpp index f3612ea2046..9b9d3fcd08f 100644 --- a/test/unit/math/mix/prob/weibull_cdf_test.cpp +++ b/test/unit/math/mix/prob/weibull_cdf_test.cpp @@ -9,14 +9,14 @@ TEST(mathMixScalFun, weibull_cdf) { using stan::math::positive_constrain; return stan::math::weibull_cdf(lb_constrain(y, 0.0), - positive_constrain(alpha), - positive_constrain(sigma)); + positive_constrain(alpha), + positive_constrain(sigma)); }; using stan::math::log; Eigen::VectorXd y(3); - y << stan::math::NEGATIVE_INFTY, 1.2, 0.0; // lb_constrain(y[0], 0.0) = 0.0 + y << stan::math::NEGATIVE_INFTY, 1.2, 0.0; // lb_constrain(y[0], 0.0) = 0.0 Eigen::VectorXd alpha(3); alpha << 2.0, 3.0, 4.0; diff --git a/test/unit/math/mix/prob/weibull_lcdf_test.cpp b/test/unit/math/mix/prob/weibull_lcdf_test.cpp index 214f14cb84b..80e9b7f6a7b 100644 --- a/test/unit/math/mix/prob/weibull_lcdf_test.cpp +++ b/test/unit/math/mix/prob/weibull_lcdf_test.cpp @@ -16,7 +16,7 @@ TEST(mathMixScalFun, weibull_lcdf) { using stan::math::log; Eigen::VectorXd y(3); - y << stan::math::NEGATIVE_INFTY, 1.2, 0.0; // lb_constrain(y[0], 0.0) = 0.0 + y << stan::math::NEGATIVE_INFTY, 1.2, 0.0; // lb_constrain(y[0], 0.0) = 0.0 Eigen::VectorXd alpha(3); alpha << 2.0, 3.0, 4.0; From f627b3750f4dafea1d1d341e05d112cb39fd5428 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Tue, 26 Mar 2024 17:12:37 +0200 Subject: [PATCH 7/7] Re-use computation --- stan/math/prim/prob/weibull_cdf.hpp | 6 +++--- stan/math/prim/prob/weibull_lcdf.hpp | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/stan/math/prim/prob/weibull_cdf.hpp b/stan/math/prim/prob/weibull_cdf.hpp index d6b87b23588..dd40454401d 100644 --- a/stan/math/prim/prob/weibull_cdf.hpp +++ b/stan/math/prim/prob/weibull_cdf.hpp @@ -69,8 +69,8 @@ return_type_t weibull_cdf(const T_y& y, constexpr bool any_derivs = !is_constant_all::value; const auto& log_y = to_ref_if(log(y_val)); const auto& log_sigma = to_ref_if(log(sigma_val)); - const auto& log_pow_n - = to_ref_if(alpha_val * log_y - alpha_val * log_sigma); + const auto& log_y_div_sigma = to_ref_if(log_y - log_sigma); + const auto& log_pow_n = to_ref_if(alpha_val * log_y_div_sigma); const auto& pow_n = to_ref_if(exp(log_pow_n)); const auto& log_cdf_n = to_ref_if(log1m_exp(-pow_n)); @@ -88,7 +88,7 @@ return_type_t weibull_cdf(const T_y& y, } } if (!is_constant_all::value) { - partials<1>(ops_partials) = exp(log_rep_deriv) * (log_y - log_sigma); + partials<1>(ops_partials) = exp(log_rep_deriv) * log_y_div_sigma; } } return ops_partials.build(exp(log_cdf)); diff --git a/stan/math/prim/prob/weibull_lcdf.hpp b/stan/math/prim/prob/weibull_lcdf.hpp index 8152b759c60..cb198894afa 100644 --- a/stan/math/prim/prob/weibull_lcdf.hpp +++ b/stan/math/prim/prob/weibull_lcdf.hpp @@ -68,8 +68,8 @@ return_type_t weibull_lcdf(const T_y& y, constexpr bool any_derivs = !is_constant_all::value; const auto& log_y = to_ref_if(log(y_val)); const auto& log_sigma = to_ref_if(log(sigma_val)); - const auto& log_pow_n - = to_ref_if(alpha_val * log_y - alpha_val * log_sigma); + const auto& log_y_div_sigma = to_ref_if(log_y - log_sigma); + const auto& log_pow_n = to_ref_if(alpha_val * log_y_div_sigma); const auto& pow_n = to_ref_if(exp(log_pow_n)); if (any_derivs) { @@ -85,7 +85,7 @@ return_type_t weibull_lcdf(const T_y& y, } } if (!is_constant_all::value) { - partials<1>(ops_partials) = exp(log_rep_deriv) * (log_y - log_sigma); + partials<1>(ops_partials) = exp(log_rep_deriv) * log_y_div_sigma; } } return ops_partials.build(sum(log1m_exp(-pow_n)));