Skip to content

Commit fb619c0

Browse files
t4c1t4c1
committed
Merge branch 'develop' into generalize_operators
2 parents cf2b72d + dbd2df2 commit fb619c0

File tree

7 files changed

+345
-0
lines changed

7 files changed

+345
-0
lines changed

‎stan/math/prim/prob.hpp

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -75,6 +75,9 @@
7575
#include <stan/math/prim/prob/dirichlet_lpdf.hpp>
7676
#include <stan/math/prim/prob/dirichlet_lpmf.hpp>
7777
#include <stan/math/prim/prob/dirichlet_rng.hpp>
78+
#include <stan/math/prim/prob/discrete_range_log.hpp>
79+
#include <stan/math/prim/prob/discrete_range_lpmf.hpp>
80+
#include <stan/math/prim/prob/discrete_range_rng.hpp>
7881
#include <stan/math/prim/prob/double_exponential_ccdf_log.hpp>
7982
#include <stan/math/prim/prob/double_exponential_cdf.hpp>
8083
#include <stan/math/prim/prob/double_exponential_cdf_log.hpp>

‎stan/math/prim/prob/discrete_range_log.hpp

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
#ifndef STAN_MATH_PRIM_PROB_DISCRETE_RANGE_LOG_HPP
2+
#define STAN_MATH_PRIM_PROB_DISCRETE_RANGE_LOG_HPP
3+
4+
#include <stan/math/prim/prob/discrete_range_lpmf.hpp>
5+
6+
namespace stan {
7+
namespace math {
8+
9+
/** \ingroup prob_dists
10+
* @deprecated use <code>discrete_range_lpmf</code>
11+
*/
12+
template <bool propto, typename T_y, typename T_lower, typename T_upper>
13+
double discrete_range_log(const T_y& y, const T_lower& lower,
14+
const T_upper& upper) {
15+
return discrete_range_lpmf<propto, T_y, T_lower, T_upper>(y, lower, upper);
16+
}
17+
18+
/** \ingroup prob_dists
19+
* @deprecated use <code>discrete_range_lpmf</code>
20+
*/
21+
template <typename T_y, typename T_lower, typename T_upper>
22+
inline double discrete_range_log(const T_y& y, const T_lower& lower,
23+
const T_upper& upper) {
24+
return discrete_range_lpmf<T_y, T_lower, T_upper>(y, lower, upper);
25+
}
26+
27+
} // namespace math
28+
} // namespace stan
29+
#endif

‎stan/math/prim/prob/discrete_range_lpmf.hpp

