diff --git a/stan/math/opencl/kernel_generator/elt_function_cl.hpp b/stan/math/opencl/kernel_generator/elt_function_cl.hpp index 96bbabbaf65..cccc7401638 100644 --- a/stan/math/opencl/kernel_generator/elt_function_cl.hpp +++ b/stan/math/opencl/kernel_generator/elt_function_cl.hpp @@ -20,6 +20,7 @@ #include #include #include +#include #include #include #include @@ -342,6 +343,11 @@ ADD_CLASSIFICATION_FUNCTION(isinf, ADD_CLASSIFICATION_FUNCTION(isnan, this->template get_arg<0>().extreme_diagonals()) +ADD_BINARY_FUNCTION_WITH_INCLUDES( + logistic_tail_deriv, opencl_kernels::log1p_exp_device_function, + opencl_kernels::log1m_inv_logit_device_function, + opencl_kernels::inv_logit_device_function, + opencl_kernels::logistic_tail_deriv_device_function) ADD_BINARY_FUNCTION_WITH_INCLUDES(fdim) ADD_BINARY_FUNCTION_WITH_INCLUDES(fmax) ADD_BINARY_FUNCTION_WITH_INCLUDES(fmin) diff --git a/stan/math/opencl/kernels/device_functions/logistic_tail_deriv.hpp b/stan/math/opencl/kernels/device_functions/logistic_tail_deriv.hpp new file mode 100644 index 00000000000..87f5ca52d90 --- /dev/null +++ b/stan/math/opencl/kernels/device_functions/logistic_tail_deriv.hpp @@ -0,0 +1,43 @@ +#ifndef STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_LOGISTIC_TAIL_DERIV_HPP +#define STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_LOGISTIC_TAIL_DERIV_HPP +#ifdef STAN_OPENCL + +#include +#include + +namespace stan { +namespace math { +namespace opencl_kernels { + +// \cond +static constexpr const char* logistic_tail_deriv_device_function + = "\n" + "#ifndef STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_LOGISTIC_TAIL_DERIV\n" + "#define " + "STAN_MATH_OPENCL_KERNELS_DEVICE_FUNCTIONS_LOGISTIC_TAIL_" + "DERIV\n" STRINGIFY( + // \endcond + /** \ingroup opencl_kernels + * + * Return inv_logit(-x) / sigma. + * + * @param x scaled difference + * @param sigma scale + * @return inv_logit(-x) / sigma + */ + double logistic_tail_deriv(double x, double sigma) { + if (x > 700.0) { + return exp(log1m_inv_logit(x) - log(sigma)); + } + return inv_logit(-x) / sigma; + } + // \cond + ) "\n#endif\n"; // NOLINT +// \endcond + +} // namespace opencl_kernels +} // namespace math +} // namespace stan + +#endif +#endif diff --git a/stan/math/opencl/prim/logistic_cdf.hpp b/stan/math/opencl/prim/logistic_cdf.hpp index 7d02536baa3..336191c827d 100644 --- a/stan/math/opencl/prim/logistic_cdf.hpp +++ b/stan/math/opencl/prim/logistic_cdf.hpp @@ -66,16 +66,16 @@ inline return_type_t logistic_cdf( auto any_y_neg_inf = colwise_max(cast(y_val == NEGATIVE_INFTY)); auto cond = y_val == INFTY; auto inv_sigma = elt_divide(1.0, sigma_val); - auto mu_minus_y_div_sigma = elt_multiply(mu_val - y_val, inv_sigma); - auto exp_scaled_diff = exp(mu_minus_y_div_sigma); - auto Pn = elt_divide(1.0, 1.0 + exp_scaled_diff); + auto scaled_diff = elt_multiply(y_val - mu_val, inv_sigma); + auto Pn = inv_logit(scaled_diff); auto P_expr = colwise_prod(select(cond, 1.0, Pn)); - auto y_deriv_tmp = select(cond, 0.0, - elt_divide(exp(mu_minus_y_div_sigma - log(sigma_val) - - 2.0 * log1p(exp_scaled_diff)), - Pn)); - auto sigma_deriv_tmp = elt_multiply(y_deriv_tmp, mu_minus_y_div_sigma); + // These are the log-scale derivatives; they are rescaled by the product P + // below. inv_logit(-scaled_diff) avoids the pdf / Pn quotient, which is + // 0 / 0 once Pn underflows. + auto deriv = logistic_tail_deriv(scaled_diff, sigma_val); + auto y_deriv_tmp = select(cond, 0.0, deriv); + auto sigma_deriv_tmp = select(cond, 0.0, elt_multiply(-deriv, scaled_diff)); matrix_cl any_y_neg_inf_cl; matrix_cl P_cl; diff --git a/stan/math/opencl/prim/logistic_lccdf.hpp b/stan/math/opencl/prim/logistic_lccdf.hpp index 2c9155969cf..cfe8b1c871b 100644 --- a/stan/math/opencl/prim/logistic_lccdf.hpp +++ b/stan/math/opencl/prim/logistic_lccdf.hpp @@ -66,16 +66,12 @@ inline return_type_t logistic_lccdf( auto any_y_neg_inf = colwise_max(cast(y_val == NEGATIVE_INFTY)); auto any_y_pos_inf = colwise_max(cast(y_val == INFTY)); auto inv_sigma = elt_divide(1.0, sigma_val); - auto mu_minus_y_div_sigma = elt_multiply(mu_val - y_val, inv_sigma); - auto exp_scaled_diff = exp(mu_minus_y_div_sigma); - auto Pn = 1.0 - elt_divide(1.0, 1.0 + exp_scaled_diff); - auto P_expr = colwise_sum(log(Pn)); - - auto mu_deriv = elt_divide( - exp(mu_minus_y_div_sigma - log(sigma_val) - 2.0 * log1p(exp_scaled_diff)), - Pn); + auto scaled_diff = elt_multiply(y_val - mu_val, inv_sigma); + auto P_expr = colwise_sum(log1m_inv_logit(scaled_diff)); + + auto mu_deriv = logistic_tail_deriv(-scaled_diff, sigma_val); auto y_deriv = -mu_deriv; - auto sigma_deriv = elt_multiply(-mu_deriv, mu_minus_y_div_sigma); + auto sigma_deriv = elt_multiply(mu_deriv, scaled_diff); matrix_cl any_y_neg_inf_cl; matrix_cl any_y_pos_inf_cl; diff --git a/stan/math/opencl/prim/logistic_lcdf.hpp b/stan/math/opencl/prim/logistic_lcdf.hpp index f16c9518bd1..841da7ae211 100644 --- a/stan/math/opencl/prim/logistic_lcdf.hpp +++ b/stan/math/opencl/prim/logistic_lcdf.hpp @@ -66,16 +66,16 @@ inline return_type_t logistic_lcdf( auto any_y_neg_inf = colwise_max(cast(y_val == NEGATIVE_INFTY)); auto cond = y_val == INFTY; auto inv_sigma = elt_divide(1.0, sigma_val); - auto mu_minus_y_div_sigma = elt_multiply(mu_val - y_val, inv_sigma); - auto exp_scaled_diff = exp(mu_minus_y_div_sigma); - auto Pn = elt_divide(1.0, 1.0 + exp_scaled_diff); - auto P_expr = colwise_sum(log(Pn)); - - auto y_deriv = elt_divide( - exp(mu_minus_y_div_sigma - log(sigma_val) - 2.0 * log1p(exp_scaled_diff)), - Pn); + auto scaled_diff = elt_multiply(y_val - mu_val, inv_sigma); + auto P_expr = colwise_sum(log_inv_logit(scaled_diff)); + + // y == INFTY contributes log(1) = 0 to P and zero to every partial; without + // the select the scale partial would be 0 * INFTY = NaN, which prim (where + // the element is skipped outright) never produces. + auto deriv = logistic_tail_deriv(scaled_diff, sigma_val); + auto y_deriv = select(cond, 0.0, deriv); auto mu_deriv = -y_deriv; - auto sigma_deriv = elt_multiply(y_deriv, mu_minus_y_div_sigma); + auto sigma_deriv = select(cond, 0.0, elt_multiply(-deriv, scaled_diff)); matrix_cl any_y_neg_inf_cl; matrix_cl P_cl; diff --git a/stan/math/opencl/prim/logistic_lpdf.hpp b/stan/math/opencl/prim/logistic_lpdf.hpp index cfcfb2b5ba7..3ef004d2df9 100644 --- a/stan/math/opencl/prim/logistic_lpdf.hpp +++ b/stan/math/opencl/prim/logistic_lpdf.hpp @@ -81,18 +81,13 @@ inline return_type_t logistic_lpdf( = colwise_sum(static_select::value>( logp1 - log(sigma_val), logp1)); - auto y_deriv = elt_multiply( - elt_divide(2.0, 1.0 + exp(y_minus_mu_div_sigma)) - 1.0, inv_sigma); - auto exp_mu_div_sigma = exp(elt_multiply(mu_val, inv_sigma)); - auto mu_deriv = elt_multiply( - 1.0 - - 2.0 - * elt_divide( - exp_mu_div_sigma, - exp_mu_div_sigma + exp(elt_multiply(y_val, inv_sigma))), - inv_sigma); + // d/dmu = tanh(z / 2) / sigma with z = (y - mu) / sigma. The y and sigma + // partials reuse it so that d/dy == -d/dmu exactly; the equivalent + // 2 / (1 + exp(z)) - 1 form loses all relative precision as z -> 0. + auto mu_deriv = elt_multiply(tanh(0.5 * y_minus_mu_div_sigma), inv_sigma); + auto y_deriv = -mu_deriv; auto sigma_deriv - = elt_multiply(-elt_multiply(y_deriv, y_minus_mu) - 1.0, inv_sigma); + = elt_multiply(elt_multiply(mu_deriv, y_minus_mu) - 1.0, inv_sigma); matrix_cl logp_cl; matrix_cl y_deriv_cl; diff --git a/stan/math/prim/prob/logistic_cdf.hpp b/stan/math/prim/prob/logistic_cdf.hpp index 3465e44eae6..3e162b24eaa 100644 --- a/stan/math/prim/prob/logistic_cdf.hpp +++ b/stan/math/prim/prob/logistic_cdf.hpp @@ -5,15 +5,16 @@ #include #include #include +#include +#include #include #include #include #include #include +#include #include -#include #include -#include namespace stan { namespace math { @@ -26,7 +27,6 @@ inline return_type_t logistic_cdf(const T_y& y, const T_loc& mu, const T_scale& sigma) { using T_partials_return = partials_return_t; - using std::exp; using T_y_ref = ref_type_t; using T_mu_ref = ref_type_t; using T_sigma_ref = ref_type_t; @@ -69,26 +69,26 @@ inline return_type_t logistic_cdf(const T_y& y, const T_partials_return y_dbl = y_vec.val(n); const T_partials_return mu_dbl = mu_vec.val(n); - const T_partials_return sigma_dbl = sigma_vec.val(n); const T_partials_return sigma_inv_vec = 1.0 / sigma_vec.val(n); - - // TODO(Andrew) Further simplify derivatives and log scale below - const T_partials_return Pn = inv_logit((y_dbl - mu_dbl) * sigma_inv_vec); + const T_partials_return scaled_diff = (y_dbl - mu_dbl) * sigma_inv_vec; + const T_partials_return Pn = inv_logit(scaled_diff); P *= Pn; - if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] - += exp(logistic_lpdf(y_dbl, mu_dbl, sigma_dbl)) / Pn; - } - if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] - += -exp(logistic_lpdf(y_dbl, mu_dbl, sigma_dbl)) / Pn; - } - if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] - += -(y_dbl - mu_dbl) * sigma_inv_vec - * exp(logistic_lpdf(y_dbl, mu_dbl, sigma_dbl)) / Pn; + if constexpr (is_any_autodiff_v) { + const T_partials_return deriv + = value_of_rec(scaled_diff) > 700.0 + ? exp(log1m_inv_logit(scaled_diff) - log(sigma_vec.val(n))) + : inv_logit(-scaled_diff) * sigma_inv_vec; + if constexpr (is_autodiff_v) { + partials<0>(ops_partials)[n] += deriv; + } + if constexpr (is_autodiff_v) { + partials<1>(ops_partials)[n] -= deriv; + } + if constexpr (is_autodiff_v) { + partials<2>(ops_partials)[n] -= scaled_diff * deriv; + } } } diff --git a/stan/math/prim/prob/logistic_lccdf.hpp b/stan/math/prim/prob/logistic_lccdf.hpp index 016da2da39d..2bc889f27ac 100644 --- a/stan/math/prim/prob/logistic_lccdf.hpp +++ b/stan/math/prim/prob/logistic_lccdf.hpp @@ -5,16 +5,17 @@ #include #include #include -#include #include +#include +#include +#include #include #include #include #include #include -#include +#include #include -#include namespace stan { namespace math { @@ -26,8 +27,6 @@ inline return_type_t logistic_lccdf(const T_y& y, const T_loc& mu, const T_scale& sigma) { using T_partials_return = partials_return_t; - using std::exp; - using std::log; using T_y_ref = ref_type_t; using T_mu_ref = ref_type_t; using T_sigma_ref = ref_type_t; @@ -54,42 +53,39 @@ inline return_type_t logistic_lccdf(const T_y& y, size_t N = max_size(y, mu, sigma); // Explicit return for extreme values - // The gradients are technically ill-defined, but treated as zero + // The gradients are technically ill-defined, but treated as zero. for (size_t i = 0; i < stan::math::size(y); i++) { if (y_vec.val(i) == NEGATIVE_INFTY) { return ops_partials.build(0.0); } } - - for (size_t n = 0; n < N; n++) { - // Explicit results for extreme values - // The gradients are technically ill-defined, but treated as zero - if (y_vec.val(n) == INFTY) { + for (size_t i = 0; i < stan::math::size(y); i++) { + if (y_vec.val(i) == INFTY) { return ops_partials.build(negative_infinity()); } + } + for (size_t n = 0; n < N; n++) { const T_partials_return y_dbl = y_vec.val(n); const T_partials_return mu_dbl = mu_vec.val(n); - const T_partials_return sigma_dbl = sigma_vec.val(n); const T_partials_return sigma_inv_vec = 1.0 / sigma_vec.val(n); + const T_partials_return scaled_diff = (y_dbl - mu_dbl) * sigma_inv_vec; + P += log1m_inv_logit(scaled_diff); - // TODO(Andrew) Further simplify derivatives and log-scale below - const T_partials_return Pn - = 1.0 - inv_logit((y_dbl - mu_dbl) * sigma_inv_vec); - P += log(Pn); - - if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] - -= exp(logistic_lpdf(y_dbl, mu_dbl, sigma_dbl)) / Pn; - } - if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] - -= -exp(logistic_lpdf(y_dbl, mu_dbl, sigma_dbl)) / Pn; - } - if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] - -= -(y_dbl - mu_dbl) * sigma_inv_vec - * exp(logistic_lpdf(y_dbl, mu_dbl, sigma_dbl)) / Pn; + if constexpr (is_any_autodiff_v) { + const T_partials_return deriv + = value_of_rec(scaled_diff) < -700.0 + ? exp(log_inv_logit(scaled_diff) - log(sigma_vec.val(n))) + : inv_logit(scaled_diff) * sigma_inv_vec; + if constexpr (is_autodiff_v) { + partials<0>(ops_partials)[n] -= deriv; + } + if constexpr (is_autodiff_v) { + partials<1>(ops_partials)[n] += deriv; + } + if constexpr (is_autodiff_v) { + partials<2>(ops_partials)[n] += scaled_diff * deriv; + } } } return ops_partials.build(P); diff --git a/stan/math/prim/prob/logistic_lcdf.hpp b/stan/math/prim/prob/logistic_lcdf.hpp index 10b9be9e78c..ee64c6bd021 100644 --- a/stan/math/prim/prob/logistic_lcdf.hpp +++ b/stan/math/prim/prob/logistic_lcdf.hpp @@ -5,16 +5,17 @@ #include #include #include -#include #include +#include +#include +#include #include #include #include #include #include -#include +#include #include -#include namespace stan { namespace math { @@ -26,8 +27,6 @@ inline return_type_t logistic_lcdf(const T_y& y, const T_loc& mu, const T_scale& sigma) { using T_partials_return = partials_return_t; - using std::exp; - using std::log; using T_y_ref = ref_type_t; using T_mu_ref = ref_type_t; using T_sigma_ref = ref_type_t; @@ -70,25 +69,25 @@ inline return_type_t logistic_lcdf(const T_y& y, const T_partials_return y_dbl = y_vec.val(n); const T_partials_return mu_dbl = mu_vec.val(n); - const T_partials_return sigma_dbl = sigma_vec.val(n); const T_partials_return sigma_inv_vec = 1.0 / sigma_vec.val(n); - // TODO(Andrew) Further simplify derivatives and log-scale below - const T_partials_return Pn = inv_logit((y_dbl - mu_dbl) * sigma_inv_vec); - P += log(Pn); + const T_partials_return scaled_diff = (y_dbl - mu_dbl) * sigma_inv_vec; + P += log_inv_logit(scaled_diff); - if constexpr (is_autodiff_v) { - partials<0>(ops_partials)[n] - += exp(logistic_lpdf(y_dbl, mu_dbl, sigma_dbl)) / Pn; - } - if constexpr (is_autodiff_v) { - partials<1>(ops_partials)[n] - += -exp(logistic_lpdf(y_dbl, mu_dbl, sigma_dbl)) / Pn; - } - if constexpr (is_autodiff_v) { - partials<2>(ops_partials)[n] - += -(y_dbl - mu_dbl) * sigma_inv_vec - * exp(logistic_lpdf(y_dbl, mu_dbl, sigma_dbl)) / Pn; + if constexpr (is_any_autodiff_v) { + const T_partials_return deriv + = value_of_rec(scaled_diff) > 700.0 + ? exp(log1m_inv_logit(scaled_diff) - log(sigma_vec.val(n))) + : inv_logit(-scaled_diff) * sigma_inv_vec; + if constexpr (is_autodiff_v) { + partials<0>(ops_partials)[n] += deriv; + } + if constexpr (is_autodiff_v) { + partials<1>(ops_partials)[n] -= deriv; + } + if constexpr (is_autodiff_v) { + partials<2>(ops_partials)[n] -= scaled_diff * deriv; + } } } return ops_partials.build(P); diff --git a/stan/math/prim/prob/logistic_lpdf.hpp b/stan/math/prim/prob/logistic_lpdf.hpp index 181ce1c5650..027daa3b5af 100644 --- a/stan/math/prim/prob/logistic_lpdf.hpp +++ b/stan/math/prim/prob/logistic_lpdf.hpp @@ -13,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -69,25 +70,20 @@ inline return_type_t logistic_lpdf(const T_y& y, logp -= sum(log(sigma_val)) * N / math::size(sigma); } - if constexpr (is_any_autodiff_v) { - const auto& exp_y_minus_mu_div_sigma = exp(y_minus_mu_div_sigma); - const auto& y_deriv - = to_ref_if<(is_autodiff_v && is_autodiff_v)>( - (2 / (1 + exp_y_minus_mu_div_sigma) - 1) * inv_sigma); + if constexpr (is_any_autodiff_v) { + // to_ref, not to_ref_if: tanh() of an Eigen argument returns a holder that + // owns its operand, so the product has to be evaluated inside this + const auto& mu_deriv = to_ref(tanh(0.5 * y_minus_mu_div_sigma) * inv_sigma); if constexpr (is_autodiff_v) { - partials<0>(ops_partials) = y_deriv; + partials<0>(ops_partials) = -mu_deriv; + } + if constexpr (is_autodiff_v) { + edge<1>(ops_partials).partials_ = mu_deriv; } if constexpr (is_autodiff_v) { - partials<2>(ops_partials) = (-y_deriv * y_minus_mu - 1) * inv_sigma; + partials<2>(ops_partials) = (mu_deriv * y_minus_mu - 1) * inv_sigma; } } - if constexpr (is_autodiff_v) { - const auto& exp_mu_div_sigma = to_ref(exp(mu_val * inv_sigma)); - edge<1>(ops_partials).partials_ - = (1 - - 2 * exp_mu_div_sigma / (exp_mu_div_sigma + exp(y_val * inv_sigma))) - * inv_sigma; - } return ops_partials.build(logp); } diff --git a/test/unit/math/opencl/rev/logistic_cdf_test.cpp b/test/unit/math/opencl/rev/logistic_cdf_test.cpp index f279c45691e..27d7c64912a 100644 --- a/test/unit/math/opencl/rev/logistic_cdf_test.cpp +++ b/test/unit/math/opencl/rev/logistic_cdf_test.cpp @@ -97,6 +97,44 @@ TEST(ProbDistributionsLogisticCdf, opencl_matches_cpu_small_y_neg_inf) { sigma.transpose().eval()); } +TEST(ProbDistributionsLogisticCdf, opencl_matches_cpu_lower_tail) { + Eigen::VectorXd y(3); + y << -700, -350, -1400; + Eigen::VectorXd mu(3); + mu << 0, 0, 0; + Eigen::VectorXd sigma(3); + sigma << 1, 0.5, 2; + + stan::math::test::compare_cpu_opencl_prim_rev(logistic_cdf_functor, y, mu, + sigma); +} + +// y == INFTY contributes a factor of 1 and zero partials; the scale partial +// used to be 0 * INFTY = NaN on the device while the CPU skipped the element. +TEST(ProbDistributionsLogisticCdf, opencl_matches_cpu_underflow_small_sigma) { + Eigen::VectorXd y(2); + y << 8e-298, 1.0; + Eigen::VectorXd mu(2); + mu << 0, 0; + Eigen::VectorXd sigma(2); + sigma << 1e-300, 1; + + stan::math::test::compare_cpu_opencl_prim_rev(logistic_cdf_functor, y, mu, + sigma); +} + +TEST(ProbDistributionsLogisticCdf, opencl_matches_cpu_y_pos_inf) { + Eigen::VectorXd y(3); + y << 0.3, INFINITY, 1.0; + Eigen::VectorXd mu(3); + mu << 0.3, 0.8, 1.0; + Eigen::VectorXd sigma(3); + sigma << 0.3, 0.8, 1.0; + + stan::math::test::compare_cpu_opencl_prim_rev(logistic_cdf_functor, y, mu, + sigma); +} + TEST(ProbDistributionsLogisticCdf, opencl_broadcast_y) { int N = 3; diff --git a/test/unit/math/opencl/rev/logistic_lccdf_test.cpp b/test/unit/math/opencl/rev/logistic_lccdf_test.cpp index 50796ef91b6..a425e65e843 100644 --- a/test/unit/math/opencl/rev/logistic_lccdf_test.cpp +++ b/test/unit/math/opencl/rev/logistic_lccdf_test.cpp @@ -79,6 +79,44 @@ TEST(ProbDistributionsLogisticLccdf, opencl_matches_cpu_small) { sigma.transpose().eval()); } +TEST(ProbDistributionsLogisticLccdf, opencl_matches_cpu_upper_tail) { + Eigen::VectorXd y(4); + y << 30, 36.5, 40, 85; + Eigen::VectorXd mu(4); + mu << 0, 0, 0, 5; + Eigen::VectorXd sigma(4); + sigma << 1, 1, 1, 2; + + stan::math::test::compare_cpu_opencl_prim_rev(logistic_lccdf_functor, y, mu, + sigma); +} + +// A single y == INFTY makes the whole result -INFTY with zero partials; the +// partials of the finite elements preceding it must not survive. +TEST(ProbDistributionsLogisticLccdf, opencl_matches_cpu_underflow_small_sigma) { + Eigen::VectorXd y(2); + y << -8e-298, -1.0; + Eigen::VectorXd mu(2); + mu << 0, 0; + Eigen::VectorXd sigma(2); + sigma << 1e-300, 1; + + stan::math::test::compare_cpu_opencl_prim_rev(logistic_lccdf_functor, y, mu, + sigma); +} + +TEST(ProbDistributionsLogisticLccdf, opencl_matches_cpu_y_pos_inf) { + Eigen::VectorXd y(3); + y << 1.5, INFINITY, 1.0; + Eigen::VectorXd mu(3); + mu << 0.3, 0.8, 1.0; + Eigen::VectorXd sigma(3); + sigma << 0.3, 0.8, 1.0; + + stan::math::test::compare_cpu_opencl_prim_rev(logistic_lccdf_functor, y, mu, + sigma); +} + TEST(ProbDistributionsLogisticLccdf, opencl_matches_cpu_small_y_neg_inf) { int N = 3; int M = 2; diff --git a/test/unit/math/opencl/rev/logistic_lcdf_test.cpp b/test/unit/math/opencl/rev/logistic_lcdf_test.cpp index c3cb1cc9d4e..67d25a4f78c 100644 --- a/test/unit/math/opencl/rev/logistic_lcdf_test.cpp +++ b/test/unit/math/opencl/rev/logistic_lcdf_test.cpp @@ -79,6 +79,42 @@ TEST(ProbDistributionsLogisticLcdf, opencl_matches_cpu_small) { sigma.transpose().eval()); } +TEST(ProbDistributionsLogisticLcdf, opencl_matches_cpu_lower_tail) { + Eigen::VectorXd y(3); + y << -745, -746, -1495; + Eigen::VectorXd mu(3); + mu << 0, 0, 5; + Eigen::VectorXd sigma(3); + sigma << 1, 1, 2; + + stan::math::test::compare_cpu_opencl_prim_rev(logistic_lcdf_functor, y, mu, + sigma); +} + +TEST(ProbDistributionsLogisticLcdf, opencl_matches_cpu_underflow_small_sigma) { + Eigen::VectorXd y(2); + y << 8e-298, 1.0; + Eigen::VectorXd mu(2); + mu << 0, 0; + Eigen::VectorXd sigma(2); + sigma << 1e-300, 1; + + stan::math::test::compare_cpu_opencl_prim_rev(logistic_lcdf_functor, y, mu, + sigma); +} + +TEST(ProbDistributionsLogisticLcdf, opencl_matches_cpu_y_pos_inf) { + Eigen::VectorXd y(2); + y << 1.5, INFINITY; + Eigen::VectorXd mu(2); + mu << 0, 0; + Eigen::VectorXd sigma(2); + sigma << 1, 2; + + stan::math::test::compare_cpu_opencl_prim_rev(logistic_lcdf_functor, y, mu, + sigma); +} + TEST(ProbDistributionsLogisticLcdf, opencl_matches_cpu_small_y_neg_inf) { int N = 3; int M = 2; diff --git a/test/unit/math/opencl/rev/logistic_lpdf_test.cpp b/test/unit/math/opencl/rev/logistic_lpdf_test.cpp index 44cc923c6ea..678eeb2381a 100644 --- a/test/unit/math/opencl/rev/logistic_lpdf_test.cpp +++ b/test/unit/math/opencl/rev/logistic_lpdf_test.cpp @@ -100,6 +100,20 @@ TEST(ProbDistributionsLogistic, opencl_matches_cpu_small) { sigma.transpose().eval()); } +TEST(ProbDistributionsLogistic, opencl_matches_cpu_large_location) { + Eigen::VectorXd y(3); + y << 711, 712, 1e-14; + Eigen::VectorXd mu(3); + mu << 710, 710, 0; + Eigen::VectorXd sigma(3); + sigma << 1, 2, 1; + + stan::math::test::compare_cpu_opencl_prim_rev(logistic_lpdf_functor, y, mu, + sigma); + stan::math::test::compare_cpu_opencl_prim_rev(logistic_lpdf_functor_propto, y, + mu, sigma); +} + TEST(ProbDistributionsLogistic, opencl_broadcast_y) { int N = 3; diff --git a/test/unit/math/rev/prob/logistic_cdf_test.cpp b/test/unit/math/rev/prob/logistic_cdf_test.cpp new file mode 100644 index 00000000000..ad1b424897f --- /dev/null +++ b/test/unit/math/rev/prob/logistic_cdf_test.cpp @@ -0,0 +1,76 @@ +#include +#include +#include +#include + +// The pdf / cdf quotient underflows to 0 / 0 in the lower tail. The partials +// are now built from inv_logit(-z) and only pick up the cdf itself as a +// factor, so they stay finite all the way down to the denormal range. +TEST_F(AgradRev, logistic_cdf_lower_tail) { + stan::math::var y = -745.0; + stan::math::var mu = 0.0; + stan::math::var sigma = 1.0; + + stan::math::var cdf = stan::math::logistic_cdf(y, mu, sigma); + cdf.grad(); + + // dF/dy = F * (1 - F) / sigma, and 1 - F == 1 to machine precision here + EXPECT_GT(cdf.val(), 0.0); + EXPECT_DOUBLE_EQ(cdf.val(), y.adj()); + EXPECT_DOUBLE_EQ(-cdf.val(), mu.adj()); + EXPECT_DOUBLE_EQ(745.0 * cdf.val(), sigma.adj()); +} + +// One ulp further out the cdf itself underflows to zero; the gradient must +// follow it to zero rather than become NaN. +TEST_F(AgradRev, logistic_cdf_lower_tail_underflow) { + stan::math::var y = -746.0; + stan::math::var mu = 0.0; + stan::math::var sigma = 1.0; + + stan::math::var cdf = stan::math::logistic_cdf(y, mu, sigma); + cdf.grad(); + + EXPECT_DOUBLE_EQ(0.0, cdf.val()); + EXPECT_DOUBLE_EQ(0.0, y.adj()); + EXPECT_DOUBLE_EQ(0.0, mu.adj()); + EXPECT_DOUBLE_EQ(0.0, sigma.adj()); +} + +// inv_logit(-z) underflows to zero above z = 745; with the cdf itself equal +// to one there, the gradient is the only thing left to get right. +TEST_F(AgradRev, logistic_cdf_underflow_rescued_by_small_sigma) { + stan::math::var y = 8e-298; + stan::math::var mu = 0.0; + stan::math::var sigma = 1e-300; + + stan::math::var cdf = stan::math::logistic_cdf(y, mu, sigma); + cdf.grad(); + + const double deriv = 3.6678745841780173e-48; + EXPECT_EQ(1.0, cdf.val()); + EXPECT_NEAR(deriv, y.adj(), 1e-10 * deriv); + EXPECT_NEAR(-deriv, mu.adj(), 1e-10 * deriv); + EXPECT_NEAR(-2.9342996673424137e-45, sigma.adj(), + 1e-10 * 2.9342996673424137e-45); +} + +TEST_F(AgradRev, logistic_cdf_location_scale) { + const double mu_value = 5.0; + const double sigma_value = 2.0; + const double scaled_diff = -1.5; + + stan::math::var y = mu_value + sigma_value * scaled_diff; + stan::math::var mu = mu_value; + stan::math::var sigma = sigma_value; + + stan::math::var cdf = stan::math::logistic_cdf(y, mu, sigma); + cdf.grad(); + + const double cdf_value = stan::math::inv_logit(scaled_diff); + const double deriv = stan::math::inv_logit(-scaled_diff) / sigma_value; + EXPECT_DOUBLE_EQ(cdf_value, cdf.val()); + EXPECT_DOUBLE_EQ(deriv * cdf_value, y.adj()); + EXPECT_DOUBLE_EQ(-deriv * cdf_value, mu.adj()); + EXPECT_DOUBLE_EQ(-scaled_diff * deriv * cdf_value, sigma.adj()); +} diff --git a/test/unit/math/rev/prob/logistic_lccdf_test.cpp b/test/unit/math/rev/prob/logistic_lccdf_test.cpp new file mode 100644 index 00000000000..b8048359a56 --- /dev/null +++ b/test/unit/math/rev/prob/logistic_lccdf_test.cpp @@ -0,0 +1,99 @@ +#include +#include +#include +#include + +TEST_F(AgradRev, logistic_lccdf_upper_tail) { + for (double y_value : {36.5, 40.0}) { + stan::math::var y = y_value; + stan::math::var mu = 0.0; + stan::math::var sigma = 1.0; + + stan::math::var log_ccdf = stan::math::logistic_lccdf(y, mu, sigma); + log_ccdf.grad(); + + const double deriv = stan::math::inv_logit(y_value); + EXPECT_DOUBLE_EQ(-stan::math::log1p_exp(y_value), log_ccdf.val()); + EXPECT_DOUBLE_EQ(-deriv, y.adj()); + EXPECT_DOUBLE_EQ(deriv, mu.adj()); + EXPECT_DOUBLE_EQ(y_value * deriv, sigma.adj()); + stan::math::recover_memory(); + } +} + +// Same tail, but with mu != 0 and sigma != 1 so that the scaled difference is +// distinguishable from y and the 1 / sigma factors in the partials are +// exercised. +TEST_F(AgradRev, logistic_lccdf_upper_tail_location_scale) { + const double mu_value = 5.0; + const double sigma_value = 2.0; + + for (double scaled_diff : {36.5, 40.0}) { + stan::math::var y = mu_value + sigma_value * scaled_diff; + stan::math::var mu = mu_value; + stan::math::var sigma = sigma_value; + + stan::math::var log_ccdf = stan::math::logistic_lccdf(y, mu, sigma); + log_ccdf.grad(); + + const double deriv = stan::math::inv_logit(scaled_diff) / sigma_value; + EXPECT_DOUBLE_EQ(-stan::math::log1p_exp(scaled_diff), log_ccdf.val()); + EXPECT_DOUBLE_EQ(-deriv, y.adj()); + EXPECT_DOUBLE_EQ(deriv, mu.adj()); + EXPECT_DOUBLE_EQ(scaled_diff * deriv, sigma.adj()); + stan::math::recover_memory(); + } +} + +// The 1 - inv_logit(z) cancellation degrades long before it returns -Inf: +// at z = 30 it gave -30.001021 for the value and -1.00102 for the y partial. +TEST_F(AgradRev, logistic_lccdf_upper_tail_moderate) { + stan::math::var y = 30.0; + stan::math::var mu = 0.0; + stan::math::var sigma = 1.0; + + stan::math::var log_ccdf = stan::math::logistic_lccdf(y, mu, sigma); + log_ccdf.grad(); + + EXPECT_NEAR(-30.000000000000092, log_ccdf.val(), 1e-12); + EXPECT_NEAR(-0.99999999999990652, y.adj(), 1e-12); + EXPECT_NEAR(0.99999999999990652, mu.adj(), 1e-12); + EXPECT_NEAR(29.999999999997197, sigma.adj(), 1e-10); +} + +// inv_logit(z) underflows to zero below z = -745, but dividing by a small +// enough sigma brings the quotient back into range. +TEST_F(AgradRev, logistic_lccdf_underflow_rescued_by_small_sigma) { + stan::math::var y = -8e-298; + stan::math::var mu = 0.0; + stan::math::var sigma = 1e-300; + + stan::math::var log_ccdf = stan::math::logistic_lccdf(y, mu, sigma); + log_ccdf.grad(); + + const double deriv = 3.6678745841780173e-48; + EXPECT_EQ(0.0, stan::math::inv_logit(-8e-298 / 1e-300)); + EXPECT_NEAR(-deriv, y.adj(), 1e-10 * deriv); + EXPECT_NEAR(deriv, mu.adj(), 1e-10 * deriv); + EXPECT_NEAR(-2.9342996673424137e-45, sigma.adj(), + 1e-10 * 2.9342996673424137e-45); +} + +// An infinite element short-circuits the result; the partials of the finite +// elements that precede it must not leak into the returned gradient. +TEST_F(AgradRev, logistic_lccdf_pos_inf_zeroes_partials) { + Eigen::Matrix y(2); + y << 1.5, stan::math::INFTY; + stan::math::var mu = 0.0; + stan::math::var sigma = 1.0; + + stan::math::var log_ccdf = stan::math::logistic_lccdf(y, mu, sigma); + log_ccdf.grad(); + + EXPECT_TRUE(std::isinf(log_ccdf.val())); + EXPECT_LT(log_ccdf.val(), 0.0); + EXPECT_DOUBLE_EQ(0.0, y(0).adj()); + EXPECT_DOUBLE_EQ(0.0, y(1).adj()); + EXPECT_DOUBLE_EQ(0.0, mu.adj()); + EXPECT_DOUBLE_EQ(0.0, sigma.adj()); +} diff --git a/test/unit/math/rev/prob/logistic_lcdf_test.cpp b/test/unit/math/rev/prob/logistic_lcdf_test.cpp new file mode 100644 index 00000000000..f0936bcc2d5 --- /dev/null +++ b/test/unit/math/rev/prob/logistic_lcdf_test.cpp @@ -0,0 +1,82 @@ +#include +#include +#include + +TEST_F(AgradRev, logistic_lcdf_lower_tail) { + for (double y_value : {-745.0, -746.0}) { + stan::math::var y = y_value; + stan::math::var mu = 0.0; + stan::math::var sigma = 1.0; + + stan::math::var log_cdf = stan::math::logistic_lcdf(y, mu, sigma); + log_cdf.grad(); + + const double deriv = stan::math::inv_logit(-y_value); + EXPECT_DOUBLE_EQ(-stan::math::log1p_exp(-y_value), log_cdf.val()); + EXPECT_DOUBLE_EQ(deriv, y.adj()); + EXPECT_DOUBLE_EQ(-deriv, mu.adj()); + EXPECT_DOUBLE_EQ(-y_value * deriv, sigma.adj()); + stan::math::recover_memory(); + } +} + +// Same tail, but with mu != 0 and sigma != 1 so that the scaled difference is +// distinguishable from y and the 1 / sigma factors in the partials are +// exercised. +TEST_F(AgradRev, logistic_lcdf_lower_tail_location_scale) { + const double mu_value = 5.0; + const double sigma_value = 2.0; + + for (double scaled_diff : {-745.0, -750.0}) { + stan::math::var y = mu_value + sigma_value * scaled_diff; + stan::math::var mu = mu_value; + stan::math::var sigma = sigma_value; + + stan::math::var log_cdf = stan::math::logistic_lcdf(y, mu, sigma); + log_cdf.grad(); + + const double deriv = stan::math::inv_logit(-scaled_diff) / sigma_value; + EXPECT_DOUBLE_EQ(-stan::math::log1p_exp(-scaled_diff), log_cdf.val()); + EXPECT_DOUBLE_EQ(deriv, y.adj()); + EXPECT_DOUBLE_EQ(-deriv, mu.adj()); + EXPECT_DOUBLE_EQ(-scaled_diff * deriv, sigma.adj()); + stan::math::recover_memory(); + } +} + +// inv_logit(-z) underflows to zero above z = 745, but dividing by a small +// enough sigma brings the quotient back into range. +TEST_F(AgradRev, logistic_lcdf_underflow_rescued_by_small_sigma) { + stan::math::var y = 8e-298; + stan::math::var mu = 0.0; + stan::math::var sigma = 1e-300; + + stan::math::var log_cdf = stan::math::logistic_lcdf(y, mu, sigma); + log_cdf.grad(); + + const double deriv = 3.6678745841780173e-48; + EXPECT_EQ(0.0, stan::math::inv_logit(-8e-298 / 1e-300)); + EXPECT_NEAR(deriv, y.adj(), 1e-10 * deriv); + EXPECT_NEAR(-deriv, mu.adj(), 1e-10 * deriv); + EXPECT_NEAR(-2.9342996673424137e-45, sigma.adj(), + 1e-10 * 2.9342996673424137e-45); +} + +// y == INFTY contributes log(1) = 0 and zero partials; it must not poison the +// gradient of the finite elements. +TEST_F(AgradRev, logistic_lcdf_pos_inf_element) { + Eigen::Matrix y(2); + y << 1.5, stan::math::INFTY; + stan::math::var mu = 0.0; + stan::math::var sigma = 1.0; + + stan::math::var log_cdf = stan::math::logistic_lcdf(y, mu, sigma); + log_cdf.grad(); + + const double deriv = stan::math::inv_logit(-1.5); + EXPECT_DOUBLE_EQ(stan::math::log_inv_logit(1.5), log_cdf.val()); + EXPECT_DOUBLE_EQ(deriv, y(0).adj()); + EXPECT_DOUBLE_EQ(0.0, y(1).adj()); + EXPECT_DOUBLE_EQ(-deriv, mu.adj()); + EXPECT_DOUBLE_EQ(-1.5 * deriv, sigma.adj()); +} diff --git a/test/unit/math/rev/prob/logistic_lpdf_test.cpp b/test/unit/math/rev/prob/logistic_lpdf_test.cpp new file mode 100644 index 00000000000..28882ac4d87 --- /dev/null +++ b/test/unit/math/rev/prob/logistic_lpdf_test.cpp @@ -0,0 +1,61 @@ +#include +#include +#include +#include + +TEST_F(AgradRev, logistic_lpdf_location_gradient_large_values) { + stan::math::var y = 711.0; + stan::math::var mu = 710.0; + stan::math::var sigma = 1.0; + + stan::math::var logp = stan::math::logistic_lpdf(y, mu, sigma); + logp.grad(); + + // z = 1: d/dmu = tanh(z / 2) / sigma, d/dy = -d/dmu, + // d/dsigma = (z * tanh(z / 2) - 1) / sigma + const double tanh_half = std::tanh(0.5); + EXPECT_DOUBLE_EQ(-1.0 - 2.0 * stan::math::log1p_exp(-1.0), logp.val()); + EXPECT_DOUBLE_EQ(tanh_half, mu.adj()); + EXPECT_DOUBLE_EQ(-tanh_half, y.adj()); + EXPECT_DOUBLE_EQ(tanh_half - 1.0, sigma.adj()); +} + +// The 2 / (1 + exp(z)) - 1 form of the y partial returns 0 for z below eps, +// while the location partial is computed from tanh(z / 2). The two must agree +// to the last bit: they are the same quantity with opposite sign. +TEST_F(AgradRev, logistic_lpdf_gradients_near_location) { + for (double scaled_diff : {1e-8, 1e-12, 1e-14, 1e-16}) { + stan::math::var y = scaled_diff; + stan::math::var mu = 0.0; + stan::math::var sigma = 1.0; + + stan::math::var logp = stan::math::logistic_lpdf(y, mu, sigma); + logp.grad(); + + EXPECT_DOUBLE_EQ(std::tanh(0.5 * scaled_diff), mu.adj()); + EXPECT_DOUBLE_EQ(-mu.adj(), y.adj()); + stan::math::recover_memory(); + } +} + +// The location partial is assigned through the container path of the edge, so +// exercise it with a vector argument and a non-unit scale. +TEST_F(AgradRev, logistic_lpdf_gradients_vectorized) { + Eigen::Matrix y(2); + y << 711.0, 712.0; + stan::math::var mu = 710.0; + stan::math::var sigma = 2.0; + + stan::math::var logp = stan::math::logistic_lpdf(y, mu, sigma); + logp.grad(); + + const double d0 = std::tanh(0.25) / 2.0; // z = 0.5 + const double d1 = std::tanh(0.5) / 2.0; // z = 1.0 + EXPECT_TRUE(std::isfinite(logp.val())); + EXPECT_DOUBLE_EQ(d0 + d1, mu.adj()); + EXPECT_DOUBLE_EQ(-d0, y(0).adj()); + EXPECT_DOUBLE_EQ(-d1, y(1).adj()); + EXPECT_DOUBLE_EQ( + (0.5 * std::tanh(0.25) - 1.0) / 2.0 + (1.0 * std::tanh(0.5) - 1.0) / 2.0, + sigma.adj()); +}