Skip to content
5 changes: 5 additions & 0 deletions stan/math/prim/prob.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -75,6 +75,11 @@
#include <stan/math/prim/prob/dirichlet_lpdf.hpp>
#include <stan/math/prim/prob/dirichlet_lpmf.hpp>
#include <stan/math/prim/prob/dirichlet_rng.hpp>
#include <stan/math/prim/prob/discrete_range_ccdf_log.hpp>
#include <stan/math/prim/prob/discrete_range_cdf.hpp>
#include <stan/math/prim/prob/discrete_range_cdf_log.hpp>
#include <stan/math/prim/prob/discrete_range_lccdf.hpp>
#include <stan/math/prim/prob/discrete_range_lcdf.hpp>
#include <stan/math/prim/prob/discrete_range_log.hpp>
#include <stan/math/prim/prob/discrete_range_lpmf.hpp>
#include <stan/math/prim/prob/discrete_range_rng.hpp>
Expand Down
20 changes: 20 additions & 0 deletions stan/math/prim/prob/discrete_range_ccdf_log.hpp
Original file line number Diff line number Diff line change
@@ -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 <stan/math/prim/prob/discrete_range_lccdf.hpp>

namespace stan {
namespace math {

/** \ingroup prob_dists
* @deprecated use <code>discrete_range_lccdf</code>
*/
template <typename T_y, typename T_lower, typename T_upper>
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
75 changes: 75 additions & 0 deletions stan/math/prim/prob/discrete_range_cdf.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,75 @@
#ifndef STAN_MATH_PRIM_PROB_DISCRETE_RANGE_CDF_HPP
#define STAN_MATH_PRIM_PROB_DISCRETE_RANGE_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/inv.hpp>
#include <stan/math/prim/fun/log.hpp>
#include <stan/math/prim/fun/max_size.hpp>
#include <stan/math/prim/fun/size_zero.hpp>
#include <cmath>

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<int>
* @tparam T_lower type of lower bound, either int or std::vector<int>
* @tparam T_upper type of upper bound, either int or std::vector<int>
*
* @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 <typename T_y, typename T_lower, typename T_upper>
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<T_y> y_vec(y);
scalar_seq_view<T_lower> lower_vec(lower);
scalar_seq_view<T_upper> 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];

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[optional]
We have not in general marked our local args as const. It's not a mistake per se, but it clutters the code up a lot.

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++) {

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[optional]
This can be done more efficiently. The point is to compute things like upper_dbl - lower_dbl ahead of time in their own vector rather than in the N loop because they might both be scalars (so it's only done once) whereas y is a vector, so N > 1. Overall, though, this isn't an important function for performance as there's no autodiff and it's unlikely to see use in many models. Hence the optionality.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I also thought that the operations that could be precomputed here were only sums, rather than more involved operations, so I'll leave this. I have a PR for other distributions that tries to clean up this kind of operations.

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
20 changes: 20 additions & 0 deletions stan/math/prim/prob/discrete_range_cdf_log.hpp
Original file line number Diff line number Diff line change
@@ -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 <stan/math/prim/prob/discrete_range_lcdf.hpp>

namespace stan {
namespace math {

/** \ingroup prob_dists
* @deprecated use <code>discrete_range_lcdf</code>
*/
template <typename T_y, typename T_lower, typename T_upper>
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
75 changes: 75 additions & 0 deletions stan/math/prim/prob/discrete_range_lccdf.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,75 @@
#ifndef STAN_MATH_PRIM_PROB_DISCRETE_RANGE_LCCDF_HPP
#define STAN_MATH_PRIM_PROB_DISCRETE_RANGE_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/inv.hpp>
#include <stan/math/prim/fun/log.hpp>
#include <stan/math/prim/fun/max_size.hpp>
#include <stan/math/prim/fun/size_zero.hpp>
#include <cmath>

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<int>
* @tparam T_lower type of lower bound, either int or std::vector<int>
* @tparam T_upper type of upper bound, either int or std::vector<int>
*
* @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 <typename T_y, typename T_lower, typename T_upper>
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<T_y> y_vec(y);
scalar_seq_view<T_lower> lower_vec(lower);
scalar_seq_view<T_upper> 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
75 changes: 75 additions & 0 deletions stan/math/prim/prob/discrete_range_lcdf.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,75 @@
#ifndef STAN_MATH_PRIM_PROB_DISCRETE_RANGE_LCDF_HPP
#define STAN_MATH_PRIM_PROB_DISCRETE_RANGE_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/inv.hpp>
#include <stan/math/prim/fun/log.hpp>
#include <stan/math/prim/fun/max_size.hpp>
#include <stan/math/prim/fun/size_zero.hpp>
#include <cmath>

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<int>
* @tparam T_lower type of lower bound, either int or std::vector<int>
* @tparam T_upper type of upper bound, either int or std::vector<int>
*
* @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 <typename T_y, typename T_lower, typename T_upper>
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<T_y> y_vec(y);
scalar_seq_view<T_lower> lower_vec(lower);
scalar_seq_view<T_upper> 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);
Comment thread
mcol marked this conversation as resolved.
}
return cdf;
}

} // namespace math
} // namespace stan
#endif
2 changes: 1 addition & 1 deletion stan/math/prim/prob/discrete_range_log.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ double discrete_range_log(const T_y& y, const T_lower& lower,
template <typename T_y, typename T_lower, typename T_upper>
inline double discrete_range_log(const T_y& y, const T_lower& lower,
const T_upper& upper) {
return discrete_range_lpmf<T_y, T_lower, T_upper>(y, lower, upper);
return discrete_range_lpmf(y, lower, upper);
}

} // namespace math
Expand Down
8 changes: 3 additions & 5 deletions stan/math/prim/prob/discrete_range_lpmf.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<propto>::value) {
return 0.0;
}
Expand Down
4 changes: 1 addition & 3 deletions stan/math/prim/prob/discrete_range_rng.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -33,11 +33,9 @@ namespace math {
template <typename T_lower, typename T_upper, class RNG>
inline typename VectorBuilder<true, int, T_lower, T_upper>::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);
Expand Down
65 changes: 65 additions & 0 deletions test/prob/discrete_range/discrete_range_ccdf_log_test.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,65 @@
// Arguments: Ints, Ints, Ints
#include <stan/math/prim.hpp>