Lines changed: 101 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,101 @@
1+
#ifndef STAN_MATH_PRIM_PROB_DISCRETE_RANGE_LPMF_HPP
2+
#define STAN_MATH_PRIM_PROB_DISCRETE_RANGE_LPMF_HPP
3+
4+
#include <stan/math/prim/meta.hpp>
5+
#include <stan/math/prim/err.hpp>
6+
#include <stan/math/prim/fun/constants.hpp>
7+
#include <stan/math/prim/fun/inv.hpp>
8+
#include <stan/math/prim/fun/log.hpp>
9+
#include <stan/math/prim/fun/max_size.hpp>
10+
#include <stan/math/prim/fun/size_zero.hpp>
11+
#include <stan/math/prim/fun/value_of.hpp>
12+
#include <cmath>
13+
14+
namespace stan {
15+
namespace math {
16+
17+
/** \ingroup prob_dists
18+
* Return the log PMF of a discrete range for the given y, lower and upper
19+
* bound (all integers).
20+
*
21+
\f{eqnarray*}{
22+
y &\sim& \mbox{\sf{discrete\_range}}(lower, upper) \\
23+
\log(p (y \, |\, lower, upper))
24+
&=& \log \left( \frac{1}{upper - lower + 1} \right) \\
25+
&=& -\log (upper - lower + 1)
26+
\f}
27+
*
28+
* `lower` and `upper` can each be a scalar or a one-dimensional container.
29+
* Any container arguments must be the same size.
30+
*
31+
* @tparam T_y type of scalar, either int or std::vector<int>
32+
* @tparam T_lower type of lower bound, either int or std::vector<int>
33+
* @tparam T_upper type of upper bound, either int or std::vector<int>
34+
*
35+
* @param y integer random variable
36+
* @param lower integer lower bound
37+
* @param upper integer upper bound
38+
* @return Log probability. If containers are supplied, returns the log sum
39+
* of the probabilities.
40+
* @throw std::domain_error if upper is smaller than lower.
41+
* @throw std::invalid_argument if non-scalar arguments are of different
42+
* sizes.
43+
*/
44+
template <bool propto, typename T_y, typename T_lower, typename T_upper>
45+
double discrete_range_lpmf(const T_y& y, const T_lower& lower,
46+
const T_upper& upper) {
47+
static const char* function = "discrete_range_lpmf";
48+
using std::log;
49+
50+
if (size_zero(y, lower, upper)) {
51+
return 0.0;
52+
}
53+
54+
check_not_nan(function, "Random variable", y);
55+
check_consistent_sizes(function, "Lower bound parameter", lower,
56+
"Upper bound parameter", upper);
57+
check_greater_or_equal(function, "Upper bound parameter", upper, lower);
58+
59+
if (!include_summand<propto>::value) {
60+
return 0.0;
61+
}
62+
63+
scalar_seq_view<T_y> y_vec(y);
64+
scalar_seq_view<T_lower> lower_vec(lower);
65+
scalar_seq_view<T_upper> upper_vec(upper);
66+
size_t size_lower_upper = max_size(lower, upper);
67+
size_t N = max_size(y, lower, upper);
68+
69+
for (size_t n = 0; n < N; ++n) {
70+
const double y_dbl = value_of(y_vec[n]);
71+
if (y_dbl < value_of(lower_vec[n]) || y_dbl > value_of(upper_vec[n])) {
72+
return LOG_ZERO;
73+
}
74+
}
75+
76+
VectorBuilder<true, double, T_lower, T_upper> log_upper_minus_lower(
77+
size_lower_upper);
78+
79+
for (size_t i = 0; i < size_lower_upper; i++) {
80+
const double lower_dbl = value_of(lower_vec[i]);
81+
const double upper_dbl = value_of(upper_vec[i]);
82+
log_upper_minus_lower[i] = log(upper_dbl - lower_dbl + 1);
83+
}
84+
85+
double logp(0.0);
86+
for (size_t n = 0; n < N; n++) {
87+
logp -= log_upper_minus_lower[n];
88+
}
89+
90+
return logp;
91+
}
92+
93+
template <typename T_y, typename T_lower, typename T_upper>
94+
inline double discrete_range_lpmf(const T_y& y, const T_lower& lower,
95+
const T_upper& upper) {
96+
return discrete_range_lpmf<false>(y, lower, upper);
97+
}
98+
99+
} // namespace math
100+
} // namespace stan
101+
#endif

‎stan/math/prim/prob/discrete_range_rng.hpp

