Compound Gamma-Poisson Distribution

stan/math/prim/prob.hpp

@@ -123,6 +123,11 @@
 #include <stan/math/prim/prob/gamma_log.hpp>
 #include <stan/math/prim/prob/gamma_lpdf.hpp>
 #include <stan/math/prim/prob/gamma_rng.hpp>
+#include <stan/math/prim/prob/gamma_poisson_cdf.hpp>
+#include <stan/math/prim/prob/gamma_poisson_lccdf.hpp>
+#include <stan/math/prim/prob/gamma_poisson_lcdf.hpp>
+#include <stan/math/prim/prob/gamma_poisson_lpmf.hpp>
+#include <stan/math/prim/prob/gamma_poisson_rng.hpp>
 #include <stan/math/prim/prob/gaussian_dlm_obs_log.hpp>
 #include <stan/math/prim/prob/gaussian_dlm_obs_lpdf.hpp>
 #include <stan/math/prim/prob/gaussian_dlm_obs_rng.hpp>

stan/math/prim/prob/gamma_poisson_cdf.hpp

@@ -0,0 +1,51 @@
+#ifndef STAN_MATH_PRIM_PROB_GAMMA_POISSON_CDF_HPP
+#define STAN_MATH_PRIM_PROB_GAMMA_POISSON_CDF_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/elt_divide.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
+#include <stan/math/prim/prob/neg_binomial_2_cdf.hpp>
+#include <cmath>
+
+namespace stan {
+namespace math {
+
+/**
+ * Returns the cumulative distribution function (CDF) for a Poisson random
+ * variate with a Gamma prior for the rate parameter
+ *
+ * @tparam T_n Type of count outcome
+ * @tparam T_shape Type of shape parameter for the Gamma prior
+ * @tparam T_inv_scale Type of rate parameter for the Gamma prior
+ * @param n Discrete count outcome
+ * @param alpha Shape parameter for the Gamma prior
+ * @param beta Rate parameter for the Gamma prior
+ * @return CDF
+ */
+template <typename T_n, typename T_shape, typename T_inv_scale>
+return_type_t<T_shape, T_inv_scale> gamma_poisson_cdf(const T_n& n,
+                                                      const T_shape& alpha,
+                                                      const T_inv_scale& beta) {
+  static const char* function = "gamma_poisson_cdf";
+  // To avoid an integer division below, the shape parameter is promoted to a
+  // double if it is an integer
+  using AlphaScalarT = scalar_type_t<T_shape>;
+  using PromotedIfIntT
+      = std::conditional_t<std::is_integral<AlphaScalarT>::value, double,
+                           AlphaScalarT>;
+
+  const auto& n_ref = to_ref(n);
+  ref_type_t<promote_scalar_t<PromotedIfIntT, T_shape>> alpha_ref = alpha;
+  const auto& beta_ref = to_ref(beta);
+
+  check_nonnegative(function, "Random variable", n_ref);
+  check_positive_finite(function, "Shape parameter", alpha_ref);
+  check_positive_finite(function, "Inverse scale parameter", beta_ref);
+
+  return neg_binomial_2_cdf(n_ref, elt_divide(alpha_ref, beta_ref), alpha_ref);
+}
+
+}  // namespace math
+}  // namespace stan
+#endif

stan/math/prim/prob/gamma_poisson_lccdf.hpp

@@ -0,0 +1,39 @@
+#ifndef STAN_MATH_PRIM_PROB_GAMMA_POISSON_LCCDF_HPP
+#define STAN_MATH_PRIM_PROB_GAMMA_POISSON_LCCDF_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/elt_divide.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
+#include <stan/math/prim/prob/neg_binomial_2_lccdf.hpp>
+#include <cmath>
+
+namespace stan {
+namespace math {
+
+template <typename T_n, typename T_shape, typename T_inv_scale>
+return_type_t<T_shape, T_inv_scale> gamma_poisson_lccdf(
+    const T_n& n, const T_shape& alpha, const T_inv_scale& beta) {
+  static const char* function = "gamma_poisson_lccdf";
+  // To avoid an integer division below, the shape parameter is promoted to a
+  // double if it is an integer
+  using AlphaScalarT = scalar_type_t<T_shape>;
+  using PromotedIfIntT
+      = std::conditional_t<std::is_integral<AlphaScalarT>::value, double,
+                           AlphaScalarT>;
+
+  const auto& n_ref = to_ref(n);
+  ref_type_t<promote_scalar_t<PromotedIfIntT, T_shape>> alpha_ref = alpha;
+  const auto& beta_ref = to_ref(beta);
+
+  check_nonnegative(function, "Random variable", n_ref);
+  check_positive_finite(function, "Shape parameter", alpha_ref);
+  check_positive_finite(function, "Inverse scale parameter", beta_ref);
+
+  return neg_binomial_2_lccdf(n_ref, elt_divide(alpha_ref, beta_ref),
+                              alpha_ref);
+}
+
+}  // namespace math
+}  // namespace stan
+#endif