using stan::math::var;
using std::vector;

class AgradCcdfLogDiscreteRange : public AgradCcdfLogTest {
public:
void valid_values(vector<vector<double>>& parameters,
vector<double>& ccdf_log) {
vector<double> 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<size_t>& /*index*/, vector<double>& /*value*/) {
// y

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These are going to need the out-of-support tests.

In general, this shouldn't have three lines of comments double spaced. A simple // all values valid comment is better because it doesn't look like someone forgot to fill in values here.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've added prim test for out-of-support values. The invalid_values section of distribution test doesn't let you express that easily, they are designed to capture values that are always invalid (say a negative rate parameter in a Poisson distribution). The formatting follows what's done in all other tests.


// lower

// upper
}

bool has_lower_bound() { return false; }

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not exactly sure how these tests work, but this distribution has a lower and upper bound.

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Those are used for distributions that have a fixed support, say the beta that is in [0,1]. There is a kludge for uniform that sets only the upper bound. So I'm not sure it's worth trying to use those here, I can add some prim tests about boundary behaviour.


bool has_upper_bound() { return false; }

template <class T_y, class T_lower, class T_upper, typename T3, typename T4,
typename T5>
stan::return_type_t<T_y, T_lower, T_upper> 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 <class T_y, class T_lower, class T_upper, typename T3, typename T4,
typename T5>
stan::return_type_t<T_y, T_lower, T_upper> 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;
Comment thread
mcol marked this conversation as resolved.
}

return log((upper - y) / (upper - lower + 1));
}
};
Loading