Lines changed: 62 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,62 @@
1+
#ifndef STAN_MATH_PRIM_PROB_DISCRETE_RANGE_RNG_HPP
2+
#define STAN_MATH_PRIM_PROB_DISCRETE_RANGE_RNG_HPP
3+
4+
#include <stan/math/prim/meta.hpp>
5+
#include <stan/math/prim/err.hpp>
6+
#include <stan/math/prim/fun/max_size.hpp>
7+
#include <boost/random/uniform_int_distribution.hpp>
8+
#include <boost/random/variate_generator.hpp>
9+
10+
namespace stan {
11+
namespace math {
12+
13+
/** \ingroup prob_dists
14+
* Return an integer random variate between the given lower and upper bounds
15+
* (inclusive) using the specified random number generator.
16+
*
17+
* `lower` and `upper` can each be a scalar or a one-dimensional container.
18+
* Any non-scalar inputs must be the same size.
19+
*
20+
* @tparam T_lower type of lower bound, either int or std::vector<int>
21+
* @tparam T_upper type of upper bound, either int or std::vector<int>
22+
* @tparam RNG type of random number generator
23+
*
24+
* @param lower lower bound
25+
* @param upper upper bound
26+
* @param rng random number generator
27+
* @return A (sequence of) integer random variate(s) between `lower` and
28+
* `upper`, both bounds included.
29+
* @throw std::domain_error if upper is smaller than lower.
30+
* @throw std::invalid_argument if non-scalar arguments are of different
31+
* sizes.
32+
*/
33+
template <typename T_lower, typename T_upper, class RNG>
34+
inline typename VectorBuilder<true, int, T_lower, T_upper>::type
35+
discrete_range_rng(const T_lower& lower, const T_upper& upper, RNG& rng) {
36+
using boost::random::uniform_int_distribution;
37+
using boost::variate_generator;
38+
39+
static const char* function = "discrete_range_rng";
40+
41+
check_consistent_sizes(function, "Lower bound parameter", lower,
42+
"Upper bound parameter", upper);
43+
check_greater_or_equal(function, "Upper bound parameter", upper, lower);
44+
45+
scalar_seq_view<T_lower> lower_vec(lower);
46+
scalar_seq_view<T_upper> upper_vec(upper);
47+
size_t N = max_size(lower, upper);
48+
VectorBuilder<true, int, T_lower, T_upper> output(N);
49+
50+
for (size_t n = 0; n < N; ++n) {
51+
variate_generator<RNG&, uniform_int_distribution<>> discrete_range_rng(
52+
rng, uniform_int_distribution<>(lower_vec[n], upper_vec[n]));
53+
54+
output[n] = discrete_range_rng();
55+
}
56+
57+
return output.data();
58+
}
59+
60+
} // namespace math
61+
} // namespace stan
62+
#endif

‎test/prob/discrete_range/discrete_range_test.hpp

Lines changed: 62 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,62 @@
1+
// Arguments: Ints, Ints, Ints
2+
#include <stan/math/prim.hpp>
3+
4+
using std::vector;
5+
6+
class AgradDistributionsDiscreteRange : public AgradDistributionTest {
7+
public:
8+
void valid_values(vector<vector<double>>& parameters,
9+
vector<double>& log_prob) {
10+
vector<double> param(3);
11+
12+
param[0] = 3; // y
13+
param[1] = 1; // lower
14+
param[2] = 10; // upper
15+
parameters.push_back(param);
16+
log_prob.push_back(-2.302585092994045901094); // expected log_prob
17+
18+
// case for lower == upper
19+
param[0] = 5; // y
20+
param[1] = 5; // lower
21+
param[2] = 5; // upper
22+
parameters.push_back(param);
23+
log_prob.push_back(0); // expected log_prob
24+
}
25+
26+
void invalid_values(vector<size_t>& /*index*/, vector<double>& /*value*/) {
27+
// y
28+
29+
// lower
30+
31+
// upper
32+
}
33+
34+
template <class T_y, class T_lower, class T_upper, class T3, typename T4,
35+
typename T5>
36+
stan::return_type_t<T_y, T_lower, T_upper> log_prob(const T_y& y,
37+
const T_lower& lower,
38+
const T_upper& upper,
39+
const T3&, const T4&,
40+
const T5&) {
41+
return stan::math::discrete_range_lpmf(y, lower, upper);
42+
}
43+
44+
template <bool propto, class T_y, class T_lower, class T_upper, class T3,
45+
typename T4, typename T5>
46+
double log_prob(const T_y& y, const T_lower& lower, const T_upper& upper,
47+
const T3&, const T4&, const T5&) {
48+
return stan::math::discrete_range_lpmf<propto>(y, lower, upper);
49+
}
50+
51+
template <class T_y, class T_lower, class T_upper, class T3, typename T4,
52+
typename T5>
53+
double log_prob_function(const T_y& y, const T_lower& lower,
54+
const T_upper& upper, const T3&, const T4&,
55+
const T5&) {
56+
if (y < lower || y > upper) {
57+
return stan::math::LOG_ZERO;
58+
}
59+
60+
return -log(upper - lower + 1);
61+
}
62+
};

