diff --git a/stan/math/opencl/prim/bernoulli_logit_glm_lpmf.hpp b/stan/math/opencl/prim/bernoulli_logit_glm_lpmf.hpp index 01d1bc0e23a..024a38823a8 100644 --- a/stan/math/opencl/prim/bernoulli_logit_glm_lpmf.hpp +++ b/stan/math/opencl/prim/bernoulli_logit_glm_lpmf.hpp @@ -92,21 +92,18 @@ inline return_type_t bernoulli_logit_glm_lpmf( auto ytheta_expr = matrix_vector_multiply(x_val, beta_val) + alpha_val; auto signs_expr = 2 * y_val - 1; auto ytheta_signs_expr = elt_multiply(ytheta_expr, signs_expr); - auto exp_m_ytheta_expr = exp(-ytheta_signs_expr); - const double cutoff = 20.0; - auto high_bound_expr = ytheta_signs_expr > cutoff; - auto low_bound_expr = ytheta_signs_expr < -cutoff; + auto positive_expr = ytheta_signs_expr > 0; + auto exp_m_abs_ytheta_expr + = exp(select(positive_expr, -ytheta_signs_expr, ytheta_signs_expr)); auto err_cond_expr = y_val < 0 || y_val > 1; auto logp_expr = colwise_sum(select(err_cond_expr, NOT_A_NUMBER, - select(high_bound_expr, -exp_m_ytheta_expr, - select(low_bound_expr, ytheta_signs_expr, - -log1p(exp_m_ytheta_expr))))); - auto theta_derivative_expr - = select(high_bound_expr, elt_multiply(signs_expr, exp_m_ytheta_expr), - select(low_bound_expr, signs_expr, - elt_divide(elt_multiply(signs_expr, exp_m_ytheta_expr), - (exp_m_ytheta_expr + 1)))); + select(positive_expr, 0.0, ytheta_signs_expr) + - log1p(exp_m_abs_ytheta_expr))); + auto theta_derivative_expr = elt_divide( + elt_multiply(signs_expr, + select(positive_expr, exp_m_abs_ytheta_expr, 1.0)), + exp_m_abs_ytheta_expr + 1); const int wgs = logp_expr.rows(); matrix_cl logp_cl(wgs, 1); diff --git a/stan/math/opencl/prim/bernoulli_logit_lpmf.hpp b/stan/math/opencl/prim/bernoulli_logit_lpmf.hpp index 2986e3ed501..3778e5d4ad0 100644 --- a/stan/math/opencl/prim/bernoulli_logit_lpmf.hpp +++ b/stan/math/opencl/prim/bernoulli_logit_lpmf.hpp @@ -55,18 +55,14 @@ inline return_type_t bernoulli_logit_lpmf(const T_n_cl& n, auto signs_expr = 2 * n - 1.0; // subtracting 1.0 converts int to double auto ntheta_expr = elt_multiply(signs_expr, theta_val); - auto exp_m_ntheta_expr = exp(-ntheta_expr); - static constexpr double cutoff = 20.0; - auto condition1_expr = ntheta_expr > cutoff; - auto condition2_expr = ntheta_expr < -cutoff; - auto logp_expr = colwise_sum( - select(condition1_expr, -exp_m_ntheta_expr, - select(condition2_expr, ntheta_expr, -log1p(exp_m_ntheta_expr)))); - auto deriv_expr = select( - condition1_expr, elt_multiply(signs_expr, exp_m_ntheta_expr), - select(condition2_expr, signs_expr, - elt_multiply(signs_expr, elt_divide(exp_m_ntheta_expr, - (exp_m_ntheta_expr + 1))))); + auto positive_expr = ntheta_expr > 0; + auto exp_m_abs_ntheta_expr + = exp(select(positive_expr, -ntheta_expr, ntheta_expr)); + auto logp_expr = colwise_sum(select(positive_expr, 0.0, ntheta_expr) + - log1p(exp_m_abs_ntheta_expr)); + auto deriv_expr = elt_multiply( + signs_expr, elt_divide(select(positive_expr, exp_m_abs_ntheta_expr, 1.0), + exp_m_abs_ntheta_expr + 1)); matrix_cl logp_cl; matrix_cl deriv_cl; diff --git a/stan/math/prim/prob/bernoulli_logit_glm_lpmf.hpp b/stan/math/prim/prob/bernoulli_logit_glm_lpmf.hpp index 29806c5b895..d1d7a1a73f6 100644 --- a/stan/math/prim/prob/bernoulli_logit_glm_lpmf.hpp +++ b/stan/math/prim/prob/bernoulli_logit_glm_lpmf.hpp @@ -11,8 +11,10 @@ #include #include #include +#include #include #include +#include #include #include #include @@ -112,15 +114,13 @@ inline return_type_t bernoulli_logit_glm_lpmf( ytheta = signs * (ytheta + as_array_or_scalar(alpha_val_vec)); } - // Compute the log-density and handle extreme values gracefully - // using Taylor approximations. - // And compute the derivatives wrt theta. - static constexpr double cutoff = 20.0; - Eigen::Array exp_m_ytheta = exp(-ytheta); - T_partials_return logp = sum( - (ytheta > cutoff) - .select(-exp_m_ytheta, - (ytheta < -cutoff).select(ytheta, -log1p(exp_m_ytheta)))); + // exp(-|ytheta|) in (0, 1], so no overflow enters the autodiff tape + Array exp_m_abs_ytheta + = exp((ytheta > 0).select(-ytheta, ytheta)); + // two sums: no per-element node for the difference in nested autodiff + T_partials_return logp + = sum((ytheta > 0).select(T_partials_return(0), ytheta)) + - sum(log1p(exp_m_abs_ytheta)); if (!isfinite(logp)) { check_finite(function, "Weight vector", beta); @@ -132,11 +132,8 @@ inline return_type_t bernoulli_logit_glm_lpmf( // Compute the necessary derivatives. if constexpr (is_any_autodiff_v) { Matrix theta_derivative - = (ytheta > cutoff) - .select(signs * exp_m_ytheta, - (ytheta < -cutoff) - .select(signs * T_partials_return(1.0), - signs * exp_m_ytheta / (exp_m_ytheta + 1))); + = signs * (ytheta > 0).select(exp_m_abs_ytheta, T_partials_return(1)) + / (1 + exp_m_abs_ytheta); if constexpr (is_autodiff_v) { if constexpr (T_x_rows == 1) { edge<2>(ops_partials).partials_ = theta_derivative.sum() * x_val; diff --git a/stan/math/prim/prob/bernoulli_logit_lpmf.hpp b/stan/math/prim/prob/bernoulli_logit_lpmf.hpp index f3eca1096a6..5fd408e1bd4 100644 --- a/stan/math/prim/prob/bernoulli_logit_lpmf.hpp +++ b/stan/math/prim/prob/bernoulli_logit_lpmf.hpp @@ -8,8 +8,8 @@ #include #include #include -#include #include +#include #include #include #include @@ -70,23 +70,18 @@ inline return_type_t bernoulli_logit_lpmf(const T_n& n, T_partials_return ntheta_s = signs * theta_val; ntheta = T_partials_array::Constant(1, 1, ntheta_s); } - T_partials_array exp_m_ntheta = exp(-ntheta); - static constexpr double cutoff = 20.0; - T_partials_return logp = sum( - (ntheta > cutoff) - .select(-exp_m_ntheta, - (ntheta < -cutoff).select(ntheta, -log1p(exp_m_ntheta)))); + // exp(-|ntheta|) in (0, 1], so no overflow enters the autodiff tape + T_partials_array exp_m_abs_ntheta = exp((ntheta > 0).select(-ntheta, ntheta)); + // two sums: no per-element node for the difference in nested autodiff + T_partials_return logp + = sum((ntheta > 0).select(T_partials_return(0), ntheta)) + - sum(log1p(exp_m_abs_ntheta)); auto ops_partials = make_partials_propagator(theta_ref); if constexpr (is_autodiff_v) { edge<0>(ops_partials).partials_ - = (ntheta > cutoff) - .select( - promote_scalar(signs * exp_m_ntheta), - (ntheta >= -cutoff) - .select(promote_scalar( - signs * exp_m_ntheta / (exp_m_ntheta + 1)), - promote_scalar(signs))); + = signs * (ntheta > 0).select(exp_m_abs_ntheta, T_partials_return(1)) + / (1 + exp_m_abs_ntheta); } return ops_partials.build(logp); } diff --git a/test/unit/math/mix/prob/bernoulli_logit_glm_lpmf_test.cpp b/test/unit/math/mix/prob/bernoulli_logit_glm_lpmf_test.cpp index a12cd29c636..c1baea62d27 100644 --- a/test/unit/math/mix/prob/bernoulli_logit_glm_lpmf_test.cpp +++ b/test/unit/math/mix/prob/bernoulli_logit_glm_lpmf_test.cpp @@ -1,5 +1,8 @@ #include #include +#include +#include +#include TEST_F(AgradRev, mathMixScalFun_bernoulli_logit_glm_lpmf) { auto f = [](const auto y) { @@ -37,3 +40,110 @@ TEST_F(AgradRev, mathMixScalFun_bernoulli_logit_glm_lpmf_x_column_vector) { stan::test::expect_ad(f, x, alpha, beta); } + +namespace bernoulli_logit_glm_test { +// theta = alpha + x beta = J v with v = (alpha, beta), J = [I, x]. With +// s = 2y - 1, t = s theta, p = inv_logit(t), q = 1 - p, the n-th summand is +// log(p_n) with theta_n derivatives s q, -p q and -s p q (q - p); the +// derivatives in v follow by the chain rule through J. Each tolerance is +// relative to the same sum taken over absolute values. +template +void expect_derivatives(const T_y& y, const T_x& x, + const Eigen::VectorXd& alpha, + const Eigen::VectorXd& beta) { + using stan::math::fvar; + using stan::math::var; + const int N = alpha.size(); + const int K = beta.size(); + const int M = N + K; + Eigen::MatrixXd J(N, M); + J << Eigen::MatrixXd::Identity(N, N), x.replicate(N / x.rows(), 1); + Eigen::VectorXd v(M); + v << alpha, beta; + Eigen::VectorXd theta = J * v; + Eigen::ArrayXd lp(N), d1(N), d2(N), d3(N), w(N); + for (int n = 0; n < N; ++n) { + double s = 2 * stan::get(y, n) - 1, t = s * theta(n); + double p = stan::math::inv_logit(t), q = stan::math::inv_logit(-t); + lp(n) = stan::math::log_inv_logit(t); + d1(n) = s * q; + d2(n) = -p * q; + d3(n) = -s * p * q * (q - p); + w(n) = p * q; + } + const double tiny = std::numeric_limits::min(); + Eigen::MatrixXd A = J.cwiseAbs(); + Eigen::VectorXd grad_exact = J.transpose() * d1.matrix(); + Eigen::VectorXd grad_tol + = (1e-13 * A.transpose() * d1.abs().matrix()).array() + tiny; + Eigen::MatrixXd H_exact = J.transpose() * d2.matrix().asDiagonal() * J; + Eigen::MatrixXd H_tol + = (1e-13 * A.transpose() * w.matrix().asDiagonal() * A).array() + tiny; + + auto f = [&](const auto& u) { + return stan::math::bernoulli_logit_glm_lpmf(y, x, u.head(N), u.tail(K)); + }; + double fx; + Eigen::VectorXd grad, grad_fwd; + Eigen::MatrixXd H, H_fwd, H_ffv; + std::vector grad_H; + stan::math::hessian(f, v, fx, grad, H); + stan::math::hessian(f, v, fx, grad_fwd, H_fwd); // fvar> + stan::math::grad_hessian(f, v, fx, H_ffv, grad_H); + EXPECT_NEAR(fx, lp.sum(), 1e-13 * lp.abs().sum() + tiny); + for (int i = 0; i < M; ++i) { + EXPECT_NEAR(grad(i), grad_exact(i), grad_tol(i)) << i; + EXPECT_NEAR(grad_fwd(i), grad_exact(i), grad_tol(i)) << i; + for (int j = 0; j < M; ++j) { + EXPECT_NEAR(H(i, j), H_exact(i, j), H_tol(i, j)) << i << ", " << j; + EXPECT_NEAR(H_fwd(i, j), H_exact(i, j), H_tol(i, j)) << i << ", " << j; + EXPECT_NEAR(H_ffv(i, j), H_exact(i, j), H_tol(i, j)) << i << ", " << j; + Eigen::ArrayXd Jij = J.col(i).array() * J.col(j).array(); + for (int k = 0; k < M; ++k) { + Eigen::ArrayXd Jijk = Jij * J.col(k).array(); + EXPECT_NEAR(grad_H[i](j, k), (d3 * Jijk).sum(), + 1e-13 * (w * Jijk.abs()).sum() + tiny) + << i << ", " << j << ", " << k; + } + } + } + + // derivative of the value itself + stan::math::nested_rev_autodiff nested; + Eigen::Matrix, Eigen::Dynamic, 1> v_fv(M); + for (int k = 0; k < M; ++k) { + v_fv(k) = fvar(v(k), 0); + } + f(v_fv).val_.grad(); + for (int i = 0; i < M; ++i) { + EXPECT_NEAR(v_fv(i).val_.adj(), grad_exact(i), grad_tol(i)) << i; + Eigen::Matrix>, Eigen::Dynamic, 1> v_ff(M); + for (int k = 0; k < M; ++k) { + v_ff(k) = fvar>(fvar(v(k), i == k), 0); + } + EXPECT_NEAR(f(v_ff).val_.d_, grad_exact(i), grad_tol(i)) << i; + } +} +} // namespace bernoulli_logit_glm_test + +// t = (2y - 1) theta in the tails (exp(-t) overflows at t < -709, curvature +// was 0 below -20) and at t = +-0, where the value's own derivative +// must not switch branch. Matrix and row-vector x, vector and scalar y. +TEST_F(AgradRev, mathMixScalFun_bernoulli_logit_glm_lpmf_tails) { + using bernoulli_logit_glm_test::expect_derivatives; + std::vector y{1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0}; + Eigen::MatrixXd x(13, 2); + x << 0.5, -1, 1, 0.25, -0.5, 2, 1.5, -0.75, 0.25, 0.125, -1, -0.25, 0.75, 0.5, + -0.625, 1.25, 0.375, -1, 1.125, 0.625, -0.5, -1.25, 0.75, 0.5, -0.25, 1; + Eigen::VectorXd beta(2); // dyadic x and beta: theta = 0 exactly + beta << 0.5, -0.25; + Eigen::VectorXd theta(13); + theta << -2500, 720, -30, 20.5, -0.25, 0, 0, -0.0, -0.0, 21, 21, 40, -19.5; + Eigen::VectorXd alpha = theta - x * beta; + expect_derivatives(y, x, alpha, beta); + expect_derivatives(1, x, alpha, beta); + + Eigen::RowVectorXd x_row = x.row(0); + alpha = theta.array() - x_row.dot(beta); + expect_derivatives(y, x_row, alpha, beta); +} diff --git a/test/unit/math/mix/prob/bernoulli_logit_lpmf_test.cpp b/test/unit/math/mix/prob/bernoulli_logit_lpmf_test.cpp new file mode 100644 index 00000000000..73649c25b2f --- /dev/null +++ b/test/unit/math/mix/prob/bernoulli_logit_lpmf_test.cpp @@ -0,0 +1,91 @@ +#include +#include +#include +#include +#include + +namespace bernoulli_logit_lpmf_test { +// With s = 2y - 1, t = s theta, p = inv_logit(t), q = 1 - p, the log density +// is log(p) and its theta derivatives are s q, -p q and -s p q (q - p). +struct exact { + double lp, d1, d2, d3, scale; + exact(int y, double theta) { + double s = 2 * y - 1, t = s * theta; + double p = stan::math::inv_logit(t), q = stan::math::inv_logit(-t); + lp = stan::math::log_inv_logit(t); + d1 = s * q; + d2 = -p * q; + d3 = -s * p * q * (q - p); + scale = p * q; + } +}; + +inline double tol(double x) { + return 1e-13 * std::fabs(x) + std::numeric_limits::min(); +} +} // namespace bernoulli_logit_lpmf_test + +// t = (2y - 1) theta in the tails (exp(-t) overflows at t < -709, curvature +// was 0 below -20) and at t = +-0, where the value's own derivative +// must not switch branch. Closed forms in every autodiff mode, scalar and +// vector theta. +TEST_F(AgradRev, mathMixScalFun_bernoulli_logit_lpmf_tails) { + using bernoulli_logit_lpmf_test::exact; + using bernoulli_logit_lpmf_test::tol; + using stan::math::fvar; + using stan::math::var; + std::vector y{1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0}; + Eigen::VectorXd theta(13); + theta << -2500, 720, -30, 20.5, -0.25, 0, 0, -0.0, -0.0, 21, 21, 40, -19.5; + const int N = theta.size(); + + for (int n = 0; n < N; ++n) { + exact e(y[n], theta(n)); + auto f = [&](const auto& v) { + return stan::math::bernoulli_logit_lpmf(y[n], v(0)); + }; + Eigen::VectorXd v = theta.segment(n, 1); + double fx; + Eigen::VectorXd grad; + Eigen::MatrixXd H; + stan::math::hessian(f, v, fx, grad, H); + EXPECT_NEAR(fx, e.lp, tol(e.lp)) << n; + EXPECT_NEAR(grad(0), e.d1, tol(e.d1)) << n; + EXPECT_NEAR(H(0, 0), e.d2, tol(e.d2)) << n; + stan::math::hessian(f, v, fx, grad, H); // fvar> + EXPECT_NEAR(grad(0), e.d1, tol(e.d1)) << n; + EXPECT_NEAR(H(0, 0), e.d2, tol(e.d2)) << n; + std::vector grad_H; + stan::math::grad_hessian(f, v, fx, H, grad_H); + EXPECT_NEAR(H(0, 0), e.d2, tol(e.d2)) << n; + EXPECT_NEAR(grad_H[0](0, 0), e.d3, tol(e.scale)) << n; + + // derivative of the value itself + fvar> theta_ff(fvar(theta(n), 1), 0); + EXPECT_NEAR(stan::math::bernoulli_logit_lpmf(y[n], theta_ff).val_.d_, e.d1, + tol(e.d1)) + << n; + stan::math::nested_rev_autodiff nested; + fvar theta_fv(theta(n), 0); + stan::math::bernoulli_logit_lpmf(y[n], theta_fv).val_.grad(); + EXPECT_NEAR(theta_fv.val_.adj(), e.d1, tol(e.d1)) << n; + } + + auto f + = [&](const auto& v) { return stan::math::bernoulli_logit_lpmf(y, v); }; + double fx; + Eigen::VectorXd grad; + Eigen::MatrixXd H; + stan::math::hessian(f, theta, fx, grad, H); + double lp = 0, lp_tol = 0; + for (int i = 0; i < N; ++i) { + exact e(y[i], theta(i)); + lp += e.lp; + lp_tol += tol(e.lp); + EXPECT_NEAR(grad(i), e.d1, tol(e.d1)) << i; + for (int j = 0; j < N; ++j) { + EXPECT_NEAR(H(i, j), i == j ? e.d2 : 0, tol(e.d2)) << i << ", " << j; + } + } + EXPECT_NEAR(fx, lp, lp_tol); +} diff --git a/test/unit/math/opencl/rev/bernoulli_logit_glm_lpmf_test.cpp b/test/unit/math/opencl/rev/bernoulli_logit_glm_lpmf_test.cpp index 389e6b2bc76..28733caef20 100644 --- a/test/unit/math/opencl/rev/bernoulli_logit_glm_lpmf_test.cpp +++ b/test/unit/math/opencl/rev/bernoulli_logit_glm_lpmf_test.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include using Eigen::Array; @@ -215,4 +216,33 @@ TEST(ProbDistributionsBernoulliLogitGLM, opencl_matches_cpu_big) { bernoulli_logit_glm_lpmf_functor_propto, y, x, alpha, beta); } +TEST(ProbDistributionsBernoulliLogitGLM, opencl_rev_tails_against_references) { + using bernoulli_logit_tail_refs::expect_rel; + using stan::math::var_value; + // theta = alpha + x beta with x = 1, alpha = 0.5 and beta = theta - 0.5, + // both exact for these theta, so d/dalpha = d/dbeta = d/dtheta + for (const auto& r : bernoulli_logit_tail_refs::refs()) { + vector y{r.y}; + Matrix x(1, 1); + x << 1.0; + Matrix alpha(1); + alpha << 0.5; + Matrix beta(1); + beta << r.theta - 0.5; + matrix_cl y_cl(y); + matrix_cl x_cl(x); + var_value> alpha_cl = stan::math::to_matrix_cl(alpha); + var_value> beta_cl = stan::math::to_matrix_cl(beta); + var lp + = stan::math::bernoulli_logit_glm_lpmf(y_cl, x_cl, alpha_cl, beta_cl); + lp.grad(); + Eigen::VectorXd d_alpha = stan::math::from_matrix_cl(alpha_cl.adj()); + Eigen::VectorXd d_beta = stan::math::from_matrix_cl(beta_cl.adj()); + expect_rel(r.val, lp.val(), "value", r); + expect_rel(r.d_theta, d_alpha(0), "d/dalpha", r); + expect_rel(r.d_theta, d_beta(0), "d/dbeta", r); + stan::math::recover_memory(); + } +} + #endif diff --git a/test/unit/math/opencl/rev/bernoulli_logit_lpmf_test.cpp b/test/unit/math/opencl/rev/bernoulli_logit_lpmf_test.cpp index 401b88349f7..526c85b43ab 100644 --- a/test/unit/math/opencl/rev/bernoulli_logit_lpmf_test.cpp +++ b/test/unit/math/opencl/rev/bernoulli_logit_lpmf_test.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include TEST(ProbDistributionsBernoulliLogit, error_checking) { @@ -109,4 +110,24 @@ TEST(ProbDistributionsBernoulliLogit, opencl_matches_cpu_big) { bernoulli_logit_lpmf_functor_propto, n, theta.transpose().eval()); } +TEST(ProbDistributionsBernoulliLogit, opencl_rev_tails_against_references) { + using bernoulli_logit_tail_refs::expect_rel; + using stan::math::matrix_cl; + using stan::math::var; + using stan::math::var_value; + for (const auto& r : bernoulli_logit_tail_refs::refs()) { + std::vector n{r.y}; + Eigen::VectorXd theta(1); + theta << r.theta; + matrix_cl n_cl(n); + var_value> theta_cl = stan::math::to_matrix_cl(theta); + var lp = stan::math::bernoulli_logit_lpmf(n_cl, theta_cl); + lp.grad(); + Eigen::VectorXd d_theta = stan::math::from_matrix_cl(theta_cl.adj()); + expect_rel(r.val, lp.val(), "value", r); + expect_rel(r.d_theta, d_theta(0), "d/dtheta", r); + stan::math::recover_memory(); + } +} + #endif diff --git a/test/unit/math/opencl/rev/bernoulli_logit_tail_refs.hpp b/test/unit/math/opencl/rev/bernoulli_logit_tail_refs.hpp new file mode 100644 index 00000000000..16db08cc05b --- /dev/null +++ b/test/unit/math/opencl/rev/bernoulli_logit_tail_refs.hpp @@ -0,0 +1,43 @@ +#ifndef TEST_UNIT_MATH_OPENCL_REV_BERNOULLI_LOGIT_TAIL_REFS_HPP +#define TEST_UNIT_MATH_OPENCL_REV_BERNOULLI_LOGIT_TAIL_REFS_HPP + +#include +#include +#include +#include + +namespace bernoulli_logit_tail_refs { + +// Outcome y and logit theta; the log probability log(inv_logit(t)), with +// t = (2y - 1) theta, and its derivative (2y - 1) inv_logit(-t) with +// respect to theta. mpmath, 60 digits, rounded to double. +struct ref { + int y; + double theta; + double val; + double d_theta; +}; + +inline std::vector refs() { + return {{1, 21.0, -7.582560425037146e-10, 7.582560422162384e-10}, + {0, 21.0, -21.000000000758256, -0.999999999241744}, + {1, -25.0, -25.000000000013888, 0.9999999999861121}, + {0, -25.0, -1.3887943864867583e-11, -1.3887943864771146e-11}, + {1, 30.0, -9.357622968839737e-14, 9.357622968839299e-14}, + {0, 30.0, -30.000000000000092, -0.9999999999999064}, + {1, 0.3, -0.5543552444685271, 0.425557483188341}}; +} + +// Relative check, absolute below the smallest normal double. +inline void expect_rel(double expected, double actual, const char* what, + const ref& r) { + EXPECT_NEAR( + expected, actual, + 1e-13 + * std::fmax(std::fabs(expected), std::numeric_limits::min())) + << what << " at y = " << r.y << ", theta = " << r.theta; +} + +} // namespace bernoulli_logit_tail_refs + +#endif diff --git a/test/unit/math/prim/prob/bernoulli_logit_test.cpp b/test/unit/math/prim/prob/bernoulli_logit_test.cpp index e608decd807..da5ac6f1534 100644 --- a/test/unit/math/prim/prob/bernoulli_logit_test.cpp +++ b/test/unit/math/prim/prob/bernoulli_logit_test.cpp @@ -65,20 +65,14 @@ TEST(ProbDistributionsBernoulliLogit, logitChiSquareGoodnessFitTest) { assert_chi_squared(counts, expected, 1e-6); } -TEST(ProbDistributionsBernoulliLogit, cutoff) { - double cutoff = 20; +TEST(ProbDistributionsBernoulliLogit, continuity) { for (int n = 0; n <= 1; ++n) { - for (int sign : {-1, 1}) { - double before_cutoff - = stan::math::bernoulli_logit_lpmf(n, sign * cutoff - 1e-14); - double after_cutoff - = stan::math::bernoulli_logit_lpmf(n, sign * cutoff + 1e-14); - double relative_error_at_cutoff = log(before_cutoff / after_cutoff); - EXPECT_NEAR(relative_error_at_cutoff, 0, 1e-8) - << "bernoulli_logit_lpmf changes too much around cutoff for n = " << n - << ", cutoff = " << (sign * cutoff) - << ", value at cutoff - 1e-14: " << before_cutoff - << ", value at cutoff + 1e-14: " << after_cutoff; + for (double theta : {-20.0, 0.0, 20.0}) { + double before = stan::math::bernoulli_logit_lpmf(n, theta - 1e-14); + double after = stan::math::bernoulli_logit_lpmf(n, theta + 1e-14); + EXPECT_NEAR(log(before / after), 0, 1e-8) + << "bernoulli_logit_lpmf changes too much around theta = " << theta + << " for n = " << n << ": " << before << ", " << after; } } } diff --git a/test/unit/math/rev/prob/bernoulli_logit_glm_lpmf_test.cpp b/test/unit/math/rev/prob/bernoulli_logit_glm_lpmf_test.cpp index cb8c2a16395..97b69e20bef 100644 --- a/test/unit/math/rev/prob/bernoulli_logit_glm_lpmf_test.cpp +++ b/test/unit/math/rev/prob/bernoulli_logit_glm_lpmf_test.cpp @@ -40,7 +40,7 @@ TEST_F(AgradRev, bernoulli_logit_glm_upper_tail_gradient) { std::vector{1}, x, alpha, beta); logp.grad(); - const double expected_theta_derivative = std::exp(-25.0); + const double expected_theta_derivative = stan::math::inv_logit(-25.0); EXPECT_DOUBLE_EQ(expected_theta_derivative, alpha.adj()); EXPECT_DOUBLE_EQ(expected_theta_derivative, beta.adj()(0)); EXPECT_DOUBLE_EQ(25.0 * expected_theta_derivative, x.adj()(0, 0)); @@ -57,7 +57,7 @@ TEST_F(AgradRev, bernoulli_logit_glm_upper_tail_gradient_broadcast_x) { std::vector{1, 1}, x, alpha, beta); logp.grad(); - const double expected_theta_derivative = std::exp(-25.0); + const double expected_theta_derivative = stan::math::inv_logit(-25.0); EXPECT_DOUBLE_EQ(2.0 * expected_theta_derivative, alpha.adj()); EXPECT_DOUBLE_EQ(2.0 * expected_theta_derivative, beta.adj()(0)); EXPECT_DOUBLE_EQ(50.0 * expected_theta_derivative, x.adj()(0, 0)); diff --git a/test/unit/math/rev/prob/bernoulli_logit_lpmf_test.cpp b/test/unit/math/rev/prob/bernoulli_logit_lpmf_test.cpp index 98179185526..8ed94532ac2 100644 --- a/test/unit/math/rev/prob/bernoulli_logit_lpmf_test.cpp +++ b/test/unit/math/rev/prob/bernoulli_logit_lpmf_test.cpp @@ -10,7 +10,7 @@ TEST_F(AgradRev, bernoulli_logit_lpmf_upper_tail_gradient_scalar) { stan::math::var logp = stan::math::bernoulli_logit_lpmf(1, theta); logp.grad(); - EXPECT_DOUBLE_EQ(std::exp(-25.0), theta.adj()); + EXPECT_DOUBLE_EQ(stan::math::inv_logit(-25.0), theta.adj()); } TEST_F(AgradRev, bernoulli_logit_lpmf_upper_tail_gradient_vector) { @@ -20,6 +20,6 @@ TEST_F(AgradRev, bernoulli_logit_lpmf_upper_tail_gradient_vector) { stan::math::var logp = stan::math::bernoulli_logit_lpmf(n, theta); logp.grad(); - EXPECT_DOUBLE_EQ(std::exp(-25.0), theta[0].adj()); - EXPECT_DOUBLE_EQ(-std::exp(-25.0), theta[1].adj()); + EXPECT_DOUBLE_EQ(stan::math::inv_logit(-25.0), theta[0].adj()); + EXPECT_DOUBLE_EQ(-stan::math::inv_logit(-25.0), theta[1].adj()); }