stan/math/prim/prob/gamma_poisson_lcdf.hpp

@@ -0,0 +1,50 @@
+#ifndef STAN_MATH_PRIM_PROB_GAMMA_POISSON_LCDF_HPP
+#define STAN_MATH_PRIM_PROB_GAMMA_POISSON_LCDF_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/elt_divide.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
+#include <stan/math/prim/prob/neg_binomial_2_lcdf.hpp>
+#include <cmath>
+
+namespace stan {
+namespace math {
+
+/**
+ * Returns the log of the cumulative distribution function (LCDF) for a Poisson
+ * random variate with a Gamma prior for the rate parameter
+ *
+ * @tparam T_n Type of count outcome
+ * @tparam T_shape Type of shape parameter for the Gamma prior
+ * @tparam T_inv_scale Type of rate parameter for the Gamma prior
+ * @param n Discrete count outcome
+ * @param alpha Shape parameter for the Gamma prior
+ * @param beta Rate parameter for the Gamma prior
+ * @return Log CDF
+ */
+template <typename T_n, typename T_shape, typename T_inv_scale>
+return_type_t<T_shape, T_inv_scale> gamma_poisson_lcdf(
+    const T_n& n, const T_shape& alpha, const T_inv_scale& beta) {
+  static const char* function = "gamma_poisson_lcdf";
+  // To avoid an integer division below, the shape parameter is promoted to a
+  // double if it is an integer
+  using AlphaScalarT = scalar_type_t<T_shape>;
+  using PromotedIfIntT
+      = std::conditional_t<std::is_integral<AlphaScalarT>::value, double,
+                           AlphaScalarT>;
+
+  const auto& n_ref = to_ref(n);
+  ref_type_t<promote_scalar_t<PromotedIfIntT, T_shape>> alpha_ref = alpha;
+  const auto& beta_ref = to_ref(beta);
+
+  check_nonnegative(function, "Random variable", n_ref);
+  check_positive_finite(function, "Shape parameter", alpha_ref);
+  check_positive_finite(function, "Inverse scale parameter", beta_ref);
+
+  return neg_binomial_2_lcdf(n_ref, elt_divide(alpha_ref, beta_ref), alpha_ref);
+}
+
+}  // namespace math
+}  // namespace stan
+#endif

stan/math/prim/prob/gamma_poisson_lpmf.hpp

@@ -0,0 +1,64 @@
+#ifndef STAN_MATH_PRIM_PROB_GAMMA_POISSON_LPMF_HPP
+#define STAN_MATH_PRIM_PROB_GAMMA_POISSON_LPMF_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/elt_divide.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
+#include <stan/math/prim/prob/neg_binomial_2_lpmf.hpp>
+#include <cmath>
+
+namespace stan {
+namespace math {
+
+/**
+ * Returns the log-probability of a Poisson distribution with a gamma prior
+ * for the rate parameter, after marginalizing the rate parameter out of the
+ * likelihood. The likelihood is implemented using the neg_binomial_2_lpmf,
+ * for more details on the derivation see:
+ * https://gregorygundersen.com/blog/2019/09/16/poisson-gamma-nb/
+ *
+ *
+ * @tparam propto
+ * @tparam T_n Type of count outcome
+ * @tparam T_shape Type of shape parameter for the Gamma prior
+ * @tparam T_inv_scale Type of rate parameter for the Gamma prior
+ * @param n Discrete count outcome
+ * @param alpha Shape parameter for the Gamma prior
+ * @param beta Rate parameter for the Gamma prior
+ * @return log-likelihood
+ */
+template <bool propto, typename T_n, typename T_shape, typename T_inv_scale,
+          require_all_not_nonscalar_prim_or_rev_kernel_expression_t<
+              T_n, T_shape, T_inv_scale>* = nullptr>
+return_type_t<T_shape, T_inv_scale> gamma_poisson_lpmf(
+    const T_n& n, const T_shape& alpha, const T_inv_scale& beta) {
+  static const char* function = "gamma_poisson_lpmf";
+  // To avoid an integer division below, the shape parameter is promoted to a
+  // double if it is an integer
+  using AlphaScalarT = scalar_type_t<T_shape>;
+  using PromotedIfIntT
+      = std::conditional_t<std::is_integral<AlphaScalarT>::value, double,
+                           AlphaScalarT>;
+
+  const auto& n_ref = to_ref(n);
+  ref_type_t<promote_scalar_t<PromotedIfIntT, T_shape>> alpha_ref = alpha;
+  const auto& beta_ref = to_ref(beta);
+
+  check_nonnegative(function, "Random variable", n_ref);
+  check_positive_finite(function, "Shape parameter", alpha_ref);
+  check_positive_finite(function, "Inverse scale parameter", beta_ref);
+
+  return neg_binomial_2_lpmf<propto>(n_ref, elt_divide(alpha_ref, beta_ref),
+                                     alpha_ref);
+}
+
+template <typename T_n, typename T_shape, typename T_inv_scale>
+inline return_type_t<T_shape, T_inv_scale> gamma_poisson_lpmf(
+    const T_n& n, const T_shape& alpha, const T_inv_scale& beta) {
+  return gamma_poisson_lpmf<false>(n, alpha, beta);
+}
+
+}  // namespace math
+}  // namespace stan
+#endif