‎test/unit/math/prim/prob/discrete_range_log_test.cpp

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,24 @@
1+
#include <stan/math/prim.hpp>
2+
#include <gtest/gtest.h>
3+
4+
TEST(ProbDiscreteRange, log_matches_lpmf) {
5+
int y = 3;
6+
int lower = -2;
7+
int upper = 5;
8+
9+
EXPECT_FLOAT_EQ((stan::math::discrete_range_lpmf(y, lower, upper)),
10+
(stan::math::discrete_range_log(y, lower, upper)));
11+
EXPECT_FLOAT_EQ((stan::math::discrete_range_lpmf<true>(y, lower, upper)),
12+
(stan::math::discrete_range_log<true>(y, lower, upper)));
13+
EXPECT_FLOAT_EQ((stan::math::discrete_range_lpmf<false>(y, lower, upper)),
14+
(stan::math::discrete_range_log<false>(y, lower, upper)));
15+
EXPECT_FLOAT_EQ(
16+
(stan::math::discrete_range_lpmf<true, int, int, int>(y, lower, upper)),
17+
(stan::math::discrete_range_log<true, int, int, int>(y, lower, upper)));
18+
EXPECT_FLOAT_EQ(
19+
(stan::math::discrete_range_lpmf<false, int, int, int>(y, lower, upper)),
20+
(stan::math::discrete_range_log<false, int, int, int>(y, lower, upper)));
21+
EXPECT_FLOAT_EQ(
22+
(stan::math::discrete_range_lpmf<int, int, int>(y, lower, upper)),
23+
(stan::math::discrete_range_log<int, int, int>(y, lower, upper)));
24+
}

‎test/unit/math/prim/prob/discrete_range_test.cpp

Lines changed: 64 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,64 @@
1+
#include <stan/math/prim.hpp>
2+
#include <boost/math/distributions.hpp>
3+
#include <boost/random/mersenne_twister.hpp>
4+
#include <gtest/gtest.h>
5+
#include <vector>
6+
7+
TEST(ProbDistributionsDiscreteRange, error_check) {
8+
using stan::math::discrete_range_rng;
9+
boost::random::mt19937 rng;
10+
11+
std::vector<int> lower{5, 11, -15};
12+
std::vector<int> upper{7, 15, -10};
13+
14+
EXPECT_NO_THROW(discrete_range_rng(lower, upper, rng));
15+
EXPECT_THROW(discrete_range_rng(lower, 10, rng), std::domain_error);
16+
EXPECT_THROW(discrete_range_rng(10, upper, rng), std::domain_error);
17+
18+
std::vector<int> vec2{1, 2};
19+
EXPECT_THROW(discrete_range_rng(lower, vec2, rng), std::invalid_argument);
20+
EXPECT_THROW(discrete_range_rng(vec2, upper, rng), std::invalid_argument);
21+
}
22+
23+
TEST(ProbDistributionsDiscreteRange, boundary_values) {
24+
using stan::math::discrete_range_rng;
25+
boost::random::mt19937 rng;
26+
27+
std::vector<int> lower{-5, 11, 17};
28+
EXPECT_EQ(lower, discrete_range_rng(lower, lower, rng));
29+
30+
std::vector<int> upper(lower);
31+
for (int i = 0; i < upper.size(); i++) {
32+
++upper[i];
33+
}
34+
35+
EXPECT_LE(lower, discrete_range_rng(lower, upper, rng));
36+
EXPECT_GE(upper, discrete_range_rng(lower, upper, rng));
37+
}
38+
39+
TEST(ProbDistributionsDiscreteRange, chiSquareGoodnessFitTest) {
40+
boost::random::mt19937 rng;
41+
42+
int N = 10000;
43+
int lower = -3;
44+
int upper = 20;
45+
46+
int K = upper - lower + 1;
47+
boost::math::chi_squared mydist(K - 1);
48+
49+
std::vector<int> bin(K, 0);
50+
std::vector<double> expect(K, static_cast<double>(N) / K);
51+
52+
for (int count = 0; count < N; ++count) {
53+
int a = stan::math::discrete_range_rng(lower, upper, rng);
54+
++bin[a - lower];
55+
}
56+
57+
double chi = 0;
58+
59+
for (int j = 0; j < K; j++) {
60+
chi += (bin[j] - expect[j]) * (bin[j] - expect[j]) / expect[j];
61+
}
62+
63+
EXPECT_TRUE(chi < quantile(complement(mydist, 1e-6)));
64+
}

0 commit comments

Comments
 (0)