Skip to content

Commit 575cff4

Browse files
rok-cesnovarrok-cesnovar
authored
Merge pull request #1938 from dirmeier/feature/issue-200-add_poisson_binomial_pmf
Feature/issue 869: add poisson binomial pmf/rng/cdf/ccdf
2 parents 3ac7b31 + 348e8c5 commit 575cff4

13 files changed

Lines changed: 839 additions & 1 deletion

‎.gitignore‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -92,3 +92,4 @@ lib/tbb/**
9292

9393
# gdb debug history
9494
.gdb_history
95+

‎stan/math/prim/fun.hpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -229,6 +229,7 @@
229229
#include <stan/math/prim/fun/owens_t.hpp>
230230
#include <stan/math/prim/fun/Phi.hpp>
231231
#include <stan/math/prim/fun/Phi_approx.hpp>
232+
#include <stan/math/prim/fun/poisson_binomial_log_probs.hpp>
232233
#include <stan/math/prim/fun/polar.hpp>
233234
#include <stan/math/prim/fun/positive_constrain.hpp>
234235
#include <stan/math/prim/fun/positive_free.hpp>

‎stan/math/prim/fun/poisson_binomial_log_probs.hpp‎

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
#ifndef STAN_MATH_PRIM_FUN_POISSON_BINOMIAL_LOG_PROBS_HPP
2+
#define STAN_MATH_PRIM_FUN_POISSON_BINOMIAL_LOG_PROBS_HPP
3+
4+
#include <stan/math/prim/meta.hpp>
5+
#include <stan/math/prim/fun/log.hpp>
6+
#include <stan/math/prim/fun/log1m.hpp>
7+
#include <stan/math/prim/fun/log_sum_exp.hpp>
8+
#include <stan/math/prim/fun/max_size.hpp>
9+
10+
namespace stan {
11+
namespace math {
12+
13+
/**
14+
* Returns the last row of the log probability matrix of the Poisson-Binomial
15+
* distribution given the number of successes and a vector of success
16+
* probabilities.
17+
*
18+
* @tparam T_theta template expression
19+
* @param y numbers of successes
20+
* @param theta N-dimensional vector of success probabilities for each trial
21+
* @return the last row of the computed log probability matrix
22+
*/
23+
template <typename T_theta, typename T_scalar = scalar_type_t<T_theta>,
24+
require_eigen_col_vector_t<T_theta>* = nullptr>
25+
plain_type_t<T_theta> poisson_binomial_log_probs(int y, const T_theta& theta) {
26+
int size_theta = theta.size();
27+
plain_type_t<T_theta> log_theta = log(theta);
28+
plain_type_t<T_theta> log1m_theta = log1m(theta);
29+
30+
Eigen::Matrix<T_scalar, Eigen::Dynamic, Eigen::Dynamic> alpha(size_theta + 1,
31+
y + 1);
32+
33+
// alpha[i, j] = log prob of j successes in first i trials
34+
alpha(0, 0) = 0.0;
35+
for (int i = 0; i < size_theta; ++i) {
36+
// no success in i trials
37+
alpha(i + 1, 0) = alpha(i, 0) + log1m_theta[i];
38+
39+
// 0 < j < i successes in i trials
40+
for (int j = 0; j < std::min(y, i); ++j) {
41+
alpha(i + 1, j + 1) = log_sum_exp(alpha(i, j) + log_theta[i],
42+
alpha(i, j + 1) + log1m_theta[i]);
43+
}
44+
45+
// i successes in i trials
46+
if (i < y) {
47+
alpha(i + 1, i + 1) = alpha(i, i) + log_theta(i);
48+
}
49+
}
50+
51+
return alpha.row(size_theta);
52+
}
53+
54+
template <typename T_theta, typename T_scalar = scalar_type_t<T_theta>>
55+
auto poisson_binomial_log_probs(const std::vector<int>& y,
56+
const T_theta& theta) {
57+
size_t max_sizes = std::max(stan::math::size(y), size_mvt(theta));
58+
std::vector<Eigen::Matrix<T_scalar, Eigen::Dynamic, 1>> result(max_sizes);
59+
vector_seq_view<T_theta> theta_vec(theta);
60+
61+
for (size_t i = 0; i < max_sizes; ++i) {
62+
result[i] = poisson_binomial_log_probs(y[i], theta_vec[i]);
63+
}
64+
65+
return result;
66+
}
67+
68+
} // namespace math
69+
} // namespace stan
70+
#endif

