diff --git a/stan/math/prim/prob.hpp b/stan/math/prim/prob.hpp index 6aad4c64441..e58a2e7769c 100644 --- a/stan/math/prim/prob.hpp +++ b/stan/math/prim/prob.hpp @@ -75,6 +75,11 @@ #include #include #include +#include +#include +#include +#include +#include #include #include #include diff --git a/stan/math/prim/prob/discrete_range_ccdf_log.hpp b/stan/math/prim/prob/discrete_range_ccdf_log.hpp new file mode 100644 index 00000000000..324192e8c72 --- /dev/null +++ b/stan/math/prim/prob/discrete_range_ccdf_log.hpp @@ -0,0 +1,20 @@ +#ifndef STAN_MATH_PRIM_PROB_DISCRETE_RANGE_CCDF_LOG_HPP +#define STAN_MATH_PRIM_PROB_DISCRETE_RANGE_CCDF_LOG_HPP + +#include + +namespace stan { +namespace math { + +/** \ingroup prob_dists + * @deprecated use discrete_range_lccdf + */ +template +double discrete_range_ccdf_log(const T_y& y, const T_lower& lower, + const T_upper& upper) { + return discrete_range_lccdf(y, lower, upper); +} + +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/prim/prob/discrete_range_cdf.hpp b/stan/math/prim/prob/discrete_range_cdf.hpp new file mode 100644 index 00000000000..9a54ec2c00d --- /dev/null +++ b/stan/math/prim/prob/discrete_range_cdf.hpp @@ -0,0 +1,75 @@ +#ifndef STAN_MATH_PRIM_PROB_DISCRETE_RANGE_CDF_HPP +#define STAN_MATH_PRIM_PROB_DISCRETE_RANGE_CDF_HPP + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** \ingroup prob_dists + * Return the CDF of a discrete range distribution for the given y, + * lower and upper bounds (all integers). + * + * `y`, `lower` and `upper` can each be a scalar or a one-dimensional container. + * Any container arguments must be the same size. + * + * @tparam T_y type of scalar, either int or std::vector + * @tparam T_lower type of lower bound, either int or std::vector + * @tparam T_upper type of upper bound, either int or std::vector + * + * @param y integer random variable + * @param lower integer lower bound + * @param upper integer upper bound + * @return The CDF evaluated at the specified arguments. If containers are + * supplied, returns the product of the CDFs. + * @throw std::domain_error if upper is smaller than lower. + * @throw std::invalid_argument if non-scalar arguments are of different + * sizes. + */ +template +double discrete_range_cdf(const T_y& y, const T_lower& lower, + const T_upper& upper) { + static const char* function = "discrete_range_cdf"; + check_consistent_sizes(function, "Lower bound parameter", lower, + "Upper bound parameter", upper); + check_greater_or_equal(function, "Upper bound parameter", upper, lower); + + if (size_zero(y, lower, upper)) { + return 1; + } + + scalar_seq_view y_vec(y); + scalar_seq_view lower_vec(lower); + scalar_seq_view upper_vec(upper); + size_t N = max_size(y, lower, upper); + + for (size_t n = 0; n < N; ++n) { + const int y_dbl = y_vec[n]; + if (y_dbl < lower_vec[n]) { + return 0; + } + if (y_dbl > upper_vec[n]) { + return 1; + } + } + + double cdf(1.0); + for (size_t n = 0; n < N; n++) { + const double y_dbl = y_vec[n]; + const double lower_dbl = lower_vec[n]; + const double upper_dbl = upper_vec[n]; + cdf *= (y_dbl - lower_dbl + 1) / (upper_dbl - lower_dbl + 1); + } + return cdf; +} + +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/prim/prob/discrete_range_cdf_log.hpp b/stan/math/prim/prob/discrete_range_cdf_log.hpp new file mode 100644 index 00000000000..c21d7e17014 --- /dev/null +++ b/stan/math/prim/prob/discrete_range_cdf_log.hpp @@ -0,0 +1,20 @@ +#ifndef STAN_MATH_PRIM_PROB_DISCRETE_RANGE_CDF_LOG_HPP +#define STAN_MATH_PRIM_PROB_DISCRETE_RANGE_CDF_LOG_HPP + +#include + +namespace stan { +namespace math { + +/** \ingroup prob_dists + * @deprecated use discrete_range_lcdf + */ +template +double discrete_range_cdf_log(const T_y& y, const T_lower& lower, + const T_upper& upper) { + return discrete_range_lcdf(y, lower, upper); +} + +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/prim/prob/discrete_range_lccdf.hpp b/stan/math/prim/prob/discrete_range_lccdf.hpp new file mode 100644 index 00000000000..6f2d002b9b8 --- /dev/null +++ b/stan/math/prim/prob/discrete_range_lccdf.hpp @@ -0,0 +1,75 @@ +#ifndef STAN_MATH_PRIM_PROB_DISCRETE_RANGE_LCCDF_HPP +#define STAN_MATH_PRIM_PROB_DISCRETE_RANGE_LCCDF_HPP + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** \ingroup prob_dists + * Return the log CCDF of a discrete range distribution for the given y, + * lower and upper bounds (all integers). + * + * `y`, `lower` and `upper` can each be a scalar or a one-dimensional container. + * Any container arguments must be the same size. + * + * @tparam T_y type of scalar, either int or std::vector + * @tparam T_lower type of lower bound, either int or std::vector + * @tparam T_upper type of upper bound, either int or std::vector + * + * @param y integer random variable + * @param lower integer lower bound + * @param upper integer upper bound + * @return The log CCDF evaluated at the specified arguments. If containers are + * supplied, returns the sum of the log CCDFs. + * @throw std::domain_error if upper is smaller than lower. + * @throw std::invalid_argument if non-scalar arguments are of different + * sizes. + */ +template +double discrete_range_lccdf(const T_y& y, const T_lower& lower, + const T_upper& upper) { + static const char* function = "discrete_range_lccdf"; + check_consistent_sizes(function, "Lower bound parameter", lower, + "Upper bound parameter", upper); + check_greater_or_equal(function, "Upper bound parameter", upper, lower); + + if (size_zero(y, lower, upper)) { + return 0; + } + + scalar_seq_view y_vec(y); + scalar_seq_view lower_vec(lower); + scalar_seq_view upper_vec(upper); + size_t N = max_size(y, lower, upper); + + for (size_t n = 0; n < N; ++n) { + const int y_dbl = y_vec[n]; + if (y_dbl < lower_vec[n]) { + return 0; + } + if (y_dbl > upper_vec[n]) { + return LOG_ZERO; + } + } + + double ccdf(0.0); + for (size_t n = 0; n < N; n++) { + const int y_dbl = y_vec[n]; + const int lower_dbl = lower_vec[n]; + const int upper_dbl = upper_vec[n]; + ccdf += log(upper_dbl - y_dbl) - log(upper_dbl - lower_dbl + 1); + } + return ccdf; +} + +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/prim/prob/discrete_range_lcdf.hpp b/stan/math/prim/prob/discrete_range_lcdf.hpp new file mode 100644 index 00000000000..88de348c5e5 --- /dev/null +++ b/stan/math/prim/prob/discrete_range_lcdf.hpp @@ -0,0 +1,75 @@ +#ifndef STAN_MATH_PRIM_PROB_DISCRETE_RANGE_LCDF_HPP +#define STAN_MATH_PRIM_PROB_DISCRETE_RANGE_LCDF_HPP + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** \ingroup prob_dists + * Return the log CDF of a discrete range distribution for the given y, + * lower and upper bounds (all integers). + * + * `y`, `lower` and `upper` can each be a scalar or a one-dimensional container. + * Any container arguments must be the same size. + * + * @tparam T_y type of scalar, either int or std::vector + * @tparam T_lower type of lower bound, either int or std::vector + * @tparam T_upper type of upper bound, either int or std::vector + * + * @param y integer random variable + * @param lower integer lower bound + * @param upper integer upper bound + * @return The log CDF evaluated at the specified arguments. If containers are + * supplied, returns the sum of the log CDFs. + * @throw std::domain_error if upper is smaller than lower. + * @throw std::invalid_argument if non-scalar arguments are of different + * sizes. + */ +template +double discrete_range_lcdf(const T_y& y, const T_lower& lower, + const T_upper& upper) { + static const char* function = "discrete_range_lcdf"; + check_consistent_sizes(function, "Lower bound parameter", lower, + "Upper bound parameter", upper); + check_greater_or_equal(function, "Upper bound parameter", upper, lower); + + if (size_zero(y, lower, upper)) { + return 0; + } + + scalar_seq_view y_vec(y); + scalar_seq_view lower_vec(lower); + scalar_seq_view upper_vec(upper); + size_t N = max_size(y, lower, upper); + + for (size_t n = 0; n < N; ++n) { + const int y_dbl = y_vec[n]; + if (y_dbl < lower_vec[n]) { + return LOG_ZERO; + } + if (y_dbl > upper_vec[n]) { + return 0; + } + } + + double cdf(0.0); + for (size_t n = 0; n < N; n++) { + const int y_dbl = y_vec[n]; + const int lower_dbl = lower_vec[n]; + const int upper_dbl = upper_vec[n]; + cdf += log(y_dbl - lower_dbl + 1) - log(upper_dbl - lower_dbl + 1); + } + return cdf; +} + +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/prim/prob/discrete_range_log.hpp b/stan/math/prim/prob/discrete_range_log.hpp index 2c47f84eaa2..bee68587ae1 100644 --- a/stan/math/prim/prob/discrete_range_log.hpp +++ b/stan/math/prim/prob/discrete_range_log.hpp @@ -21,7 +21,7 @@ double discrete_range_log(const T_y& y, const T_lower& lower, template inline double discrete_range_log(const T_y& y, const T_lower& lower, const T_upper& upper) { - return discrete_range_lpmf(y, lower, upper); + return discrete_range_lpmf(y, lower, upper); } } // namespace math diff --git a/stan/math/prim/prob/discrete_range_lpmf.hpp b/stan/math/prim/prob/discrete_range_lpmf.hpp index 28adbf64a63..c691e5cde65 100644 --- a/stan/math/prim/prob/discrete_range_lpmf.hpp +++ b/stan/math/prim/prob/discrete_range_lpmf.hpp @@ -46,16 +46,14 @@ double discrete_range_lpmf(const T_y& y, const T_lower& lower, const T_upper& upper) { static const char* function = "discrete_range_lpmf"; using std::log; - - if (size_zero(y, lower, upper)) { - return 0.0; - } - check_not_nan(function, "Random variable", y); check_consistent_sizes(function, "Lower bound parameter", lower, "Upper bound parameter", upper); check_greater_or_equal(function, "Upper bound parameter", upper, lower); + if (size_zero(y, lower, upper)) { + return 0.0; + } if (!include_summand::value) { return 0.0; } diff --git a/stan/math/prim/prob/discrete_range_rng.hpp b/stan/math/prim/prob/discrete_range_rng.hpp index 9d2a5f56e3a..4348694ece9 100644 --- a/stan/math/prim/prob/discrete_range_rng.hpp +++ b/stan/math/prim/prob/discrete_range_rng.hpp @@ -33,11 +33,9 @@ namespace math { template inline typename VectorBuilder::type discrete_range_rng(const T_lower& lower, const T_upper& upper, RNG& rng) { + static const char* function = "discrete_range_rng"; using boost::random::uniform_int_distribution; using boost::variate_generator; - - static const char* function = "discrete_range_rng"; - check_consistent_sizes(function, "Lower bound parameter", lower, "Upper bound parameter", upper); check_greater_or_equal(function, "Upper bound parameter", upper, lower); diff --git a/test/prob/discrete_range/discrete_range_ccdf_log_test.hpp b/test/prob/discrete_range/discrete_range_ccdf_log_test.hpp new file mode 100644 index 00000000000..ebdf201079f --- /dev/null +++ b/test/prob/discrete_range/discrete_range_ccdf_log_test.hpp @@ -0,0 +1,65 @@ +// Arguments: Ints, Ints, Ints +#include + +using stan::math::var; +using std::vector; + +class AgradCcdfLogDiscreteRange : public AgradCcdfLogTest { + public: + void valid_values(vector>& parameters, + vector& ccdf_log) { + vector param(3); + + param[0] = 3; // y + param[1] = 1; // lower + param[2] = 5; // upper + parameters.push_back(param); + ccdf_log.push_back(log(2.0 / 5)); // expected ccdf_log + + param[0] = 9; // y + param[1] = 5; // lower + param[2] = 15; // upper + parameters.push_back(param); + ccdf_log.push_back(log(6.0 / 11)); // expected ccdf_log + + param[0] = 0; // y + param[1] = -4; // lower + param[2] = 5; // upper + parameters.push_back(param); + ccdf_log.push_back(log(5.0 / 10)); // expected ccdf_log + } + + void invalid_values(vector& /*index*/, vector& /*value*/) { + // y + + // lower + + // upper + } + + bool has_lower_bound() { return false; } + + bool has_upper_bound() { return false; } + + template + stan::return_type_t ccdf_log(const T_y& y, + const T_lower& lower, + const T_upper& upper, + const T3&, const T4&, + const T5&) { + return stan::math::discrete_range_lccdf(y, lower, upper); + } + + template + stan::return_type_t ccdf_log_function( + const T_y& y, const T_lower& lower, const T_upper& upper, const T3&, + const T4&, const T5&) { + if (y < lower || y > upper) { + return stan::math::LOG_ZERO; + } + + return log((upper - y) / (upper - lower + 1)); + } +}; diff --git a/test/prob/discrete_range/discrete_range_cdf_log_test.hpp b/test/prob/discrete_range/discrete_range_cdf_log_test.hpp new file mode 100644 index 00000000000..426564cf56f --- /dev/null +++ b/test/prob/discrete_range/discrete_range_cdf_log_test.hpp @@ -0,0 +1,63 @@ +// Arguments: Ints, Ints, Ints +#include + +using stan::math::var; +using std::numeric_limits; +using std::vector; + +class AgradCdfLogDiscreteRange : public AgradCdfLogTest { + public: + void valid_values(vector>& parameters, + vector& cdf_log) { + vector param(3); + + param[0] = 3; // y + param[1] = 1; // lower + param[2] = 5; // upper + parameters.push_back(param); + cdf_log.push_back(log(3.0 / 5)); // expected cdf_log + + param[0] = 9; // y + param[1] = 5; // lower + param[2] = 15; // upper + parameters.push_back(param); + cdf_log.push_back(log(5.0 / 11)); // expected cdf_log + + param[0] = 0; // y + param[1] = -4; // lower + param[2] = 5; // upper + parameters.push_back(param); + cdf_log.push_back(log(5.0 / 10)); // expected cdf_log + } + + void invalid_values(vector& /*index*/, vector& /*value*/) { + // y + + // lower + + // upper + } + + bool has_lower_bound() { return false; } + + bool has_upper_bound() { return false; } + + template + double cdf_log(const T_y& y, const T_lower& lower, const T_upper& upper, + const T3&, const T4&, const T5&) { + return stan::math::discrete_range_lcdf(y, lower, upper); + } + + template + double cdf_log_function(const T_y& y, const T_lower& lower, + const T_upper& upper, const T3&, const T4&, + const T5&) { + if (y < lower || y > upper) { + return stan::math::LOG_ZERO; + } + + return log((y - lower + 1.0) / (upper - lower + 1.0)); + } +}; diff --git a/test/prob/discrete_range/discrete_range_cdf_test.hpp b/test/prob/discrete_range/discrete_range_cdf_test.hpp new file mode 100644 index 00000000000..c0d74b41c51 --- /dev/null +++ b/test/prob/discrete_range/discrete_range_cdf_test.hpp @@ -0,0 +1,61 @@ +// Arguments: Ints, Ints, Ints +#include + +using stan::math::var; +using std::numeric_limits; +using std::vector; + +class AgradCdfDiscreteRange : public AgradCdfTest { + public: + void valid_values(vector>& parameters, vector& cdf) { + vector param(3); + + param[0] = 3; // y + param[1] = 1; // lower + param[2] = 5; // upper + parameters.push_back(param); + cdf.push_back(3.0 / 5); // expected cdf + + param[0] = 9; // y + param[1] = 5; // lower + param[2] = 15; // upper + parameters.push_back(param); + cdf.push_back(5.0 / 11); // expected cdf + + param[0] = 0; // y + param[1] = -4; // lower + param[2] = 5; // upper + parameters.push_back(param); + cdf.push_back(5.0 / 10); // expected cdf + } + + void invalid_values(vector& /*index*/, vector& /*value*/) { + // y + + // lower + + // upper + } + + bool has_lower_bound() { return false; } + + bool has_upper_bound() { return false; } + + template + double cdf(const T_y& y, const T_lower& lower, const T_upper& upper, + const T3&, const T4&, const T5&) { + return stan::math::discrete_range_cdf(y, lower, upper); + } + + template + double cdf_function(const T_y& y, const T_lower& lower, const T_upper& upper, + const T3&, const T4&, const T5&) { + if (y < lower || y > upper) { + return 0; + } + + return (y - lower + 1.0) / (upper - lower + 1.0); + } +}; diff --git a/test/unit/math/prim/prob/discrete_range_ccdf_log_test.cpp b/test/unit/math/prim/prob/discrete_range_ccdf_log_test.cpp new file mode 100644 index 00000000000..e4e8d900aba --- /dev/null +++ b/test/unit/math/prim/prob/discrete_range_ccdf_log_test.cpp @@ -0,0 +1,38 @@ +#include +#include + +TEST(ProbDiscreteRange, cdf_log_matches_lccdf) { + using stan::math::discrete_range_ccdf_log; + using stan::math::discrete_range_lccdf; + + for (int lower = 0; lower < 5; ++lower) { + for (int upper = lower; upper < 5; ++upper) { + for (int y = lower; y <= upper; ++y) { + EXPECT_FLOAT_EQ((discrete_range_lccdf(y, lower, upper)), + (discrete_range_ccdf_log(y, lower, upper))); + EXPECT_FLOAT_EQ((discrete_range_lccdf(y, lower, upper)), + (discrete_range_ccdf_log(y, lower, upper))); + } + } + } +} + +TEST(ProbDiscreteRange, lccdf_boundaries) { + using stan::math::discrete_range_lccdf; + int lower = 1; + int upper = 5; + + EXPECT_FLOAT_EQ(std::log(4.0 / 5), discrete_range_lccdf(lower, lower, upper)); + EXPECT_FLOAT_EQ(stan::math::LOG_ZERO, + discrete_range_lccdf(upper, lower, upper)); +} + +TEST(ProbDiscreteRange, lccdf_out_of_support) { + using stan::math::discrete_range_lccdf; + int lower = 1; + int upper = 5; + + EXPECT_FLOAT_EQ(0.0, discrete_range_lccdf(lower - 1, lower, upper)); + EXPECT_FLOAT_EQ(stan::math::LOG_ZERO, + discrete_range_lccdf(upper + 1, lower, upper)); +} diff --git a/test/unit/math/prim/prob/discrete_range_cdf_log_test.cpp b/test/unit/math/prim/prob/discrete_range_cdf_log_test.cpp new file mode 100644 index 00000000000..3c796f29070 --- /dev/null +++ b/test/unit/math/prim/prob/discrete_range_cdf_log_test.cpp @@ -0,0 +1,37 @@ +#include +#include + +TEST(ProbDiscreteRange, cdf_log_matches_lcdf) { + using stan::math::discrete_range_cdf_log; + using stan::math::discrete_range_lcdf; + + for (int lower = 0; lower < 5; ++lower) { + for (int upper = lower; upper < 5; ++upper) { + for (int y = lower; y <= upper; ++y) { + EXPECT_FLOAT_EQ((discrete_range_lcdf(y, lower, upper)), + (discrete_range_cdf_log(y, lower, upper))); + EXPECT_FLOAT_EQ((discrete_range_lcdf(y, lower, upper)), + (discrete_range_cdf_log(y, lower, upper))); + } + } + } +} + +TEST(ProbDiscreteRange, lcdf_boundaries) { + using stan::math::discrete_range_lcdf; + int lower = 1; + int upper = 5; + + EXPECT_FLOAT_EQ(std::log(1.0 / 5), discrete_range_lcdf(lower, lower, upper)); + EXPECT_FLOAT_EQ(0, discrete_range_lcdf(upper, lower, upper)); +} + +TEST(ProbDiscreteRange, lcdf_out_of_support) { + using stan::math::discrete_range_lcdf; + int lower = 1; + int upper = 5; + + EXPECT_FLOAT_EQ(stan::math::LOG_ZERO, + discrete_range_lcdf(lower - 1, lower, upper)); + EXPECT_FLOAT_EQ(0.0, discrete_range_lcdf(upper + 1, lower, upper)); +} diff --git a/test/unit/math/prim/prob/discrete_range_cdf_test.cpp b/test/unit/math/prim/prob/discrete_range_cdf_test.cpp new file mode 100644 index 00000000000..e5c251210d6 --- /dev/null +++ b/test/unit/math/prim/prob/discrete_range_cdf_test.cpp @@ -0,0 +1,20 @@ +#include +#include + +TEST(ProbDiscreteRange, cdf_boundaries) { + using stan::math::discrete_range_cdf; + int lower = 1; + int upper = 5; + + EXPECT_FLOAT_EQ(1.0 / 5, discrete_range_cdf(lower, lower, upper)); + EXPECT_FLOAT_EQ(5.0 / 5, discrete_range_cdf(upper, lower, upper)); +} + +TEST(ProbDiscreteRange, cdf_out_of_support) { + using stan::math::discrete_range_cdf; + int lower = 1; + int upper = 5; + + EXPECT_FLOAT_EQ(0.0, discrete_range_cdf(lower - 1, lower, upper)); + EXPECT_FLOAT_EQ(1.0, discrete_range_cdf(upper + 1, lower, upper)); +}