stan/math/prim/prob/gamma_poisson_rng.hpp

@@ -0,0 +1,57 @@
+#ifndef STAN_MATH_PRIM_PROB_GAMMA_POISSON_RNG_HPP
+#define STAN_MATH_PRIM_PROB_GAMMA_POISSON_RNG_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/elt_divide.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
+#include <stan/math/prim/prob/neg_binomial_2_rng.hpp>
+
+namespace stan {
+namespace math {
+
+/** \ingroup prob_dists
+ * Return a poisson random variate where the rate parameter follows the
+ * specified shape and inverse scale parameters using the given
+ * random number generator.
+ *
+ * alpha and beta can each be a scalar or a one-dimensional container. Any
+ * non-scalar inputs must be the same size.
+ *
+ * @tparam T_shape type of shape parameter
+ * @tparam T_inv type of inverse scale parameter
+ * @tparam RNG type of random number generator
+ * @param alpha (Sequence of) positive shape parameter(s)
+ * @param beta (Sequence of) positive inverse scale parameter(s)
+ * @param rng random number generator
+ * @return (Sequence of) poisson-gamma random variate(s)
+ * @throw std::domain_error if alpha or beta are nonpositive
+ * @throw std::invalid_argument if non-scalar arguments are of different
+ * sizes
+ */
+template <typename T_shape, typename T_inv, class RNG>
+inline typename VectorBuilder<true, int, T_shape, T_inv>::type
+gamma_poisson_rng(const T_shape& alpha, const T_inv& beta, RNG& rng) {
+  static const char* function = "gamma_poisson_rng";
+  // To avoid an integer division below, the shape parameter is promoted to a
+  // double if it is an integer
+  using AlphaScalarT = scalar_type_t<T_shape>;
+  using PromotedIfIntT
+      = std::conditional_t<std::is_integral<AlphaScalarT>::value, double,
+                           AlphaScalarT>;
+
+  check_consistent_sizes(function, "Shape parameter", alpha,
+                         "Inverse scale Parameter", beta);
+
+  ref_type_t<promote_scalar_t<PromotedIfIntT, T_shape>> alpha_ref = alpha;
+  const auto& beta_ref = to_ref(beta);
+
+  check_positive_finite(function, "Shape parameter", alpha_ref);
+  check_positive_finite(function, "Inverse scale parameter", beta_ref);
+
+  return neg_binomial_2_rng(elt_divide(alpha_ref, beta_ref), alpha_ref, rng);
+}
+
+}  // namespace math
+}  // namespace stan
+#endif

test/unit/math/prim/prob/gamma_poisson_cdf_test.cpp

