From c2f2bad463dc99f0e6669a381c015db329fc1df0 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Thu, 22 Oct 2026 08:53:44 +0300 Subject: [PATCH] Add poisson_log_cdf, poisson_log_lcdf and poisson_log_lccdf --- stan/math/prim/prob.hpp | 3 + stan/math/prim/prob/poisson_log_cdf.hpp | 94 ++++++++++++++++++ stan/math/prim/prob/poisson_log_lccdf.hpp | 89 +++++++++++++++++ stan/math/prim/prob/poisson_log_lcdf.hpp | 87 ++++++++++++++++ .../poisson/poisson_log_ccdf_log_test.hpp | 73 ++++++++++++++ .../prob/poisson/poisson_log_cdf_log_test.hpp | 69 +++++++++++++ test/prob/poisson/poisson_log_cdf_test.hpp | 67 +++++++++++++ .../math/mix/prob/poisson_log_cdfs_test.cpp | 27 +++++ .../math/prim/prob/poisson_log_cdfs_test.cpp | 62 ++++++++++++ .../rev/prob/poisson_log_cdfs_reference.hpp | 99 +++++++++++++++++++ .../math/rev/prob/poisson_log_cdfs_test.cpp | 54 ++++++++++ 11 files changed, 724 insertions(+) create mode 100644 stan/math/prim/prob/poisson_log_cdf.hpp create mode 100644 stan/math/prim/prob/poisson_log_lccdf.hpp create mode 100644 stan/math/prim/prob/poisson_log_lcdf.hpp create mode 100644 test/prob/poisson/poisson_log_ccdf_log_test.hpp create mode 100644 test/prob/poisson/poisson_log_cdf_log_test.hpp create mode 100644 test/prob/poisson/poisson_log_cdf_test.hpp create mode 100644 test/unit/math/mix/prob/poisson_log_cdfs_test.cpp create mode 100644 test/unit/math/prim/prob/poisson_log_cdfs_test.cpp create mode 100644 test/unit/math/rev/prob/poisson_log_cdfs_reference.hpp create mode 100644 test/unit/math/rev/prob/poisson_log_cdfs_test.cpp diff --git a/stan/math/prim/prob.hpp b/stan/math/prim/prob.hpp index 4fae64c1ace..bcebd40c900 100644 --- a/stan/math/prim/prob.hpp +++ b/stan/math/prim/prob.hpp @@ -240,7 +240,10 @@ #include #include #include +#include #include +#include +#include #include #include #include diff --git a/stan/math/prim/prob/poisson_log_cdf.hpp b/stan/math/prim/prob/poisson_log_cdf.hpp new file mode 100644 index 00000000000..f03b604e688 --- /dev/null +++ b/stan/math/prim/prob/poisson_log_cdf.hpp @@ -0,0 +1,94 @@ +#ifndef STAN_MATH_PRIM_PROB_POISSON_LOG_CDF_HPP +#define STAN_MATH_PRIM_PROB_POISSON_LOG_CDF_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** \ingroup prob_dists + * The cumulative distribution function of the Poisson distribution with + * log rate alpha, Q(n + 1, exp(alpha)) with the regularized upper + * incomplete gamma function Q. It is computed from log Q, so its gradient + * stays finite where the probability underflows. + * + * @tparam T_n type of n + * @tparam T_log_rate type of alpha + * @param n random variable + * @param alpha log rate parameter + * @return probability, or product of probabilities for containers + * @throw std::domain_error if alpha is NaN + * @throw std::invalid_argument if container sizes mismatch + */ +template +inline return_type_t poisson_log_cdf(const T_n& n, + const T_log_rate& alpha) { + using T_partials_return = partials_return_t; + using T_n_ref = ref_type_if_not_constant_t; + using T_alpha_ref = ref_type_if_not_constant_t; + static constexpr const char* function = "poisson_log_cdf"; + check_consistent_sizes(function, "Random variable", n, "Log rate parameter", + alpha); + T_n_ref n_ref = n; + T_alpha_ref alpha_ref = alpha; + check_not_nan(function, "Log rate parameter", alpha_ref); + + if (size_zero(n, alpha)) { + return 1.0; + } + + auto ops_partials = make_partials_propagator(alpha_ref); + scalar_seq_view n_vec(n_ref); + scalar_seq_view alpha_vec(alpha_ref); + const size_t N = max_size(n, alpha); + + // An element with probability 0 (n < 0, or a rate exp(alpha) that is + // infinite) makes the result 0; the gradients are treated as zero + for (size_t i = 0; i < N; i++) { + if (n_vec.val(i) < 0 || std::exp(value_of_rec(alpha_vec.val(i))) == INFTY) { + return ops_partials.build(0.0); + } + } + + T_partials_return log_P(0.0); + for (size_t i = 0; i < N; i++) { + const T_partials_return lambda = exp(T_partials_return(alpha_vec.val(i))); + // An element with probability 1 (n = INT_MAX, which stands for an + // infinite count, or a rate of 0) contributes nothing + if (n_vec.val(i) == std::numeric_limits::max() + || value_of_rec(lambda) == 0) { + continue; + } + const auto tail + = internal::log_gamma_tail>( + n_vec.val(i) + 1.0, lambda); + log_P += tail.log_prob; + if constexpr (is_autodiff_v) { + partials<0>(ops_partials)[i] += tail.d_z * lambda; + } + } + const T_partials_return P = exp(log_P); + if constexpr (is_autodiff_v) { + for (size_t i = 0; i < stan::math::size(alpha); ++i) { + partials<0>(ops_partials)[i] *= P; + } + } + return ops_partials.build(P); +} + +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/prim/prob/poisson_log_lccdf.hpp b/stan/math/prim/prob/poisson_log_lccdf.hpp new file mode 100644 index 00000000000..472ff00cf1b --- /dev/null +++ b/stan/math/prim/prob/poisson_log_lccdf.hpp @@ -0,0 +1,89 @@ +#ifndef STAN_MATH_PRIM_PROB_POISSON_LOG_LCCDF_HPP +#define STAN_MATH_PRIM_PROB_POISSON_LOG_LCCDF_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** \ingroup prob_dists + * The log complementary cumulative distribution function of the Poisson + * distribution with log rate alpha, log P(n + 1, exp(alpha)) with the + * regularized lower incomplete gamma function P. It is computed in log + * space, so it stays finite where the probability underflows. + * + * @tparam T_n type of n + * @tparam T_log_rate type of alpha + * @param n random variable + * @param alpha log rate parameter + * @return log probability, or sum of log probabilities for containers + * @throw std::domain_error if alpha is NaN + * @throw std::invalid_argument if container sizes mismatch + */ +template +inline return_type_t poisson_log_lccdf(const T_n& n, + const T_log_rate& alpha) { + using T_partials_return = partials_return_t; + using T_n_ref = ref_type_if_not_constant_t; + using T_alpha_ref = ref_type_if_not_constant_t; + static constexpr const char* function = "poisson_log_lccdf"; + check_consistent_sizes(function, "Random variable", n, "Log rate parameter", + alpha); + T_n_ref n_ref = n; + T_alpha_ref alpha_ref = alpha; + check_not_nan(function, "Log rate parameter", alpha_ref); + + if (size_zero(n, alpha)) { + return 0; + } + + auto ops_partials = make_partials_propagator(alpha_ref); + scalar_seq_view n_vec(n_ref); + scalar_seq_view alpha_vec(alpha_ref); + const size_t N = max_size(n, alpha); + + // An element with probability 0 (n >= 0 with a rate exp(alpha) of 0, or + // n = INT_MAX, which stands for an infinite count) makes the result -inf; + // the gradients are treated as zero + for (size_t i = 0; i < N; i++) { + if (n_vec.val(i) >= 0 + && (std::exp(value_of_rec(alpha_vec.val(i))) == 0 + || n_vec.val(i) == std::numeric_limits::max())) { + return ops_partials.build(negative_infinity()); + } + } + + T_partials_return log_P(0.0); + for (size_t i = 0; i < N; i++) { + const T_partials_return lambda = exp(T_partials_return(alpha_vec.val(i))); + // An element with probability 1 (n < 0, or an infinite rate) + // contributes nothing + if (n_vec.val(i) < 0 || value_of_rec(lambda) == INFTY) { + continue; + } + const auto tail + = internal::log_gamma_tail>( + n_vec.val(i) + 1.0, lambda); + log_P += tail.log_prob; + if constexpr (is_autodiff_v) { + partials<0>(ops_partials)[i] += tail.d_z * lambda; + } + } + return ops_partials.build(log_P); +} + +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/prim/prob/poisson_log_lcdf.hpp b/stan/math/prim/prob/poisson_log_lcdf.hpp new file mode 100644 index 00000000000..52d09aa5fcc --- /dev/null +++ b/stan/math/prim/prob/poisson_log_lcdf.hpp @@ -0,0 +1,87 @@ +#ifndef STAN_MATH_PRIM_PROB_POISSON_LOG_LCDF_HPP +#define STAN_MATH_PRIM_PROB_POISSON_LOG_LCDF_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** \ingroup prob_dists + * The log cumulative distribution function of the Poisson distribution + * with log rate alpha, log Q(n + 1, exp(alpha)) with the regularized upper + * incomplete gamma function Q. It is computed in log space, so it stays + * finite where the probability underflows. + * + * @tparam T_n type of n + * @tparam T_log_rate type of alpha + * @param n random variable + * @param alpha log rate parameter + * @return log probability, or sum of log probabilities for containers + * @throw std::domain_error if alpha is NaN + * @throw std::invalid_argument if container sizes mismatch + */ +template +inline return_type_t poisson_log_lcdf(const T_n& n, + const T_log_rate& alpha) { + using T_partials_return = partials_return_t; + using T_n_ref = ref_type_if_not_constant_t; + using T_alpha_ref = ref_type_if_not_constant_t; + static constexpr const char* function = "poisson_log_lcdf"; + check_consistent_sizes(function, "Random variable", n, "Log rate parameter", + alpha); + T_n_ref n_ref = n; + T_alpha_ref alpha_ref = alpha; + check_not_nan(function, "Log rate parameter", alpha_ref); + + if (size_zero(n, alpha)) { + return 0; + } + + auto ops_partials = make_partials_propagator(alpha_ref); + scalar_seq_view n_vec(n_ref); + scalar_seq_view alpha_vec(alpha_ref); + const size_t N = max_size(n, alpha); + + // An element with probability 0 (n < 0, or a rate exp(alpha) that is + // infinite) makes the result -inf; the gradients are treated as zero + for (size_t i = 0; i < N; i++) { + if (n_vec.val(i) < 0 || std::exp(value_of_rec(alpha_vec.val(i))) == INFTY) { + return ops_partials.build(negative_infinity()); + } + } + + T_partials_return log_P(0.0); + for (size_t i = 0; i < N; i++) { + const T_partials_return lambda = exp(T_partials_return(alpha_vec.val(i))); + // An element with probability 1 (n = INT_MAX, which stands for an + // infinite count, or a rate of 0) contributes nothing + if (n_vec.val(i) == std::numeric_limits::max() + || value_of_rec(lambda) == 0) { + continue; + } + const auto tail + = internal::log_gamma_tail>( + n_vec.val(i) + 1.0, lambda); + log_P += tail.log_prob; + if constexpr (is_autodiff_v) { + partials<0>(ops_partials)[i] += tail.d_z * lambda; + } + } + return ops_partials.build(log_P); +} + +} // namespace math +} // namespace stan +#endif diff --git a/test/prob/poisson/poisson_log_ccdf_log_test.hpp b/test/prob/poisson/poisson_log_ccdf_log_test.hpp new file mode 100644 index 00000000000..1c38b0a67a5 --- /dev/null +++ b/test/prob/poisson/poisson_log_ccdf_log_test.hpp @@ -0,0 +1,73 @@ +// Arguments: Ints, Doubles +#include +#include +#include +#include +#include +#include + +using stan::math::var; +using std::numeric_limits; +using std::vector; + +class AgradCcdfLogPoissonLog : public AgradCcdfLogTest { + public: + void valid_values(vector >& parameters, + vector& ccdf_log) { + vector param(2); + + param[0] = 17; // n + param[1] = std::log(13.0); // alpha + parameters.push_back(param); + ccdf_log.push_back( + std::log(1.0 - 0.8904649795242025600572)); // expected ccdf_log + + // ccdf_log_function below loses the digits of 1 - cdf where the ccdf is + // small, so this point keeps the ccdf near 0.6 + param[0] = 40; // n + param[1] = std::log(42.0); // alpha + parameters.push_back(param); + ccdf_log.push_back(-0.5413568944089692353016477); // expected ccdf_log + + param[0] = 0.0; // n + param[1] = std::log(3.0); // alpha + parameters.push_back(param); + ccdf_log.push_back( + std::log(1.0 - 0.04978706836786394446248)); // expected ccdf_log + } + + void invalid_values(vector& index, vector& value) { + // alpha + index.push_back(1U); + value.push_back(numeric_limits::quiet_NaN()); + } + + bool has_lower_bound() { return false; } + + bool has_upper_bound() { return false; } + + template + stan::return_type_t ccdf_log(const T_n& n, + const T_log_rate& alpha, const T2&, + const T3&, const T4&, const T5&) { + return stan::math::poisson_log_lccdf(n, alpha); + } + + template + stan::return_type_t ccdf_log_function(const T_n& n, + const T_log_rate& alpha, + const T2&, const T3&, + const T4&, const T5&) { + using stan::math::exp; + using stan::math::lgamma; + using stan::math::log; + stan::return_type_t cdf(0); + for (int i = 0; i <= n; i++) { + cdf += exp(i * alpha - lgamma(i + 1)); + } + cdf *= exp(-exp(alpha)); + return log(1.0 - cdf); + } +}; diff --git a/test/prob/poisson/poisson_log_cdf_log_test.hpp b/test/prob/poisson/poisson_log_cdf_log_test.hpp new file mode 100644 index 00000000000..31aa5542f53 --- /dev/null +++ b/test/prob/poisson/poisson_log_cdf_log_test.hpp @@ -0,0 +1,69 @@ +// Arguments: Ints, Doubles +#include +#include +#include +#include +#include +#include + +using stan::math::var; +using std::numeric_limits; +using std::vector; + +class AgradCdfLogPoissonLog : public AgradCdfLogTest { + public: + void valid_values(vector >& parameters, + vector& cdf_log) { + vector param(2); + + param[0] = 17; // n + param[1] = std::log(13.0); // alpha + parameters.push_back(param); + cdf_log.push_back(std::log(0.8904649795242025600572)); // expected cdf_log + + param[0] = 82; // n + param[1] = std::log(42.0); // alpha + parameters.push_back(param); + cdf_log.push_back(std::log(0.9999999845303266798879)); // expected cdf_log + + param[0] = 0.0; // n + param[1] = std::log(3.0); // alpha + parameters.push_back(param); + cdf_log.push_back(std::log(0.04978706836786394446248)); // expected cdf_log + } + + void invalid_values(vector& index, vector& value) { + // alpha + index.push_back(1U); + value.push_back(numeric_limits::quiet_NaN()); + } + + bool has_lower_bound() { return false; } + + bool has_upper_bound() { return false; } + + template + stan::return_type_t cdf_log(const T_n& n, const T_log_rate& alpha, + const T2&, const T3&, const T4&, + const T5&) { + return stan::math::poisson_log_lcdf(n, alpha); + } + + template + stan::return_type_t cdf_log_function(const T_n& n, + const T_log_rate& alpha, + const T2&, const T3&, + const T4&, const T5&) { + using stan::math::exp; + using stan::math::lgamma; + using stan::math::log; + stan::return_type_t cdf(0); + for (int i = 0; i <= n; i++) { + cdf += exp(i * alpha - lgamma(i + 1)); + } + cdf *= exp(-exp(alpha)); + return log(cdf); + } +}; diff --git a/test/prob/poisson/poisson_log_cdf_test.hpp b/test/prob/poisson/poisson_log_cdf_test.hpp new file mode 100644 index 00000000000..7248b08fde3 --- /dev/null +++ b/test/prob/poisson/poisson_log_cdf_test.hpp @@ -0,0 +1,67 @@ +// Arguments: Ints, Doubles +#include +#include +#include +#include +#include +#include + +using stan::math::var; +using std::numeric_limits; +using std::vector; + +class AgradCdfPoissonLog : public AgradCdfTest { + public: + void valid_values(vector >& parameters, vector& cdf) { + vector param(2); + + param[0] = 17; // n + param[1] = std::log(13.0); // alpha + parameters.push_back(param); + cdf.push_back(0.8904649795242025600572); // expected cdf + + param[0] = 82; // n + param[1] = std::log(42.0); // alpha + parameters.push_back(param); + cdf.push_back(0.9999999845303266798879); // expected cdf + + param[0] = 0.0; // n + param[1] = std::log(3.0); // alpha + parameters.push_back(param); + cdf.push_back(0.04978706836786394446248); // expected cdf + } + + void invalid_values(vector& index, vector& value) { + // alpha + index.push_back(1U); + value.push_back(numeric_limits::quiet_NaN()); + } + + bool has_lower_bound() { return false; } + + bool has_upper_bound() { return false; } + + template + stan::return_type_t cdf(const T_n& n, const T_log_rate& alpha, + const T2&, const T3&, const T4&, + const T5&) { + return stan::math::poisson_log_cdf(n, alpha); + } + + template + stan::return_type_t cdf_function(const T_n& n, + const T_log_rate& alpha, + const T2&, const T3&, const T4&, + const T5&) { + using stan::math::exp; + using stan::math::lgamma; + stan::return_type_t cdf(0); + for (int i = 0; i <= n; i++) { + cdf += exp(i * alpha - lgamma(i + 1)); + } + cdf *= exp(-exp(alpha)); + return cdf; + } +}; diff --git a/test/unit/math/mix/prob/poisson_log_cdfs_test.cpp b/test/unit/math/mix/prob/poisson_log_cdfs_test.cpp new file mode 100644 index 00000000000..4a8d35fbbe7 --- /dev/null +++ b/test/unit/math/mix/prob/poisson_log_cdfs_test.cpp @@ -0,0 +1,27 @@ +#include +#include + +// Each function at two points: exp(alpha) > n + 2 (continued fraction) and +// exp(alpha) < n + 2 (power series side) of the incomplete gamma function. +// The count is bound into the functor; only alpha is differentiated. +TEST_F(AgradRev, mathMixScalFun_poisson_log_cdfs) { + auto check = [](const auto& f) { + stan::test::expect_ad(f(3), 2.0); + stan::test::expect_ad(f(5), 0.5); + }; + check([](int n) { + return [n](const auto& alpha) { + return stan::math::poisson_log_cdf(n, alpha); + }; + }); + check([](int n) { + return [n](const auto& alpha) { + return stan::math::poisson_log_lcdf(n, alpha); + }; + }); + check([](int n) { + return [n](const auto& alpha) { + return stan::math::poisson_log_lccdf(n, alpha); + }; + }); +} diff --git a/test/unit/math/prim/prob/poisson_log_cdfs_test.cpp b/test/unit/math/prim/prob/poisson_log_cdfs_test.cpp new file mode 100644 index 00000000000..f5b1ab9e0b4 --- /dev/null +++ b/test/unit/math/prim/prob/poisson_log_cdfs_test.cpp @@ -0,0 +1,62 @@ +#include +#include +#include +#include +#include + +TEST(ProbPoissonLog, cdfs_match_poisson_at_exp_alpha) { + using stan::math::poisson_cdf; + using stan::math::poisson_lccdf; + using stan::math::poisson_lcdf; + using stan::math::poisson_log_cdf; + using stan::math::poisson_log_lccdf; + using stan::math::poisson_log_lcdf; + std::vector n{0, 3, 17, 82}; + std::vector alpha{-2.0, 0.5, std::log(13.0), std::log(42.0)}; + std::vector lambda; + for (double a : alpha) { + lambda.push_back(std::exp(a)); + } + EXPECT_NEAR(poisson_lcdf(n, lambda), poisson_log_lcdf(n, alpha), 1e-13); + EXPECT_NEAR(poisson_lccdf(n, lambda), poisson_log_lccdf(n, alpha), 1e-13); + EXPECT_NEAR(poisson_cdf(n, lambda), poisson_log_cdf(n, alpha), 1e-15); +} + +TEST(ProbPoissonLog, cdfs_edge_values) { + using stan::math::poisson_log_cdf; + using stan::math::poisson_log_lccdf; + using stan::math::poisson_log_lcdf; + const double inf = std::numeric_limits::infinity(); + const int n_max = std::numeric_limits::max(); + // rate exp(alpha) of 0: P(N <= n) = 1 for n >= 0 + EXPECT_EQ(0.0, poisson_log_lcdf(3, -inf)); + EXPECT_EQ(-inf, poisson_log_lccdf(3, -inf)); + EXPECT_EQ(1.0, poisson_log_cdf(3, -inf)); + EXPECT_EQ(-inf, poisson_log_lccdf(3, -800.0)); + // infinite rate, also from overflow of exp(alpha) + EXPECT_EQ(-inf, poisson_log_lcdf(3, inf)); + EXPECT_EQ(-inf, poisson_log_lcdf(3, 800.0)); + EXPECT_EQ(0.0, poisson_log_lccdf(3, inf)); + EXPECT_EQ(0.0, poisson_log_cdf(3, 800.0)); + // a negative count has P(N <= n) = 0 + EXPECT_EQ(-inf, poisson_log_lcdf(-1, 1.0)); + EXPECT_EQ(0.0, poisson_log_lccdf(-1, 1.0)); + EXPECT_EQ(0.0, poisson_log_cdf(-1, 1.0)); + // n = INT_MAX stands for an infinite count + EXPECT_EQ(0.0, poisson_log_lcdf(n_max, 1.0)); + EXPECT_EQ(-inf, poisson_log_lccdf(n_max, 1.0)); + EXPECT_EQ(1.0, poisson_log_cdf(n_max, 1.0)); + // an element with probability 1 keeps the other elements + std::vector n{-1, 3}; + EXPECT_DOUBLE_EQ(poisson_log_lccdf(3, 1.0), poisson_log_lccdf(n, 1.0)); + std::vector m{n_max, 3}; + EXPECT_DOUBLE_EQ(poisson_log_lcdf(3, 1.0), poisson_log_lcdf(m, 1.0)); + EXPECT_DOUBLE_EQ(poisson_log_cdf(3, 1.0), poisson_log_cdf(m, 1.0)); +} + +TEST(ProbPoissonLog, cdfs_throw_for_nan_alpha) { + const double nan = std::numeric_limits::quiet_NaN(); + EXPECT_THROW(stan::math::poisson_log_cdf(3, nan), std::domain_error); + EXPECT_THROW(stan::math::poisson_log_lcdf(3, nan), std::domain_error); + EXPECT_THROW(stan::math::poisson_log_lccdf(3, nan), std::domain_error); +} diff --git a/test/unit/math/rev/prob/poisson_log_cdfs_reference.hpp b/test/unit/math/rev/prob/poisson_log_cdfs_reference.hpp new file mode 100644 index 00000000000..2d5c11f5e85 --- /dev/null +++ b/test/unit/math/rev/prob/poisson_log_cdfs_reference.hpp @@ -0,0 +1,99 @@ +#ifndef TEST_UNIT_MATH_REV_PROB_POISSON_LOG_CDFS_REFERENCE_HPP +#define TEST_UNIT_MATH_REV_PROB_POISSON_LOG_CDFS_REFERENCE_HPP + +#include +#include + +namespace stan { +namespace test { +namespace poisson_log_cdfs { + +/** + * Values and gradients of poisson_log_cdf, poisson_log_lcdf and + * poisson_log_lccdf, computed with mpmath at 60 digits by two routes: the + * closed form in the regularized incomplete gamma function with numerical + * differentiation, and a quadrature of the scaled density with the chain + * rule. Only rows where the routes agree to 1e-25 are kept. Each row is one + * point (a, z) = (n + 1, exp(alpha)) of the incomplete gamma function; the + * comment names the region. The field y holds n, p1 holds alpha, and the + * gradient is with respect to alpha. + */ +struct Row { + std::string function; + double y; + double p1; + double p2; + double value; + double grad[3]; + int n_grad; +}; + +inline const std::vector& rows() { + // clang-format off + // NOLINTBEGIN(whitespace/line_length) + static const std::vector r = { + {"poisson_log_cdf", 9.0, 6.907755278982137, 0.0, 1.4114888035657156e-413, {-1.3987981970763469e-410, 0, 0}, 1}, // upper_deep + {"poisson_log_cdf", 149.0, 7.090076835776092, 0.0, 1.3223466032723597e-323, {-1.3899733428888607e-320, 0, 0}, 1}, // upper_deep + {"poisson_log_cdf", 199.0, 7.3777589082278725, 0.0, 1.6236213327688525e-430, {-2.274923733659286e-427, 0, 0}, 1}, // upper_deep + {"poisson_log_cdf", 9.0, 5.703782474656201, 0.0, 2.8784992655371365e-114, {-8.3773168797734638e-112, 0, 0}, 1}, // upper_tail + {"poisson_log_cdf", 149.0, 6.396929655216146, 0.0, 8.1414836560592776e-108, {-3.6744831799506609e-105, 0, 0}, 1}, // upper_tail + {"poisson_log_cdf", 4.0, 2.0149030205422647, 0.0, 0.13206185628772066, {-0.54687297341277503, 0, 0}, 1}, // bulk + {"poisson_log_cdf", 49.0, 4.0943445622221, 0.0, 0.084406681093692144, {-1.1635598764140689, 0, 0}, 1}, // bulk + {"poisson_log_cdf", 199.0, 5.298317366548036, 0.0, 0.49059658199276574, {-5.6395455371841644, 0, 0}, 1}, // bulk_big + {"poisson_log_cdf", 199.0, 5.703782474656201, 0.0, 3.3711032555259333e-10, {-3.4677495862069145e-8, 0, 0}, 1}, // bulk_big + {"poisson_log_cdf", 399.0, 5.991464547107982, 0.0, 0.49335087016109675, {-7.9771835220132198, 0, 0}, 1}, // bulk_big + {"poisson_log_cdf", 999.0, 7.003065458786462, 0.0, 0.0010593232539300154, {-0.11590731420658176, 0, 0}, 1}, // bulk_big + {"poisson_log_cdf", 999.0, 6.802394763324311, 0.0, 0.99945009773428821, {-0.059269532854190015, 0, 0}, 1}, // bulk_big + {"poisson_log_cdf", 399.0, 5.703782474656201, 0.0, 0.99999997793347986, {-2.2688361353646397e-6, 0, 0}, 1}, // bulk_big + {"poisson_log_cdf", 199.0, 3.912023005428146, 0.0, 1.0, {-3.0438141221254986e-55, 0, 0}, 1}, // lower_tail + {"poisson_log_cdf", 49.0, 1.6094379124341003, 0.0, 1.0, {-9.8383651913938528e-31, 0, 0}, 1}, // lower_tail + {"poisson_log_cdf", 999.0, 4.605170185988092, 0.0, 1.0, {-9.2450143061426858e-609, 0, 0}, 1}, // lower_deep + {"poisson_log_cdf", 499.0, 3.4011973816621555, 0.0, 1.0, {-1.3942940431553544e-406, 0, 0}, 1}, // lower_deep + {"poisson_log_cdf", 99.0, -6.907755278982137, 0.0, 1.0, {-1.0704393134139705e-456, 0, 0}, 1}, // lower_deep + {"poisson_log_lcdf", 9.0, 6.907755278982137, 0.0, -950.62299837015584, {-991.00906329734278, 0, 0}, 1}, // upper_deep + {"poisson_log_lcdf", 149.0, 7.090076835776092, 0.0, -743.45557714901798, {-1051.1414627973844, 0, 0}, 1}, // upper_deep + {"poisson_log_lcdf", 199.0, 7.3777589082278725, 0.0, -989.62693094236569, {-1401.1418104366312, 0, 0}, 1}, // upper_deep + {"poisson_log_lcdf", 9.0, 5.703782474656201, 0.0, -261.43743153134223, {-291.03071103998464, 0, 0}, 1}, // upper_tail + {"poisson_log_lcdf", 149.0, 6.396929655216146, 0.0, -246.58221761262201, {-451.32844763692875, 0, 0}, 1}, // upper_tail + {"poisson_log_lcdf", 4.0, 2.0149030205422647, 0.0, -2.0244848579180105, {-4.1410365474593457, 0, 0}, 1}, // bulk + {"poisson_log_lcdf", 49.0, 4.0943445622221, 0.0, -2.4721087206353344, {-13.785163228044796, 0, 0}, 1}, // bulk + {"poisson_log_lcdf", 199.0, 5.298317366548036, 0.0, -0.71213311417597261, {-11.495280937907807, 0, 0}, 1}, // bulk_big + {"poisson_log_lcdf", 199.0, 5.703782474656201, 0.0, -21.810610863624879, {-102.86690508582192, 0, 0}, 1}, // bulk_big + {"poisson_log_lcdf", 399.0, 5.991464547107982, 0.0, -0.70653465390004304, {-16.1693918152174, 0, 0}, 1}, // bulk_big + {"poisson_log_lcdf", 999.0, 7.003065458786462, 0.0, -6.8501250144362293, {-109.41637859508294, 0, 0}, 1}, // bulk_big + {"poisson_log_lcdf", 999.0, 6.802394763324311, 0.0, -0.0005500535174143549, {-0.059302143237117671, 0, 0}, 1}, // bulk_big + {"poisson_log_lcdf", 399.0, 5.703782474656201, 0.0, -2.2066520386349792e-8, {-2.2688361854299591e-6, 0, 0}, 1}, // bulk_big + {"poisson_log_lcdf", 199.0, 3.912023005428146, 0.0, -2.0247590148473267e-57, {-3.0438141221254986e-55, 0, 0}, 1}, // lower_tail + {"poisson_log_lcdf", 49.0, 1.6094379124341003, 0.0, -2.1810592140784796e-32, {-9.8383651913938528e-31, 0, 0}, 1}, // lower_tail + {"poisson_log_lcdf", 999.0, 4.605170185988092, 0.0, -1.0270971815586292e-611, {-9.2450143061426858e-609, 0, 0}, 1}, // lower_deep + {"poisson_log_lcdf", 499.0, 3.4011973816621555, 0.0, -2.9661812029195561e-409, {-1.3942940431553544e-406, 0, 0}, 1}, // lower_deep + {"poisson_log_lcdf", 99.0, -6.907755278982137, 0.0, -1.0704499119269212e-458, {-1.0704393134139705e-456, 0, 0}, 1}, // lower_deep + {"poisson_log_lccdf", 9.0, 6.907755278982137, 0.0, -1.4114888035657156e-413, {1.3987981970763469e-410, 0, 0}, 1}, // upper_deep + {"poisson_log_lccdf", 149.0, 7.090076835776092, 0.0, -1.3223466032723597e-323, {1.3899733428888607e-320, 0, 0}, 1}, // upper_deep + {"poisson_log_lccdf", 199.0, 7.3777589082278725, 0.0, -1.6236213327688525e-430, {2.274923733659286e-427, 0, 0}, 1}, // upper_deep + {"poisson_log_lccdf", 9.0, 5.703782474656201, 0.0, -2.8784992655371365e-114, {8.3773168797734638e-112, 0, 0}, 1}, // upper_tail + {"poisson_log_lccdf", 149.0, 6.396929655216146, 0.0, -8.1414836560592776e-108, {3.6744831799506609e-105, 0, 0}, 1}, // upper_tail + {"poisson_log_lccdf", 4.0, 2.0149030205422647, 0.0, -0.14163482986540188, {0.63008288940238452, 0, 0}, 1}, // bulk + {"poisson_log_lccdf", 49.0, 4.0943445622221, 0.0, -0.08818298788719865, {1.2708260888185176, 0, 0}, 1}, // bulk + {"poisson_log_lccdf", 199.0, 5.298317366548036, 0.0, -0.6745150065903589, {11.070882797068462, 0, 0}, 1}, // bulk_big + {"poisson_log_lccdf", 199.0, 5.703782474656201, 0.0, -3.3711032560941502e-10, {3.4677495873759287e-8, 0, 0}, 1}, // bulk_big + {"poisson_log_lccdf", 399.0, 5.991464547107982, 0.0, -0.67993656656899868, {15.744986129848256, 0, 0}, 1}, // bulk_big + {"poisson_log_lccdf", 999.0, 7.003065458786462, 0.0, -0.0010598847333686817, {0.11603022772496961, 0, 0}, 1}, // bulk_big + {"poisson_log_lccdf", 999.0, 6.802394763324311, 0.0, -7.5057699942338739, {107.78193971881817, 0, 0}, 1}, // bulk_big + {"poisson_log_lccdf", 399.0, 5.703782474656201, 0.0, -17.629204302744172, {102.81803024099743, 0, 0}, 1}, // bulk_big + {"poisson_log_lccdf", 199.0, 3.912023005428146, 0.0, -130.54189961219743, {150.32969848784755, 0, 0}, 1}, // lower_tail + {"poisson_log_lccdf", 49.0, 1.6094379124341003, 0.0, -72.902912338991316, {45.108198474797784, 0, 0}, 1}, // lower_tail + {"poisson_log_lccdf", 999.0, 4.605170185988092, 0.0, -1406.8527552662529, {900.11096049482814, 0, 0}, 1}, // lower_deep + {"poisson_log_lccdf", 499.0, 3.4011973816621555, 0.0, -940.67002769935042, {470.06367708856666, 0, 0}, 1}, // lower_deep + {"poisson_log_lccdf", 99.0, -6.907755278982137, 0.0, -1054.515893552739, {99.999009901086207, 0, 0}, 1}, // lower_deep + }; + // NOLINTEND + // clang-format on + return r; +} + +} // namespace poisson_log_cdfs +} // namespace test +} // namespace stan + +#endif diff --git a/test/unit/math/rev/prob/poisson_log_cdfs_test.cpp b/test/unit/math/rev/prob/poisson_log_cdfs_test.cpp new file mode 100644 index 00000000000..d31094dc1f2 --- /dev/null +++ b/test/unit/math/rev/prob/poisson_log_cdfs_test.cpp @@ -0,0 +1,54 @@ +#include +#include +#include +#include +#include +#include + +namespace { + +/** + * Relative error; a reference below 1e-300 (zero or subnormal, with few + * significant bits) needs an absolute error below 1e-300. + */ +void expect_close(double expected, double actual, double tol, + const std::string& what) { + ASSERT_TRUE(std::isfinite(actual)) << what << ": not finite, " << actual; + if (std::fabs(expected) < 1e-300) { + EXPECT_LE(std::fabs(actual - expected), 1e-300) + << what << ": expected " << expected << ", got " << actual; + return; + } + const double err = std::fabs(actual - expected) / std::fabs(expected); + EXPECT_LT(err, tol) << what << ": relative error " << err << ", expected " + << expected << ", got " << actual; +} + +} // namespace + +TEST_F(AgradRev, ProbDistributionsPoissonLogCdfs_reference) { + using stan::math::var; + namespace ref = stan::test::poisson_log_cdfs; + for (const auto& r : ref::rows()) { + const int n = static_cast(r.y); + const std::string at = r.function + " at n = " + std::to_string(n) + + ", alpha = " + std::to_string(r.p1); + ASSERT_EQ(1, r.n_grad) << at; + var alpha = r.p1; + var lp; + if (r.function == "poisson_log_cdf") { + lp = stan::math::poisson_log_cdf(n, alpha); + } else if (r.function == "poisson_log_lcdf") { + lp = stan::math::poisson_log_lcdf(n, alpha); + } else { + ASSERT_EQ("poisson_log_lccdf", r.function); + lp = stan::math::poisson_log_lccdf(n, alpha); + } + lp.grad(); + // exp(alpha) is rounded, and near the mean of a large count the value + // changes by thousands of times the relative change of the rate + expect_close(r.value, lp.val(), 1e-11, "value of " + at); + expect_close(r.grad[0], alpha.adj(), 1e-9, "gradient of " + at); + stan::math::recover_memory(); + } +}