‎stan/math/prim/prob.hpp‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -261,6 +261,10 @@
261261
#include <stan/math/prim/prob/pareto_type_2_log.hpp>
262262
#include <stan/math/prim/prob/pareto_type_2_lpdf.hpp>
263263
#include <stan/math/prim/prob/pareto_type_2_rng.hpp>
264+
#include <stan/math/prim/prob/poisson_binomial_lccdf.hpp>
265+
#include <stan/math/prim/prob/poisson_binomial_lcdf.hpp>
266+
#include <stan/math/prim/prob/poisson_binomial_lpmf.hpp>
267+
#include <stan/math/prim/prob/poisson_binomial_rng.hpp>
264268
#include <stan/math/prim/prob/poisson_ccdf_log.hpp>
265269
#include <stan/math/prim/prob/poisson_cdf.hpp>
266270
#include <stan/math/prim/prob/poisson_cdf_log.hpp>

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

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@ namespace math {
2222
*
2323
* @tparam T_n type of successes parameter
2424
* @tparam T_N type of population size parameter
25-
* @tparam theta type of chance of success parameter
25+
* @tparam T_prob type of chance of success parameter
2626
* @param n successes parameter
2727
* @param N population size parameter
2828
* @param theta chance of success parameter

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

Lines changed: 67 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,67 @@
1+
#ifndef STAN_MATH_PRIM_PROB_POISSON_BINOMIAL_LCCDF_HPP
2+
#define STAN_MATH_PRIM_PROB_POISSON_BINOMIAL_LCCDF_HPP
3+
4+
#include <stan/math/prim/meta.hpp>
5+
#include <stan/math/prim/err.hpp>
6+
#include <stan/math/prim/fun/binomial_coefficient_log.hpp>
7+
#include <stan/math/prim/fun/log.hpp>
8+
#include <stan/math/prim/fun/log1m.hpp>
9+
#include <stan/math/prim/fun/log_sum_exp.hpp>
10+
#include <stan/math/prim/fun/log1m_exp.hpp>
11+
#include <stan/math/prim/fun/max_size.hpp>
12+
#include <stan/math/prim/fun/multiply_log.hpp>
13+
#include <stan/math/prim/fun/size.hpp>
14+
#include <stan/math/prim/fun/size_zero.hpp>
15+
#include <stan/math/prim/fun/value_of.hpp>
16+
#include <stan/math/prim/fun/poisson_binomial_log_probs.hpp>
17+
18+
namespace stan {
19+
namespace math {
20+
21+
/** \ingroup prob_dists
22+
* Returns the log CCDF for the Poisson-binomial distribution evaluated at the
23+
* specified number of successes and probabilities of successes.
24+
*
25+
* @tparam T_y type of number of successes parameter
26+
* @tparam T_theta type of chance of success parameters
27+
* @param y input scalar, vector, or nested vector of numbers of successes
28+
* @param theta array of chances of success parameters
29+
* @return sum of log probabilities
30+
* @throw std::domain_error if y is out of bounds
31+
* @throw std::domain_error if theta is not a valid vector of probabilities
32+
* @throw std::invalid_argument If y and theta are different lengths
33+
*/
34+
template <bool propto, typename T_y, typename T_theta>
35+
return_type_t<T_theta> poisson_binomial_lccdf(const T_y& y,
36+
const T_theta& theta) {
37+
static const char* function = "poisson_binomial_lccdf";
38+
39+
size_t size_theta = size_mvt(theta);
40+
if (size_theta > 1) {
41+
check_consistent_sizes(function, "Successes variables", y,
42+
"Probability parameters", theta);
43+
}
44+
45+
size_t max_sz = std::max(stan::math::size(y), size_mvt(theta));
46+
scalar_seq_view<T_y> y_vec(y);
47+
vector_seq_view<T_theta> theta_vec(theta);
48+
49+
for (size_t i = 0; i < max_sz; ++i) {
50+
check_bounded(function, "Successes variable", y_vec[i], 0,
51+
theta_vec[i].size());
52+
check_finite(function, "Probability parameters", theta_vec[i]);
53+
check_bounded(function, "Probability parameters", theta_vec[i], 0.0, 1.0);
54+
}
55+
56+
return sum(log1m_exp(log_sum_exp(poisson_binomial_log_probs(y, theta))));
57+
}
58+
59+
template <typename T_y, typename T_theta>
60+
return_type_t<T_theta> poisson_binomial_lccdf(const T_y& y,
61+
const T_theta& theta) {
62+
return poisson_binomial_lccdf<false>(y, theta);
63+
}
64+
65+
} // namespace math
66+
} // namespace stan
67+
#endif

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

Lines changed: 66 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,66 @@
1+
#ifndef STAN_MATH_PRIM_PROB_POISSON_BINOMIAL_LCDF_HPP
2+
#define STAN_MATH_PRIM_PROB_POISSON_BINOMIAL_LCDF_HPP
3+
4+
#include <stan/math/prim/meta.hpp>
5+
#include <stan/math/prim/err.hpp>
6+
#include <stan/math/prim/fun/binomial_coefficient_log.hpp>
7+
#include <stan/math/prim/fun/log.hpp>
8+
#include <stan/math/prim/fun/log1m.hpp>
9+
#include <stan/math/prim/fun/log_sum_exp.hpp>
10+
#include <stan/math/prim/fun/max_size.hpp>
11+
#include <stan/math/prim/fun/multiply_log.hpp>
12+
#include <stan/math/prim/fun/size.hpp>
13+
#include <stan/math/prim/fun/size_zero.hpp>
14+
#include <stan/math/prim/fun/value_of.hpp>
15+
#include <stan/math/prim/fun/poisson_binomial_log_probs.hpp>
16+
17+
namespace stan {
18+
namespace math {
19+
20+
/** \ingroup prob_dists
21+
* Returns the log CDF for the Poisson-binomial distribution evaluated at the
22+
* specified number of successes and probabilities of successes.
23+
*
24+
* @tparam T_y type of number of successes parameter
25+
* @tparam T_theta type of chance of success parameters
26+
* @param y input scalar, vector, or nested vector of numbers of successes
27+
* @param theta array of chances of success parameters
28+
* @return sum of log probabilities
29+
* @throw std::domain_error if y is out of bounds
30+
* @throw std::domain_error if theta is not a valid vector of probabilities
31+
* @throw std::invalid_argument If y and theta are different lengths
32+
*/
33+
template <bool propto, typename T_y, typename T_theta>
34+
return_type_t<T_theta> poisson_binomial_lcdf(const T_y& y,
35+
const T_theta& theta) {
36+
static const char* function = "poisson_binomial_lcdf";
37+
38+
size_t size_theta = size_mvt(theta);
39+
if (size_theta > 1) {
40+
check_consistent_sizes(function, "Successes variables", y,
41+
"Probability parameters", theta);
42+
}
43+
44+
size_t max_sz = std::max(stan::math::size(y), size_theta);
45+
scalar_seq_view<T_y> y_vec(y);
46+
vector_seq_view<T_theta> theta_vec(theta);
47+
48+
for (size_t i = 0; i < max_sz; ++i) {
49+
check_bounded(function, "Successes variable", y_vec[i], 0,
50+
theta_vec[i].size());
51+
check_finite(function, "Probability parameters", theta_vec[i]);
52+
check_bounded(function, "Probability parameters", theta_vec[i], 0.0, 1.0);
53+
}
54+
55+
return sum(log_sum_exp(poisson_binomial_log_probs(y, theta)));
56+
}
57+
58+
template <typename T_y, typename T_theta>
59+
return_type_t<T_theta> poisson_binomial_lcdf(const T_y& y,
60+
const T_theta& theta) {
61+
return poisson_binomial_lcdf<false>(y, theta);
62+
}
63+
64+
} // namespace math
65+
} // namespace stan
66+
#endif

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

Lines changed: 64 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,64 @@
1+
#ifndef STAN_MATH_PRIM_PROB_POISSON_BINOMIAL_LPMF_HPP
2+
#define STAN_MATH_PRIM_PROB_POISSON_BINOMIAL_LPMF_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 <stan/math/prim/fun/poisson_binomial_log_probs.hpp>
8+
9+
namespace stan {
10+
namespace math {
11+
12+
/** \ingroup prob_dists
13+
* Returns the log PMF for the Poisson-binomial distribution evaluated at an
14+
* specified array of numbers of successes and probabilities of successes.
15+
*
16+
* @tparam T_y type of number of successes parameter
17+
* @tparam T_theta type of chance of success parameters
18+
* @param y input scalar, vector, or nested vector of numbers of successes
19+
* @param theta array of chances of success parameters
20+
* @return sum of log probabilities
21+
* @throw std::domain_error if y is out of bounds
22+
* @throw std::domain_error if theta is not a valid vector of probabilities
23+
* @throw std::invalid_argument If y and theta are different lengths
24+
*/
25+
template <bool propto, typename T_y, typename T_theta>
26+
return_type_t<T_theta> poisson_binomial_lpmf(const T_y& y,
27+
const T_theta& theta) {
28+
static const char* function = "poisson_binomial_lpmf";
29+
30+
size_t size_theta = size_mvt(theta);
31+
if (size_theta > 1) {
32+
check_consistent_sizes(function, "Successes variables", y,
33+
"Probability parameters", theta);
34+
}
35+
36+
size_t max_sz = std::max(stan::math::size(y), size_theta);
37+
scalar_seq_view<T_y> y_vec(y);
38+
vector_seq_view<T_theta> theta_vec(theta);
39+
40+
for (size_t i = 0; i < max_sz; ++i) {
41+
check_bounded(function, "Successes variable", y_vec[i], 0,
42+
theta_vec[i].size());
43+
check_finite(function, "Probability parameters", theta_vec[i]);
44+
check_bounded(function, "Probability parameters", theta_vec[i], 0.0, 1.0);
45+
}
46+
47+
return_type_t<T_theta> log_prob = 0.0;
48+
for (size_t i = 0; i < max_sz; ++i) {
49+
auto x = poisson_binomial_log_probs(y_vec[i], theta_vec[i]);
50+
log_prob += x(y_vec[i]);
51+
}
52+
53+
return log_prob;
54+
}
55+
56+
template <typename T_y, typename T_theta>
57+
return_type_t<T_theta> poisson_binomial_lpmf(const T_y& y,
58+
const T_theta& theta) {
59+
return poisson_binomial_lpmf<false>(y, theta);
60+
}
61+
62+
} // namespace math
63+
} // namespace stan
64+
#endif

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

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
#ifndef STAN_MATH_PRIM_PROB_POISSON_BINOMIAL_RNG_HPP
2+
#define STAN_MATH_PRIM_PROB_POISSON_BINOMIAL_RNG_HPP
3+
4+
#include <stan/math/prim/meta.hpp>
5+
#include <stan/math/prim/err.hpp>
6+
#include <stan/math/prim/fun/size.hpp>
7+
#include <boost/random/bernoulli_distribution.hpp>
8+
#include <boost/random/variate_generator.hpp>
9+
10+
namespace stan {
11+
namespace math {
12+
13+
/** \ingroup prob_dists
14+
* Return a pseudorandom Poisson binomial random variable for the given vector
15+
* of success parameters using the specified random number
16+
* generator.
17+
*
18+
* @tparam T_theta Type of change of success parameter
19+
* @tparam RNG class of rng
20+
* @param theta (Sequence of) chance of success parameter(s)
21+
* @param rng random number generator
22+
* @return a Poisson binomial distribution random variable
23+
* @throw std::domain_error if theta is not a valid probability
24+
*/
25+
template <typename T_theta, typename RNG>
26+
inline int poisson_binomial_rng(
27+
const Eigen::Matrix<T_theta, Eigen::Dynamic, 1>& theta, RNG& rng) {
28+
static const char* function = "poisson_binomial_rng";
29+
check_finite(function, "Probability parameters", theta);
30+
check_bounded(function, "Probability parameters", theta, 0.0, 1.0);
31+
32+
int y = 0;
33+
for (size_t i = 0; i < theta.size(); ++i) {
34+
boost::variate_generator<RNG&, boost::bernoulli_distribution<> >
35+
bernoulli_rng(rng, boost::bernoulli_distribution<>(theta(i)));
36+
y += bernoulli_rng();
37+
}
38+
39+
return y;
40+
}
41+
42+
} // namespace math
43+
} // namespace stan
44+
#endif

0 commit comments

Comments
 (0)