From 05fe23884359ea1c229160ec1c086580dd028f17 Mon Sep 17 00:00:00 2001 From: gregorp90 Date: Mon, 3 May 2021 17:08:49 +0200 Subject: [PATCH 01/16] Started on loglogistic_cdf. --- stan/math/prim/prob/loglogistic_cdf.hpp | 119 ++++++++++++++++++ .../prob/loglogistic/loglogistic_cdf_test.hpp | 63 ++++++++++ 2 files changed, 182 insertions(+) create mode 100644 stan/math/prim/prob/loglogistic_cdf.hpp create mode 100644 test/prob/loglogistic/loglogistic_cdf_test.hpp diff --git a/stan/math/prim/prob/loglogistic_cdf.hpp b/stan/math/prim/prob/loglogistic_cdf.hpp new file mode 100644 index 00000000000..8d2ab321192 --- /dev/null +++ b/stan/math/prim/prob/loglogistic_cdf.hpp @@ -0,0 +1,119 @@ +#ifndef STAN_MATH_PRIM_PROB_LOGLOGISTIC_CDF_HPP +#define STAN_MATH_PRIM_PROB_LOGLOGISTIC_CDF_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** \ingroup prob_dists + * The log of the loglogistic density for the specified scalar(s) + * given the specified scales(s) and shape(s). y, alpha, or beta + * can each be either a scalar or a vector. Any vector inputs + * must be the same length. + * + *

