Add poisson_log_cdf, poisson_log_lcdf and poisson_log_lccdf

stan/math/prim/prob.hpp

@@ -240,7 +240,10 @@
 #include <stan/math/prim/prob/poisson_cdf_log.hpp>
 #include <stan/math/prim/prob/poisson_lccdf.hpp>
 #include <stan/math/prim/prob/poisson_lcdf.hpp>
+#include <stan/math/prim/prob/poisson_log_cdf.hpp>
 #include <stan/math/prim/prob/poisson_log_glm_lpmf.hpp>
+#include <stan/math/prim/prob/poisson_log_lccdf.hpp>
+#include <stan/math/prim/prob/poisson_log_lcdf.hpp>
 #include <stan/math/prim/prob/poisson_log_lpmf.hpp>
 #include <stan/math/prim/prob/poisson_log_rng.hpp>
 #include <stan/math/prim/prob/poisson_lpmf.hpp>

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 <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/constants.hpp>
+#include <stan/math/prim/fun/exp.hpp>
+#include <stan/math/prim/fun/log_gamma_q_dgamma.hpp>
+#include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/scalar_seq_view.hpp>
+#include <stan/math/prim/fun/size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
+#include <stan/math/prim/fun/value_of.hpp>
+#include <stan/math/prim/fun/value_of_rec.hpp>
+#include <stan/math/prim/functor/partials_propagator.hpp>
+#include <cmath>
+#include <limits>
+
+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 <typename T_n, typename T_log_rate>
+inline return_type_t<T_log_rate> poisson_log_cdf(const T_n& n,
+                                                 const T_log_rate& alpha) {
+  using T_partials_return = partials_return_t<T_n, T_log_rate>;
+  using T_n_ref = ref_type_if_not_constant_t<T_n>;
+  using T_alpha_ref = ref_type_if_not_constant_t<T_log_rate>;
+  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<T_n_ref> n_vec(n_ref);
+  scalar_seq_view<T_alpha_ref> 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<int>::max()
+        || value_of_rec(lambda) == 0) {
+      continue;
+    }
+    const auto tail
+        = internal::log_gamma_tail<true, false, is_autodiff_v<T_log_rate>>(
+            n_vec.val(i) + 1.0, lambda);
+    log_P += tail.log_prob;
+    if constexpr (is_autodiff_v<T_log_rate>) {
+      partials<0>(ops_partials)[i] += tail.d_z * lambda;
+    }
+  }
+  const T_partials_return P = exp(log_P);
+  if constexpr (is_autodiff_v<T_log_rate>) {
+    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

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 <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/constants.hpp>
+#include <stan/math/prim/fun/exp.hpp>
+#include <stan/math/prim/fun/log_gamma_q_dgamma.hpp>
+#include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/scalar_seq_view.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
+#include <stan/math/prim/fun/value_of.hpp>
+#include <stan/math/prim/fun/value_of_rec.hpp>
+#include <stan/math/prim/functor/partials_propagator.hpp>
+#include <cmath>
+#include <limits>
+
+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 <typename T_n, typename T_log_rate>
+inline return_type_t<T_log_rate> poisson_log_lccdf(const T_n& n,
+                                                   const T_log_rate& alpha) {
+  using T_partials_return = partials_return_t<T_n, T_log_rate>;
+  using T_n_ref = ref_type_if_not_constant_t<T_n>;
+  using T_alpha_ref = ref_type_if_not_constant_t<T_log_rate>;
+  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<T_n_ref> n_vec(n_ref);
+  scalar_seq_view<T_alpha_ref> 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<int>::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<false, false, is_autodiff_v<T_log_rate>>(
+            n_vec.val(i) + 1.0, lambda);
+    log_P += tail.log_prob;
+    if constexpr (is_autodiff_v<T_log_rate>) {
+      partials<0>(ops_partials)[i] += tail.d_z * lambda;
+    }
+  }
+  return ops_partials.build(log_P);
+}
+
+}  // namespace math
+}  // namespace stan
+#endif

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 <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/constants.hpp>
+#include <stan/math/prim/fun/exp.hpp>
+#include <stan/math/prim/fun/log_gamma_q_dgamma.hpp>
+#include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/scalar_seq_view.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
+#include <stan/math/prim/fun/value_of.hpp>
+#include <stan/math/prim/fun/value_of_rec.hpp>
+#include <stan/math/prim/functor/partials_propagator.hpp>
+#include <cmath>
+#include <limits>
+
+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 <typename T_n, typename T_log_rate>
+inline return_type_t<T_log_rate> poisson_log_lcdf(const T_n& n,
+                                                  const T_log_rate& alpha) {
+  using T_partials_return = partials_return_t<T_n, T_log_rate>;
+  using T_n_ref = ref_type_if_not_constant_t<T_n>;
+  using T_alpha_ref = ref_type_if_not_constant_t<T_log_rate>;
+  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<T_n_ref> n_vec(n_ref);
+  scalar_seq_view<T_alpha_ref> 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<int>::max()
+        || value_of_rec(lambda) == 0) {
+      continue;
+    }
+    const auto tail
+        = internal::log_gamma_tail<true, false, is_autodiff_v<T_log_rate>>(
+            n_vec.val(i) + 1.0, lambda);
+    log_P += tail.log_prob;
+    if constexpr (is_autodiff_v<T_log_rate>) {
+      partials<0>(ops_partials)[i] += tail.d_z * lambda;
+    }
+  }
+  return ops_partials.build(log_P);
+}
+
+}  // namespace math
+}  // namespace stan
+#endif

test/prob/poisson/poisson_log_ccdf_log_test.hpp

@@ -0,0 +1,73 @@
+// Arguments: Ints, Doubles
+#include <stan/math/prim/prob/poisson_log_lccdf.hpp>
+#include <stan/math/prim/fun/exp.hpp>
+#include <stan/math/prim/fun/log.hpp>
+#include <stan/math/prim/fun/lgamma.hpp>
+#include <cmath>
+#include <limits>
+
+using stan::math::var;
+using std::numeric_limits;
+using std::vector;
+
+class AgradCcdfLogPoissonLog : public AgradCcdfLogTest {
+ public:
+  void valid_values(vector<vector<double> >& parameters,
+                    vector<double>& ccdf_log) {
+    vector<double> 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<size_t>& index, vector<double>& value) {
+    // alpha
+    index.push_back(1U);
+    value.push_back(numeric_limits<double>::quiet_NaN());
+  }
+
+  bool has_lower_bound() { return false; }
+
+  bool has_upper_bound() { return false; }
+
+  template <typename T_n, typename T_log_rate, typename T2, typename T3,
+            typename T4, typename T5>
+  stan::return_type_t<T_log_rate> 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 <typename T_n, typename T_log_rate, typename T2, typename T3,
+            typename T4, typename T5>
+  stan::return_type_t<T_log_rate> 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<T_log_rate> 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);
+  }
+};