@@ -0,0 +1,44 @@
+#include <stan/math/prim.hpp>
+#include <gtest/gtest.h>
+#include <vector>
+
+TEST(ProbDistributionsPoissonGammaCDF, values) {
+  using stan::math::gamma_poisson_cdf;
+
+  // Reference values calculated by extraDistr::pgpois
+  EXPECT_FLOAT_EQ(gamma_poisson_cdf(5, 1, 6), 0.999991500140248);
+  EXPECT_FLOAT_EQ(gamma_poisson_cdf(7, 8, 3), 0.982700161635877);
+  EXPECT_FLOAT_EQ(gamma_poisson_cdf(1, 6, 1), 0.0625);
+  EXPECT_FLOAT_EQ(gamma_poisson_cdf(4, 0.1, 0.1), 0.932876334310129);
+
+  std::vector<int> y{5, 7, 1, 4};
+
+  Eigen::VectorXd alpha(4);
+  alpha << 1, 8, 6, 0.1;
+
+  Eigen::VectorXd beta(4);
+  beta << 6, 3, 1, 0.1;
+
+  double cdf
+      = 0.999991500140248 * 0.982700161635877 * 0.0625 * 0.932876334310129;
+
+  EXPECT_FLOAT_EQ(gamma_poisson_cdf(y, alpha, beta), cdf);
+}
+TEST(ProbDistributionsPoissonGammaCDF, errors) {
+  EXPECT_NO_THROW(stan::math::gamma_poisson_cdf(0, 6, 5));
+  EXPECT_NO_THROW(stan::math::gamma_poisson_cdf(10, 0.5, 15));
+  EXPECT_NO_THROW(stan::math::gamma_poisson_cdf(31, 2, 1e9));
+
+  EXPECT_THROW(stan::math::gamma_poisson_cdf(4, 0.1, -2), std::domain_error);
+  EXPECT_THROW(stan::math::gamma_poisson_cdf(5, 6, -2), std::domain_error);
+  EXPECT_THROW(stan::math::gamma_poisson_cdf(6, -6, -0.1), std::domain_error);
+  EXPECT_THROW(
+      stan::math::gamma_poisson_cdf(7, stan::math::positive_infinity(), 2),
+      std::domain_error);
+  EXPECT_THROW(
+      stan::math::gamma_poisson_cdf(8, stan::math::positive_infinity(), 6),
+      std::domain_error);
+  EXPECT_THROW(
+      stan::math::gamma_poisson_cdf(9, 2, stan::math::positive_infinity()),
+      std::domain_error);
+}

test/unit/math/prim/prob/gamma_poisson_lcdf_test.cpp

@@ -0,0 +1,44 @@
+#include <stan/math/prim.hpp>
+#include <gtest/gtest.h>
+#include <vector>
+
+TEST(ProbDistributionsPoissonGammaLCDF, values) {
+  using stan::math::gamma_poisson_lcdf;
+
+  // Reference values calculated by extraDistr::pgpois
+  EXPECT_FLOAT_EQ(gamma_poisson_lcdf(5, 1, 6), -8.49989587616611e-06);
+  EXPECT_FLOAT_EQ(gamma_poisson_lcdf(7, 8, 3), -0.0174512291323641);
+  EXPECT_FLOAT_EQ(gamma_poisson_lcdf(1, 6, 1), -2.77258872223978);
+  EXPECT_FLOAT_EQ(gamma_poisson_lcdf(4, 0.1, 0.1), -0.0694826332112239);
+
+  std::vector<int> y{5, 7, 1, 4};
+
+  Eigen::VectorXd alpha(4);
+  alpha << 1, 8, 6, 0.1;
+
+  Eigen::VectorXd beta(4);
+  beta << 6, 3, 1, 0.1;
+
+  double cdf = -8.49989587616611e-06 - 0.0174512291323641 - 2.77258872223978
+               - 0.0694826332112239;
+
+  EXPECT_FLOAT_EQ(gamma_poisson_lcdf(y, alpha, beta), cdf);
+}
+TEST(ProbDistributionsPoissonGammaLCDF, errors) {
+  EXPECT_NO_THROW(stan::math::gamma_poisson_lcdf(0, 6, 5));
+  EXPECT_NO_THROW(stan::math::gamma_poisson_lcdf(10, 0.5, 15));
+  EXPECT_NO_THROW(stan::math::gamma_poisson_lcdf(31, 2, 1e9));
+
+  EXPECT_THROW(stan::math::gamma_poisson_lcdf(4, 0.1, -2), std::domain_error);
+  EXPECT_THROW(stan::math::gamma_poisson_lcdf(5, 6, -2), std::domain_error);
+  EXPECT_THROW(stan::math::gamma_poisson_lcdf(6, -6, -0.1), std::domain_error);
+  EXPECT_THROW(
+      stan::math::gamma_poisson_lcdf(7, stan::math::positive_infinity(), 2),
+      std::domain_error);
+  EXPECT_THROW(
+      stan::math::gamma_poisson_lcdf(8, stan::math::positive_infinity(), 6),
+      std::domain_error);
+  EXPECT_THROW(
+      stan::math::gamma_poisson_lcdf(9, 2, stan::math::positive_infinity()),
+      std::domain_error);
+}