The result log probability is defined to be the sum of the + * log probabilities for each observation/scale/shape triple. + * + * @tparam T_y type of scalar. + * @tparam T_scale type of scale parameter. + * @tparam T_shape type of shape parameter. + * @param y (Sequence of) scalar(s). + * @param alpha (Sequence of) scale parameter(s) + * for the loglogistic distribution. + * @param beta (Sequence of) shape parameter(s) for the + * loglogistic distribution. + * @return The log of the product of the densities. + * @throw std::domain_error if any of the inputs are not positive and finite. + */ +template * = nullptr> +return_type_t loglogistic_cdf(const T_y& y, + const T_scale& alpha, + const T_shape& beta) { + using T_partials_return = partials_return_t; + using T_y_ref = ref_type_t; + using T_alpha_ref = ref_type_t; + using T_beta_ref = ref_type_t; + using std::pow; + static const char* function = "loglogistic_cdf"; + check_consistent_sizes(function, "Random variable", y, "Scale parameter", + alpha, "Shape parameter", beta); + T_y_ref y_ref = y; + T_alpha_ref alpha_ref = alpha; + T_beta_ref beta_ref = beta; + + 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) beta_val = to_ref(as_value_column_array_or_scalar(beta_ref)); + + + check_positive(function, "Random variable", y_val); + check_positive_finite(function, "Scale parameter", alpha_val); + check_positive_finite(function, "Shape parameter", beta_val); + + if (size_zero(y, alpha, beta)) { + return 1.0; + } + + operands_and_partials ops_partials( + y_ref, alpha_ref, beta_ref); + + T_partials_return cdf(1.0); + cdf *= prod(1 / (1 + pow(alpha_val / y_val, beta_val))); + + if (!is_constant_all::value) { + const auto& y_deriv = (pow(alpha_val, beta_val) * beta_val * + pow(y_val, -beta_val - 1)) / pow(1 + pow(alpha_val / y_val, beta_val), 2); + ops_partials.edge1_.partials_ = y_deriv; + // std::cout << "Partial 1: " << y_deriv << std::endl; + } + if (!is_constant_all::value) { + const auto& alpha_deriv = (-pow(y_val, -beta_val) * beta_val * + pow(alpha_val, beta_val - 1)) / pow(1 + pow(alpha_val / y_val, beta_val), 2); + ops_partials.edge2_.partials_ = alpha_deriv; + // std::cout << "Partial 2: " << alpha_deriv << std::endl; + // std::cout << "Partial 2: " << (alpha_deriv == 0) << std::endl; + } + if (!is_constant_all::value) { + const auto& beta_deriv = (-pow(alpha_val / y_val, beta_val) * + log(alpha_val / y_val)) / + pow(1 + pow(alpha_val / y_val, beta_val), 2); + ops_partials.edge3_.partials_ = beta_deriv; + // Is it the product here, rather than sum, if it is a scalar??? + + + // std::cout << "Partial 3: " << beta_deriv << std::endl; + } + if (!is_constant_all::value) { + for (size_t n = 0; n < stan::math::size(beta); ++n) { + // std::cout << "ii1" << std::endl; + if (ops_partials.edge3_.partials_[n] != ops_partials.edge3_.partials_[n]) { + ops_partials.edge3_.partials_[n] = 0; + // std::cout << "ii2" << std::endl; + } + // ops_partials.edge3_.partials_[n] *= cdf; + } + // std::cout << "Partial 3: " << ops_partials.edge3_.partials_ << std::endl; + } + return ops_partials.build(cdf); +} + +} // namespace math +} // namespace stan +#endif diff --git a/test/prob/loglogistic/loglogistic_cdf_test.hpp b/test/prob/loglogistic/loglogistic_cdf_test.hpp new file mode 100644 index 00000000000..3429ef6e26d --- /dev/null +++ b/test/prob/loglogistic/loglogistic_cdf_test.hpp @@ -0,0 +1,63 @@ +// Arguments: Doubles, Doubles, Doubles +#include + +using stan::math::var; +using std::numeric_limits; +using std::vector; +using std::pow; + +class AgradCdfLogistic : public AgradCdfTest { + public: + void valid_values(vector >& parameters, vector& cdf) { + vector param(3); + + param[0] = 2.0; // y + param[1] = 3.0; // Scale + param[2] = 2.0; // Shape + parameters.push_back(param); + cdf.push_back(0.30769230769230765388); // expected cdf + } + + void invalid_values(vector& index, vector& value) { + // mu + index.push_back(1U); + value.push_back(-numeric_limits::infinity()); + + index.push_back(1U); + value.push_back(numeric_limits::infinity()); + + // sigma + index.push_back(2U); + value.push_back(-1.0); + + index.push_back(2U); + value.push_back(0.0); + + index.push_back(2U); + value.push_back(numeric_limits::infinity()); + } + + bool has_lower_bound() { return false; } + + bool has_upper_bound() { return false; } + + template + stan::return_type_t cdf(const T_y& y, + const T_scale& alpha, + const T_shape& beta, + const T3&, + const T4&, const T5&) { + return stan::math::loglogistic_cdf(y, alpha, beta); + } + + template + stan::return_type_t cdf_function(const T_y& y, + const T_scale& alpha, + const T_shape& beta, + const T3&, const T4&, + const T5&) { + return 1.0 / (1.0 + pow(alpha / y, beta)); + } +}; From 5c283860a3f524ef5fa0f9320400cda4fc22c6f6 Mon Sep 17 00:00:00 2001 From: gregorp90 Date: Tue, 4 May 2021 16:11:59 +0200 Subject: [PATCH 02/16] Added lower bound to test. --- .../prob/loglogistic/loglogistic_cdf_test.hpp | 28 ++++++++++++++++--- 1 file changed, 24 insertions(+), 4 deletions(-) diff --git a/test/prob/loglogistic/loglogistic_cdf_test.hpp b/test/prob/loglogistic/loglogistic_cdf_test.hpp index 3429ef6e26d..7fb2a2d9074 100644 --- a/test/prob/loglogistic/loglogistic_cdf_test.hpp +++ b/test/prob/loglogistic/loglogistic_cdf_test.hpp @@ -16,17 +16,35 @@ class AgradCdfLogistic : public AgradCdfTest { param[2] = 2.0; // Shape parameters.push_back(param); cdf.push_back(0.30769230769230765388); // expected cdf + + + param[0] = 20.0; // y + param[1] = 11.5; // Scale + param[2] = 15.0; // Shape + parameters.push_back(param); + cdf.push_back(0.99975173823523311167); // expected cdf } void invalid_values(vector& index, vector& value) { - // mu - index.push_back(1U); + // y + index.push_back(0U); + value.push_back(-1.0); + + index.push_back(0U); value.push_back(-numeric_limits::infinity()); + + // alpha + index.push_back(1U); + value.push_back(-1); + + index.push_back(1U); + value.push_back(0); + index.push_back(1U); value.push_back(numeric_limits::infinity()); - // sigma + // beta index.push_back(2U); value.push_back(-1.0); @@ -37,7 +55,9 @@ class AgradCdfLogistic : public AgradCdfTest { value.push_back(numeric_limits::infinity()); } - bool has_lower_bound() { return false; } + bool has_lower_bound() { return true; } + + double lower_bound() { return 0.0; } bool has_upper_bound() { return false; } From 04ebbc5efab096a83679f6d644cd378b5b24b3ca Mon Sep 17 00:00:00 2001 From: gregorp90 Date: Tue, 4 May 2021 16:12:24 +0200 Subject: [PATCH 03/16] Fixed derivatives. --- stan/math/prim/prob.hpp | 1 + stan/math/prim/prob/loglogistic_cdf.hpp | 13 ++++++++----- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/stan/math/prim/prob.hpp b/stan/math/prim/prob.hpp index 04416a966ec..e7455a88ad7 100644 --- a/stan/math/prim/prob.hpp +++ b/stan/math/prim/prob.hpp @@ -175,6 +175,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/prob/loglogistic_cdf.hpp b/stan/math/prim/prob/loglogistic_cdf.hpp index 8d2ab321192..2c1a5ea3105 100644 --- a/stan/math/prim/prob/loglogistic_cdf.hpp +++ b/stan/math/prim/prob/loglogistic_cdf.hpp @@ -63,7 +63,7 @@ return_type_t loglogistic_cdf(const T_y& y, decltype(auto) beta_val = to_ref(as_value_column_array_or_scalar(beta_ref)); - check_positive(function, "Random variable", y_val); + check_nonnegative(function, "Random variable", y_val); check_positive_finite(function, "Scale parameter", alpha_val); check_positive_finite(function, "Shape parameter", beta_val); @@ -74,19 +74,22 @@ return_type_t loglogistic_cdf(const T_y& y, operands_and_partials ops_partials( y_ref, alpha_ref, beta_ref); + const auto& cdf_elementwise = to_ref_if::value>(1 / (1 + pow(alpha_val / y_val, beta_val))); + T_partials_return cdf(1.0); - cdf *= prod(1 / (1 + pow(alpha_val / y_val, beta_val))); + cdf *= prod(cdf_elementwise); if (!is_constant_all::value) { const auto& y_deriv = (pow(alpha_val, beta_val) * beta_val * pow(y_val, -beta_val - 1)) / pow(1 + pow(alpha_val / y_val, beta_val), 2); - ops_partials.edge1_.partials_ = y_deriv; + ops_partials.edge1_.partials_ = (y_deriv / cdf_elementwise) * cdf; // std::cout << "Partial 1: " << y_deriv << std::endl; } if (!is_constant_all::value) { const auto& alpha_deriv = (-pow(y_val, -beta_val) * beta_val * pow(alpha_val, beta_val - 1)) / pow(1 + pow(alpha_val / y_val, beta_val), 2); - ops_partials.edge2_.partials_ = alpha_deriv; + ops_partials.edge2_.partials_ = (alpha_deriv / cdf_elementwise) * cdf; // std::cout << "Partial 2: " << alpha_deriv << std::endl; // std::cout << "Partial 2: " << (alpha_deriv == 0) << std::endl; } @@ -94,7 +97,7 @@ return_type_t loglogistic_cdf(const T_y& y, const auto& beta_deriv = (-pow(alpha_val / y_val, beta_val) * log(alpha_val / y_val)) / pow(1 + pow(alpha_val / y_val, beta_val), 2); - ops_partials.edge3_.partials_ = beta_deriv; + ops_partials.edge3_.partials_ = (beta_deriv / cdf_elementwise) * cdf; // Is it the product here, rather than sum, if it is a scalar??? From 1af99c1f7c69846b77fc3664d61628d1a059c62b Mon Sep 17 00:00:00 2001 From: gregorp90 Date: Wed, 5 May 2021 15:09:17 +0200 Subject: [PATCH 04/16] Temporary debugging. --- stan/math/prim/prob/loglogistic_cdf.hpp | 50 ++++++++++++++++++------- 1 file changed, 36 insertions(+), 14 deletions(-) diff --git a/stan/math/prim/prob/loglogistic_cdf.hpp b/stan/math/prim/prob/loglogistic_cdf.hpp index 2c1a5ea3105..d2a60984a3d 100644 --- a/stan/math/prim/prob/loglogistic_cdf.hpp +++ b/stan/math/prim/prob/loglogistic_cdf.hpp @@ -81,39 +81,61 @@ return_type_t loglogistic_cdf(const T_y& y, cdf *= prod(cdf_elementwise); if (!is_constant_all::value) { + std::cout << "Partial 1: " << ops_partials.edge1_.partials_[0] << std::endl; const auto& y_deriv = (pow(alpha_val, beta_val) * beta_val * pow(y_val, -beta_val - 1)) / pow(1 + pow(alpha_val / y_val, beta_val), 2); ops_partials.edge1_.partials_ = (y_deriv / cdf_elementwise) * cdf; - // std::cout << "Partial 1: " << y_deriv << std::endl; + std::cout << "Partial 1: " << (y_deriv / cdf_elementwise) * cdf << std::endl; } if (!is_constant_all::value) { const auto& alpha_deriv = (-pow(y_val, -beta_val) * beta_val * pow(alpha_val, beta_val - 1)) / pow(1 + pow(alpha_val / y_val, beta_val), 2); ops_partials.edge2_.partials_ = (alpha_deriv / cdf_elementwise) * cdf; - // std::cout << "Partial 2: " << alpha_deriv << std::endl; + std::cout << "Partial 2: " << (alpha_deriv / cdf_elementwise) * cdf << std::endl; // std::cout << "Partial 2: " << (alpha_deriv == 0) << std::endl; } if (!is_constant_all::value) { - const auto& beta_deriv = (-pow(alpha_val / y_val, beta_val) * - log(alpha_val / y_val)) / + const auto& beta_deriv = multiply_log(-pow(alpha_val / y_val, beta_val), + (alpha_val / y_val)) / pow(1 + pow(alpha_val / y_val, beta_val), 2); + // const auto& beta_deriv = (-pow(alpha_val / y_val, beta_val) * + // log(alpha_val / y_val)) / + // pow(1 + pow(alpha_val / y_val, beta_val), 2); ops_partials.edge3_.partials_ = (beta_deriv / cdf_elementwise) * cdf; // Is it the product here, rather than sum, if it is a scalar??? - // std::cout << "Partial 3: " << beta_deriv << std::endl; + std::cout << "Partial 3: " << (beta_deriv / cdf_elementwise) * cdf << std::endl; } + std::cout << "Partial edge y 1: " << ops_partials.edge1_.partials_[0] << std::endl; + std::cout << "Partial edge y 2: " << ops_partials.edge1_.partials_[1] << std::endl; + if (!is_constant_all::value) { - for (size_t n = 0; n < stan::math::size(beta); ++n) { - // std::cout << "ii1" << std::endl; - if (ops_partials.edge3_.partials_[n] != ops_partials.edge3_.partials_[n]) { - ops_partials.edge3_.partials_[n] = 0; - // std::cout << "ii2" << std::endl; - } - // ops_partials.edge3_.partials_[n] *= cdf; - } - // std::cout << "Partial 3: " << ops_partials.edge3_.partials_ << std::endl; + // std::cout << "Partial edge y 1: " << ops_partials.edge1_.partials_[0] << std::endl; + // std::cout << "Partial edge y 2: " << ops_partials.edge1_.partials_[1] << std::endl; + std::cout << "Partial edge alpha 1 " << ops_partials.edge2_.partials_[0] << std::endl; + std::cout << "Partial edge beta 1 " << ops_partials.edge3_.partials_[0] << std::endl; + std::cout << "Partial edge beta 2 " << ops_partials.edge3_.partials_[1] << std::endl; } + + // We could move this to above if? Probably. Above ops_partials, so that + // we don't have to access ops_partials. + // if (!is_constant_all::value) { + // for (size_t n = 0; n < stan::math::size(beta); ++n) { // Beta is not right here? + // // We should have something of larger size? + // // std::cout << "ii1" << std::endl; + // std::cout << n << std::endl; + // std::cout << "Partial edge: " << ops_partials.edge1_.partials_[n] << std::endl; + // std::cout << "Partial edge: " << ops_partials.edge2_.partials_[n] << std::endl; + // std::cout << "Partial edge: " << ops_partials.edge3_.partials_[n] << std::endl; + // // if (ops_partials.edge3_.partials_[n] != ops_partials.edge3_.partials_[n]) { + // // ops_partials.edge3_.partials_[n] = 0; + // // // std::cout << "ii2" << std::endl; + // // } + // // ops_partials.edge3_.partials_[n] *= cdf; + // } + // std::cout << "Partial 3: " << ops_partials.edge3_.partials_ << std::endl; + // } return ops_partials.build(cdf); } From d9694d43c9c334da3caebff59726390a02f6a878 Mon Sep 17 00:00:00 2001 From: gregorp90 Date: Wed, 5 May 2021 16:05:45 +0200 Subject: [PATCH 05/16] Optimized calculations. --- stan/math/prim/prob/loglogistic_cdf.hpp | 97 ++++++++++--------------- 1 file changed, 40 insertions(+), 57 deletions(-) diff --git a/stan/math/prim/prob/loglogistic_cdf.hpp b/stan/math/prim/prob/loglogistic_cdf.hpp index d2a60984a3d..393625f0b78 100644 --- a/stan/math/prim/prob/loglogistic_cdf.hpp +++ b/stan/math/prim/prob/loglogistic_cdf.hpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include @@ -74,68 +75,50 @@ return_type_t loglogistic_cdf(const T_y& y, operands_and_partials ops_partials( y_ref, alpha_ref, beta_ref); - const auto& cdf_elementwise = to_ref_if::value>(1 / (1 + pow(alpha_val / y_val, beta_val))); - - T_partials_return cdf(1.0); - cdf *= prod(cdf_elementwise); - - if (!is_constant_all::value) { - std::cout << "Partial 1: " << ops_partials.edge1_.partials_[0] << std::endl; - const auto& y_deriv = (pow(alpha_val, beta_val) * beta_val * - pow(y_val, -beta_val - 1)) / pow(1 + pow(alpha_val / y_val, beta_val), 2); - ops_partials.edge1_.partials_ = (y_deriv / cdf_elementwise) * cdf; - std::cout << "Partial 1: " << (y_deriv / cdf_elementwise) * cdf << std::endl; - } - if (!is_constant_all::value) { - const auto& alpha_deriv = (-pow(y_val, -beta_val) * beta_val * - pow(alpha_val, beta_val - 1)) / pow(1 + pow(alpha_val / y_val, beta_val), 2); - ops_partials.edge2_.partials_ = (alpha_deriv / cdf_elementwise) * cdf; - std::cout << "Partial 2: " << (alpha_deriv / cdf_elementwise) * cdf << std::endl; - // std::cout << "Partial 2: " << (alpha_deriv == 0) << std::endl; + if (sum(promote_scalar(y_val == 0))) { + return ops_partials.build(0.0); } - if (!is_constant_all::value) { - const auto& beta_deriv = multiply_log(-pow(alpha_val / y_val, beta_val), - (alpha_val / y_val)) / - pow(1 + pow(alpha_val / y_val, beta_val), 2); - // const auto& beta_deriv = (-pow(alpha_val / y_val, beta_val) * - // log(alpha_val / y_val)) / - // pow(1 + pow(alpha_val / y_val, beta_val), 2); - ops_partials.edge3_.partials_ = (beta_deriv / cdf_elementwise) * cdf; - // Is it the product here, rather than sum, if it is a scalar??? + const auto& alpha_div_y = to_ref_if::value>( + alpha_val / y_val); + const auto& alpha_div_y_pow_beta = to_ref_if::value>(pow(alpha_div_y, beta_val)); + const auto& prod_all = to_ref_if::value>(1 / (1 + alpha_div_y_pow_beta)); - std::cout << "Partial 3: " << (beta_deriv / cdf_elementwise) * cdf << std::endl; - } - std::cout << "Partial edge y 1: " << ops_partials.edge1_.partials_[0] << std::endl; - std::cout << "Partial edge y 2: " << ops_partials.edge1_.partials_[1] << std::endl; - - if (!is_constant_all::value) { - // std::cout << "Partial edge y 1: " << ops_partials.edge1_.partials_[0] << std::endl; - // std::cout << "Partial edge y 2: " << ops_partials.edge1_.partials_[1] << std::endl; - std::cout << "Partial edge alpha 1 " << ops_partials.edge2_.partials_[0] << std::endl; - std::cout << "Partial edge beta 1 " << ops_partials.edge3_.partials_[0] << std::endl; - std::cout << "Partial edge beta 2 " << ops_partials.edge3_.partials_[1] << std::endl; + T_partials_return cdf(1.0); + cdf *= prod(prod_all); + + if (!is_constant_all::value) { + const auto& prod_all_sq = to_ref_if::value + + !is_constant_all::value + + !is_constant_all::value + >= 2>(prod_all * prod_all); + const auto& cdf_div_elt = to_ref_if::value + + !is_constant_all::value + + !is_constant_all::value + >= 2>(cdf / prod_all); + if (!is_constant_all::value) { + const auto& alpha_div_times_beta = + to_ref_if::value + + !is_constant_all::value == 2>(alpha_div_y_pow_beta * beta_val); + if (!is_constant_all::value) { + const auto& y_deriv = alpha_div_times_beta * inv(y_val) * prod_all_sq; + ops_partials.edge1_.partials_ = y_deriv * cdf_div_elt; + } + if (!is_constant_all::value) { + const auto& alpha_deriv = - alpha_div_times_beta * inv(alpha_val) * + prod_all_sq; + ops_partials.edge2_.partials_ = alpha_deriv * cdf_div_elt; + } + } + if (!is_constant_all::value) { + const auto& beta_deriv = -multiply_log(alpha_div_y_pow_beta, + alpha_div_y) * prod_all_sq; + ops_partials.edge3_.partials_ = beta_deriv * cdf_div_elt; + } } - // We could move this to above if? Probably. Above ops_partials, so that - // we don't have to access ops_partials. - // if (!is_constant_all::value) { - // for (size_t n = 0; n < stan::math::size(beta); ++n) { // Beta is not right here? - // // We should have something of larger size? - // // std::cout << "ii1" << std::endl; - // std::cout << n << std::endl; - // std::cout << "Partial edge: " << ops_partials.edge1_.partials_[n] << std::endl; - // std::cout << "Partial edge: " << ops_partials.edge2_.partials_[n] << std::endl; - // std::cout << "Partial edge: " << ops_partials.edge3_.partials_[n] << std::endl; - // // if (ops_partials.edge3_.partials_[n] != ops_partials.edge3_.partials_[n]) { - // // ops_partials.edge3_.partials_[n] = 0; - // // // std::cout << "ii2" << std::endl; - // // } - // // ops_partials.edge3_.partials_[n] *= cdf; - // } - // std::cout << "Partial 3: " << ops_partials.edge3_.partials_ << std::endl; - // } return ops_partials.build(cdf); } From 7cd9f2cb5675d7c30c0562f256c1b14c6d3adeff Mon Sep 17 00:00:00 2001 From: gregorp90 Date: Wed, 5 May 2021 16:11:15 +0200 Subject: [PATCH 06/16] Added test case. --- test/prob/loglogistic/loglogistic_cdf_test.hpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/test/prob/loglogistic/loglogistic_cdf_test.hpp b/test/prob/loglogistic/loglogistic_cdf_test.hpp index 7fb2a2d9074..d6aa7434dd4 100644 --- a/test/prob/loglogistic/loglogistic_cdf_test.hpp +++ b/test/prob/loglogistic/loglogistic_cdf_test.hpp @@ -23,6 +23,12 @@ class AgradCdfLogistic : public AgradCdfTest { param[2] = 15.0; // Shape parameters.push_back(param); cdf.push_back(0.99975173823523311167); // expected cdf + + param[0] = 0.0001; // y + param[1] = 1.2; // Scale + param[2] = 1.0; // Shape + parameters.push_back(param); + cdf.push_back(0.00008332638946754438); // expected cdf } void invalid_values(vector& index, vector& value) { From 305e34332edd51a546b47a036277f60ab9a93b3e Mon Sep 17 00:00:00 2001 From: gregorp90 Date: Thu, 6 May 2021 17:00:22 +0200 Subject: [PATCH 07/16] Added loglogistic.rng. --- stan/math/prim/prob.hpp | 1 + stan/math/prim/prob/loglogistic_rng.hpp | 67 +++++++++++++++++++++++++ 2 files changed, 68 insertions(+) create mode 100644 stan/math/prim/prob/loglogistic_rng.hpp diff --git a/stan/math/prim/prob.hpp b/stan/math/prim/prob.hpp index e7455a88ad7..2fce0b34e1d 100644 --- a/stan/math/prim/prob.hpp +++ b/stan/math/prim/prob.hpp @@ -177,6 +177,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/prob/loglogistic_rng.hpp b/stan/math/prim/prob/loglogistic_rng.hpp new file mode 100644 index 00000000000..62c212b1a69 --- /dev/null +++ b/stan/math/prim/prob/loglogistic_rng.hpp @@ -0,0 +1,67 @@ +#ifndef STAN_MATH_PRIM_PROB_LOGLOGISTIC_RNG_HPP +#define STAN_MATH_PRIM_PROB_LOGLOGISTIC_RNG_HPP + +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** \ingroup prob_dists + * Return a gamma random variate for the given shape and inverse + * scale parameters using the specified random number generator. + * + * alpha and beta can each be a scalar or a one-dimensional container. Any + * non-scalar inputs must be the same size. + * + * @tparam T_shape type of shape parameter + * @tparam T_inv type of inverse scale parameter + * @tparam RNG type of random number generator + * @param alpha (Sequence of) positive shape parameter(s) + * @param beta (Sequence of) positive inverse scale parameter(s) + * @param rng random number generator + * @return (Sequence of) gamma random variate(s) + * @throw std::domain_error if alpha or beta are nonpositive + * @throw std::invalid_argument if non-scalar arguments are of different + * sizes + */ +template +inline typename VectorBuilder::type loglogistic_rng( + const T_scale& alpha, const T_shape& beta, RNG& rng) { + using boost::uniform_01; + using boost::variate_generator; + using std::pow; + using T_alpha_ref = ref_type_t; + using T_beta_ref = ref_type_t; + static const char* function = "loglogistic_rng"; + check_consistent_sizes(function, "Scale parameter", alpha, + "Shape Parameter", beta); + T_alpha_ref alpha_ref = alpha; + T_beta_ref beta_ref = beta; + check_positive_finite(function, "Scale parameter", alpha_ref); + check_positive_finite(function, "Shape parameter", beta_ref); + + scalar_seq_view alpha_vec(alpha_ref); + scalar_seq_view beta_vec(beta_ref); + size_t N = max_size(alpha, beta); + VectorBuilder output(N); + + float tmp; + for (size_t n = 0; n < N; ++n) { + variate_generator > uniform01_rng( + rng, uniform_01<>()); + tmp = uniform01_rng(); + output[n] = alpha_vec[n] * pow(tmp / (1 - tmp), + 1 / beta_vec[n]); + } + return output.data(); +} + +} // namespace math +} // namespace stan +#endif From d207e3306e96030ed9fec9ddf17b86ef44ca7111 Mon Sep 17 00:00:00 2001 From: gregorp90 Date: Thu, 6 May 2021 17:16:00 +0200 Subject: [PATCH 08/16] Started on rng test for loglogistic. --- test/unit/math/prim/prob/loglogistic_test.cpp | 83 +++++++++++++++++++ 1 file changed, 83 insertions(+) create mode 100644 test/unit/math/prim/prob/loglogistic_test.cpp diff --git a/test/unit/math/prim/prob/loglogistic_test.cpp b/test/unit/math/prim/prob/loglogistic_test.cpp new file mode 100644 index 00000000000..f470822b211 --- /dev/null +++ b/test/unit/math/prim/prob/loglogistic_test.cpp @@ -0,0 +1,83 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +class LoglogisticTestRig : public VectorRNGTestRig { + public: + LoglogisticTestRig() + : VectorRNGTestRig(10000, 10, + {2.5, 1.7, 0.2, 0.1, 2.0}, + {3, 2, 1, 5, 10, 6}, + {-2.5, -1.7, -0.2, -0.1, 0.0}, + {-3, -2, -1, -4, -10, 0}, + {0.1, 1.0, 2.5, 4.0}, + {1, 2, 3, 4}, + {-2.7, -1.5, -0.5, 0.0}, + {-3, -2, -1, 0}) {} + + template + auto generate_samples(const T1& alpha, const T2& beta) const { + return stan::math::loglogistic_rng(alpha, beta, rng); + } +}; + +double icdf(double x, double alpha, double beta) { + return alpha * pow(x / (1 - x), 1 / beta); +} + +TEST(ProbDistributionsLoglogistic, errorCheck) { + check_dist_throws_all_types(LoglogisticTestRig()); +} + +TEST(ProbDistributionsLoglogistic, error_check) { + boost::random::mt19937 rng; + EXPECT_NO_THROW(stan::math::loglogistic_rng(10.0, 2.0, rng)); + + EXPECT_THROW(stan::math::loglogistic_rng(2.0, -1.0, rng), + std::domain_error); + EXPECT_THROW(stan::math::loglogistic_rng(-2.0, 1.0, rng), + std::domain_error); + EXPECT_THROW(stan::math::loglogistic_rng( + 10, stan::math::positive_infinity(), rng), + std::domain_error); + EXPECT_THROW(stan::math::loglogistic_rng( + stan::math::positive_infinity(), 2, rng), + std::domain_error); +} + +TEST(ProbDistributionsLoglogistic, test_sampling_icdf) { + for (double p : {0.0, 0.1, 0.2, 0.5, 0.7, 0.9, 0.99}) { + for (double alpha : {1.11, 0.13, 1.2, 4.67}) { + for (double beta : {0.11, 1.33, 2.0, 3.2}) { + double x = icdf(p, alpha, beta); + EXPECT_FLOAT_EQ(stan::math::loglogistic_cdf(x, alpha, beta), p); + } + } + } +} + +TEST(ProbDistributionsLoglogistic, chiSquareGoodnessFitTest) { + boost::random::mt19937 rng; + int N = 10000; + int K = stan::math::round(2 * std::pow(N, 0.4)); + + std::vector samples; + for (int i = 0; i < N; ++i) { + samples.push_back( + stan::math::loglogistic_rng(1.2, 2.0, rng)); + } + std::vector quantiles; + for (int i = 1; i < K; ++i) { + double frac = static_cast(i) / K; + quantiles.push_back(icdf(frac, 1.2, 2.0)); + } + quantiles.push_back(std::numeric_limits::max()); + + // Assert that they match + assert_matches_quantiles(samples, quantiles, 1e-6); +} From 8cd982f5b911c64776bdf7a6183562407625381c Mon Sep 17 00:00:00 2001 From: gregorp90 Date: Fri, 7 May 2021 10:34:34 +0200 Subject: [PATCH 09/16] Fixed documentation for loglogistic_cdf. --- stan/math/prim/prob/loglogistic_cdf.hpp | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/stan/math/prim/prob/loglogistic_cdf.hpp b/stan/math/prim/prob/loglogistic_cdf.hpp index 393625f0b78..00b58aae4a5 100644 --- a/stan/math/prim/prob/loglogistic_cdf.hpp +++ b/stan/math/prim/prob/loglogistic_cdf.hpp @@ -22,13 +22,11 @@ namespace stan { namespace math { /** \ingroup prob_dists - * The log of the loglogistic density for the specified scalar(s) - * given the specified scales(s) and shape(s). y, alpha, or beta - * can each be either a scalar or a vector. Any vector inputs + * The loglogistic cumulative distribution function for the specified + * scalar(s) given the specified scales(s) and shape(s). y, alpha, or + * beta can each be either a scalar or a vector. Any vector inputs * must be the same length. * - *

The result log probability is defined to be the sum of the - * log probabilities for each observation/scale/shape triple. * * @tparam T_y type of scalar. * @tparam T_scale type of scale parameter. @@ -38,8 +36,9 @@ namespace math { * for the loglogistic distribution. * @param beta (Sequence of) shape parameter(s) for the * loglogistic distribution. - * @return The log of the product of the densities. - * @throw std::domain_error if any of the inputs are not positive and finite. + * @return The loglogistic cdf evaluated at the specified arguments. + * @throw std::domain_error if any of the inputs are not positive or + * if and of the parameters are not finite. */ template Date: Fri, 7 May 2021 10:35:25 +0200 Subject: [PATCH 10/16] loglogistic_rng tests working. --- test/unit/math/prim/prob/loglogistic_test.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/unit/math/prim/prob/loglogistic_test.cpp b/test/unit/math/prim/prob/loglogistic_test.cpp index f470822b211..d6dbe207244 100644 --- a/test/unit/math/prim/prob/loglogistic_test.cpp +++ b/test/unit/math/prim/prob/loglogistic_test.cpp @@ -21,7 +21,8 @@ class LoglogisticTestRig : public VectorRNGTestRig { {-3, -2, -1, 0}) {} template - auto generate_samples(const T1& alpha, const T2& beta) const { + auto generate_samples(const T1& alpha, const T2& beta, + const T3& unused, T_rng& rng) const { return stan::math::loglogistic_rng(alpha, beta, rng); } }; From 04012c6e03249fc4ceb52f9350b015acaa97f761 Mon Sep 17 00:00:00 2001 From: gregorp90 Date: Fri, 7 May 2021 12:07:07 +0200 Subject: [PATCH 11/16] Minor changes and documentation for loglogistic_rng. --- stan/math/prim/prob/loglogistic_rng.hpp | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/stan/math/prim/prob/loglogistic_rng.hpp b/stan/math/prim/prob/loglogistic_rng.hpp index 62c212b1a69..b3dd0c2b604 100644 --- a/stan/math/prim/prob/loglogistic_rng.hpp +++ b/stan/math/prim/prob/loglogistic_rng.hpp @@ -12,19 +12,19 @@ namespace stan { namespace math { /** \ingroup prob_dists - * Return a gamma random variate for the given shape and inverse - * scale parameters using the specified random number generator. + * Return a loglogistic random variate for the given scale and + * shape parameters using the specified random number generator. * * alpha and beta can each be a scalar or a one-dimensional container. Any * non-scalar inputs must be the same size. * + * @tparam T_scale type of scale parameter * @tparam T_shape type of shape parameter - * @tparam T_inv type of inverse scale parameter * @tparam RNG type of random number generator - * @param alpha (Sequence of) positive shape parameter(s) - * @param beta (Sequence of) positive inverse scale parameter(s) + * @param alpha (Sequence of) positive scale parameter(s) + * @param beta (Sequence of) positive shape parameter(s) * @param rng random number generator - * @return (Sequence of) gamma random variate(s) + * @return (Sequence of) loglogistic random variate(s) * @throw std::domain_error if alpha or beta are nonpositive * @throw std::invalid_argument if non-scalar arguments are of different * sizes @@ -51,11 +51,10 @@ inline typename VectorBuilder output(N); - float tmp; for (size_t n = 0; n < N; ++n) { variate_generator > uniform01_rng( rng, uniform_01<>()); - tmp = uniform01_rng(); + const double tmp = uniform01_rng(); output[n] = alpha_vec[n] * pow(tmp / (1 - tmp), 1 / beta_vec[n]); } From 70782e6c1653027665970ab9ae153eee26bc8dc9 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Fri, 7 May 2021 12:03:36 +0000 Subject: [PATCH 12/16] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/prob/loglogistic_cdf.hpp | 45 ++++++++++--------- stan/math/prim/prob/loglogistic_rng.hpp | 15 +++---- .../prob/loglogistic/loglogistic_cdf_test.hpp | 11 ++--- test/unit/math/prim/prob/loglogistic_test.cpp | 37 +++++++-------- 4 files changed, 48 insertions(+), 60 deletions(-) diff --git a/stan/math/prim/prob/loglogistic_cdf.hpp b/stan/math/prim/prob/loglogistic_cdf.hpp index 00b58aae4a5..9991c932c28 100644 --- a/stan/math/prim/prob/loglogistic_cdf.hpp +++ b/stan/math/prim/prob/loglogistic_cdf.hpp @@ -41,8 +41,8 @@ namespace math { * if and of the parameters are not finite. */ template * = nullptr> + require_all_not_nonscalar_prim_or_rev_kernel_expression_t< + T_y, T_scale, T_shape>* = nullptr> return_type_t loglogistic_cdf(const T_y& y, const T_scale& alpha, const T_shape& beta) { @@ -62,7 +62,6 @@ return_type_t loglogistic_cdf(const T_y& y, decltype(auto) alpha_val = to_ref(as_value_column_array_or_scalar(alpha_ref)); decltype(auto) beta_val = to_ref(as_value_column_array_or_scalar(beta_ref)); - check_nonnegative(function, "Random variable", y_val); check_positive_finite(function, "Scale parameter", alpha_val); check_positive_finite(function, "Shape parameter", beta_val); @@ -78,42 +77,44 @@ return_type_t loglogistic_cdf(const T_y& y, return ops_partials.build(0.0); } - const auto& alpha_div_y = to_ref_if::value>( - alpha_val / y_val); - const auto& alpha_div_y_pow_beta = to_ref_if::value>(pow(alpha_div_y, beta_val)); - const auto& prod_all = to_ref_if::value>(1 / (1 + alpha_div_y_pow_beta)); + const auto& alpha_div_y + = to_ref_if::value>(alpha_val / y_val); + const auto& alpha_div_y_pow_beta + = to_ref_if::value>( + pow(alpha_div_y, beta_val)); + const auto& prod_all + = to_ref_if::value>( + 1 / (1 + alpha_div_y_pow_beta)); T_partials_return cdf(1.0); cdf *= prod(prod_all); if (!is_constant_all::value) { const auto& prod_all_sq = to_ref_if::value - + !is_constant_all::value - + !is_constant_all::value - >= 2>(prod_all * prod_all); + + !is_constant_all::value + + !is_constant_all::value + >= 2>(prod_all * prod_all); const auto& cdf_div_elt = to_ref_if::value - + !is_constant_all::value - + !is_constant_all::value - >= 2>(cdf / prod_all); + + !is_constant_all::value + + !is_constant_all::value + >= 2>(cdf / prod_all); if (!is_constant_all::value) { - const auto& alpha_div_times_beta = - to_ref_if::value + - !is_constant_all::value == 2>(alpha_div_y_pow_beta * beta_val); + const auto& alpha_div_times_beta = to_ref_if< + !is_constant_all::value + !is_constant_all::value == 2>( + alpha_div_y_pow_beta * beta_val); if (!is_constant_all::value) { const auto& y_deriv = alpha_div_times_beta * inv(y_val) * prod_all_sq; ops_partials.edge1_.partials_ = y_deriv * cdf_div_elt; } if (!is_constant_all::value) { - const auto& alpha_deriv = - alpha_div_times_beta * inv(alpha_val) * - prod_all_sq; + const auto& alpha_deriv + = -alpha_div_times_beta * inv(alpha_val) * prod_all_sq; ops_partials.edge2_.partials_ = alpha_deriv * cdf_div_elt; } } if (!is_constant_all::value) { - const auto& beta_deriv = -multiply_log(alpha_div_y_pow_beta, - alpha_div_y) * prod_all_sq; + const auto& beta_deriv + = -multiply_log(alpha_div_y_pow_beta, alpha_div_y) * prod_all_sq; ops_partials.edge3_.partials_ = beta_deriv * cdf_div_elt; } } diff --git a/stan/math/prim/prob/loglogistic_rng.hpp b/stan/math/prim/prob/loglogistic_rng.hpp index b3dd0c2b604..06b43e7e3a1 100644 --- a/stan/math/prim/prob/loglogistic_rng.hpp +++ b/stan/math/prim/prob/loglogistic_rng.hpp @@ -30,17 +30,16 @@ namespace math { * sizes */ template -inline typename VectorBuilder::type loglogistic_rng( - const T_scale& alpha, const T_shape& beta, RNG& rng) { +inline typename VectorBuilder::type +loglogistic_rng(const T_scale& alpha, const T_shape& beta, RNG& rng) { using boost::uniform_01; using boost::variate_generator; using std::pow; using T_alpha_ref = ref_type_t; using T_beta_ref = ref_type_t; static const char* function = "loglogistic_rng"; - check_consistent_sizes(function, "Scale parameter", alpha, - "Shape Parameter", beta); + check_consistent_sizes(function, "Scale parameter", alpha, "Shape Parameter", + beta); T_alpha_ref alpha_ref = alpha; T_beta_ref beta_ref = beta; check_positive_finite(function, "Scale parameter", alpha_ref); @@ -52,11 +51,9 @@ inline typename VectorBuilder output(N); for (size_t n = 0; n < N; ++n) { - variate_generator > uniform01_rng( - rng, uniform_01<>()); + variate_generator > uniform01_rng(rng, uniform_01<>()); const double tmp = uniform01_rng(); - output[n] = alpha_vec[n] * pow(tmp / (1 - tmp), - 1 / beta_vec[n]); + output[n] = alpha_vec[n] * pow(tmp / (1 - tmp), 1 / beta_vec[n]); } return output.data(); } diff --git a/test/prob/loglogistic/loglogistic_cdf_test.hpp b/test/prob/loglogistic/loglogistic_cdf_test.hpp index d6aa7434dd4..4f36f6316eb 100644 --- a/test/prob/loglogistic/loglogistic_cdf_test.hpp +++ b/test/prob/loglogistic/loglogistic_cdf_test.hpp @@ -3,8 +3,8 @@ using stan::math::var; using std::numeric_limits; -using std::vector; using std::pow; +using std::vector; class AgradCdfLogistic : public AgradCdfTest { public: @@ -17,7 +17,6 @@ class AgradCdfLogistic : public AgradCdfTest { parameters.push_back(param); cdf.push_back(0.30769230769230765388); // expected cdf - param[0] = 20.0; // y param[1] = 11.5; // Scale param[2] = 15.0; // Shape @@ -25,8 +24,8 @@ class AgradCdfLogistic : public AgradCdfTest { cdf.push_back(0.99975173823523311167); // expected cdf param[0] = 0.0001; // y - param[1] = 1.2; // Scale - param[2] = 1.0; // Shape + param[1] = 1.2; // Scale + param[2] = 1.0; // Shape parameters.push_back(param); cdf.push_back(0.00008332638946754438); // expected cdf } @@ -39,7 +38,6 @@ class AgradCdfLogistic : public AgradCdfTest { index.push_back(0U); value.push_back(-numeric_limits::infinity()); - // alpha index.push_back(1U); value.push_back(-1); @@ -71,8 +69,7 @@ class AgradCdfLogistic : public AgradCdfTest { typename T4, typename T5> stan::return_type_t cdf(const T_y& y, const T_scale& alpha, - const T_shape& beta, - const T3&, + const T_shape& beta, const T3&, const T4&, const T5&) { return stan::math::loglogistic_cdf(y, alpha, beta); } diff --git a/test/unit/math/prim/prob/loglogistic_test.cpp b/test/unit/math/prim/prob/loglogistic_test.cpp index d6dbe207244..94a38ea4cc7 100644 --- a/test/unit/math/prim/prob/loglogistic_test.cpp +++ b/test/unit/math/prim/prob/loglogistic_test.cpp @@ -10,19 +10,15 @@ class LoglogisticTestRig : public VectorRNGTestRig { public: LoglogisticTestRig() - : VectorRNGTestRig(10000, 10, - {2.5, 1.7, 0.2, 0.1, 2.0}, - {3, 2, 1, 5, 10, 6}, - {-2.5, -1.7, -0.2, -0.1, 0.0}, - {-3, -2, -1, -4, -10, 0}, - {0.1, 1.0, 2.5, 4.0}, - {1, 2, 3, 4}, - {-2.7, -1.5, -0.5, 0.0}, + : VectorRNGTestRig(10000, 10, {2.5, 1.7, 0.2, 0.1, 2.0}, + {3, 2, 1, 5, 10, 6}, {-2.5, -1.7, -0.2, -0.1, 0.0}, + {-3, -2, -1, -4, -10, 0}, {0.1, 1.0, 2.5, 4.0}, + {1, 2, 3, 4}, {-2.7, -1.5, -0.5, 0.0}, {-3, -2, -1, 0}) {} template - auto generate_samples(const T1& alpha, const T2& beta, - const T3& unused, T_rng& rng) const { + auto generate_samples(const T1& alpha, const T2& beta, const T3& unused, + T_rng& rng) const { return stan::math::loglogistic_rng(alpha, beta, rng); } }; @@ -39,16 +35,14 @@ TEST(ProbDistributionsLoglogistic, error_check) { boost::random::mt19937 rng; EXPECT_NO_THROW(stan::math::loglogistic_rng(10.0, 2.0, rng)); - EXPECT_THROW(stan::math::loglogistic_rng(2.0, -1.0, rng), - std::domain_error); - EXPECT_THROW(stan::math::loglogistic_rng(-2.0, 1.0, rng), - std::domain_error); - EXPECT_THROW(stan::math::loglogistic_rng( - 10, stan::math::positive_infinity(), rng), - std::domain_error); - EXPECT_THROW(stan::math::loglogistic_rng( - stan::math::positive_infinity(), 2, rng), - std::domain_error); + EXPECT_THROW(stan::math::loglogistic_rng(2.0, -1.0, rng), std::domain_error); + EXPECT_THROW(stan::math::loglogistic_rng(-2.0, 1.0, rng), std::domain_error); + EXPECT_THROW( + stan::math::loglogistic_rng(10, stan::math::positive_infinity(), rng), + std::domain_error); + EXPECT_THROW( + stan::math::loglogistic_rng(stan::math::positive_infinity(), 2, rng), + std::domain_error); } TEST(ProbDistributionsLoglogistic, test_sampling_icdf) { @@ -69,8 +63,7 @@ TEST(ProbDistributionsLoglogistic, chiSquareGoodnessFitTest) { std::vector samples; for (int i = 0; i < N; ++i) { - samples.push_back( - stan::math::loglogistic_rng(1.2, 2.0, rng)); + samples.push_back(stan::math::loglogistic_rng(1.2, 2.0, rng)); } std::vector quantiles; for (int i = 1; i < K; ++i) { From 1eb10bc95343cbabd962b440d396d13f37a7ea5b Mon Sep 17 00:00:00 2001 From: gregorp90 Date: Mon, 17 May 2021 17:24:35 +0200 Subject: [PATCH 13/16] Update stan/math/prim/prob/loglogistic_cdf.hpp MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Tadej Ciglarič --- stan/math/prim/prob/loglogistic_cdf.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/stan/math/prim/prob/loglogistic_cdf.hpp b/stan/math/prim/prob/loglogistic_cdf.hpp index 9991c932c28..97b6948c218 100644 --- a/stan/math/prim/prob/loglogistic_cdf.hpp +++ b/stan/math/prim/prob/loglogistic_cdf.hpp @@ -86,8 +86,7 @@ return_type_t loglogistic_cdf(const T_y& y, = to_ref_if::value>( 1 / (1 + alpha_div_y_pow_beta)); - T_partials_return cdf(1.0); - cdf *= prod(prod_all); + T_partials_return cdf = prod(prod_all); if (!is_constant_all::value) { const auto& prod_all_sq = to_ref_if::value From a894707dd691dfecf63bb54844c624bd6b7ff642 Mon Sep 17 00:00:00 2001 From: gregorp90 Date: Mon, 17 May 2021 17:24:46 +0200 Subject: [PATCH 14/16] Update stan/math/prim/prob/loglogistic_cdf.hpp MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Tadej Ciglarič --- stan/math/prim/prob/loglogistic_cdf.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/prob/loglogistic_cdf.hpp b/stan/math/prim/prob/loglogistic_cdf.hpp index 97b6948c218..e3d149d1988 100644 --- a/stan/math/prim/prob/loglogistic_cdf.hpp +++ b/stan/math/prim/prob/loglogistic_cdf.hpp @@ -92,7 +92,7 @@ return_type_t loglogistic_cdf(const T_y& y, const auto& prod_all_sq = to_ref_if::value + !is_constant_all::value + !is_constant_all::value - >= 2>(prod_all * prod_all); + >= 2>(square(prod_all)); const auto& cdf_div_elt = to_ref_if::value + !is_constant_all::value + !is_constant_all::value From 2ea763ae60f01ad1da14d09e544afe87edb80a02 Mon Sep 17 00:00:00 2001 From: gregorp90 Date: Tue, 18 May 2021 11:47:31 +0200 Subject: [PATCH 15/16] Update stan/math/prim/prob/loglogistic_cdf.hpp MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Tadej Ciglarič --- stan/math/prim/prob/loglogistic_cdf.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/prob/loglogistic_cdf.hpp b/stan/math/prim/prob/loglogistic_cdf.hpp index e3d149d1988..d658e66a298 100644 --- a/stan/math/prim/prob/loglogistic_cdf.hpp +++ b/stan/math/prim/prob/loglogistic_cdf.hpp @@ -102,7 +102,7 @@ return_type_t loglogistic_cdf(const T_y& y, !is_constant_all::value + !is_constant_all::value == 2>( alpha_div_y_pow_beta * beta_val); if (!is_constant_all::value) { - const auto& y_deriv = alpha_div_times_beta * inv(y_val) * prod_all_sq; + const auto& y_deriv = alpha_div_times_beta / y_val * prod_all_sq; ops_partials.edge1_.partials_ = y_deriv * cdf_div_elt; } if (!is_constant_all::value) { From b09414ed0482dbadc48de20f514cecf15153464f Mon Sep 17 00:00:00 2001 From: gregorp90 Date: Tue, 18 May 2021 11:47:37 +0200 Subject: [PATCH 16/16] Update stan/math/prim/prob/loglogistic_cdf.hpp MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Tadej Ciglarič --- stan/math/prim/prob/loglogistic_cdf.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/prob/loglogistic_cdf.hpp b/stan/math/prim/prob/loglogistic_cdf.hpp index d658e66a298..33d4ebb161c 100644 --- a/stan/math/prim/prob/loglogistic_cdf.hpp +++ b/stan/math/prim/prob/loglogistic_cdf.hpp @@ -107,7 +107,7 @@ return_type_t loglogistic_cdf(const T_y& y, } if (!is_constant_all::value) { const auto& alpha_deriv - = -alpha_div_times_beta * inv(alpha_val) * prod_all_sq; + = -alpha_div_times_beta / alpha_val * prod_all_sq; ops_partials.edge2_.partials_ = alpha_deriv * cdf_div_elt; } }