negative_binomial_distribution operator() doesn't properly handle special cases k == 0, p == 0, p == 1
- Dominant language
- C++
- Stars
- 42
- Forks
- 76
- Avg merge
- 4d 8h
- Merged PRs (30d)
- 5
Description
Currently negative_binomial_distribution constructor [requires](https://github.com/boostorg/random/blob/8c30753bf49504b5f677b8b410a3787907c86bd0/include/boost/random/negative_binomial_distribution.hpp#L97) k >=0 && 0 <= p <= 1, but [operator()](https://github.com/boostorg/random/blob/8c30753bf49504b5f677b8b410a3787907c86bd0/include/boost/random/negative_binomial_distribution.hpp#L117) is implemented as follows:
```
template
IntType operator()(URNG& urng) const
{
gamma_distribution gamma(_k, (1-_p)/_p);
poisson_distribution poisson(gamma(urng));
return poisson(urng);
}
```
Since gamma_distribution constructor [requires](https://github.com/boostorg/random/blob/21d66a0a2cd2864b2c592b6118370b3203d31e31/include/boost/random/gamma_distribution.hpp#L111) both input parameters to be > 0, the code violates preconditions when having _k == 0 or _p == 1; also case _p == 0 leads to invalid results.
- Case p == 1: Would seem easy to fix by adding `if (_p == 1) return 0;` given the trivial distribution.
- Case p == 0: is almost redundant and for example C++ standard [1] and many other definitions [2] [3] require p > 0: p can be zero only if k == 0, otherwise the distribution is not well defined (P(i|k,0) is zero for all i >= 0 when k > 0).
- Case k == 0: for example C++ standard [1] and many other definitions [2] [3] require k > 0.
**Example program**
```
#include
#include
int main()
{
std::mt19937 randEng;
boost::random::negative_binomial_distribution(0, 0.5)(randEng); // BOOST_ASSERT() fails in gamma_distribution constructor
boost::random::negative_binomial_distribution(10, 0)(randEng); // Returns bogus (gamma_distribution gets inf as second parameter)
boost::random::negative_binomial_distribution(10, 1)(randEng); // BOOST_ASSERT() fails in gamma_distribution constructor
boost::random::negative_binomial_distribution(0, 0)(randEng); // BOOST_ASSERT() fails in gamma_distribution constructor
return 0;
}
```
[1]: Checked from draft N4842 (2019-11-27)
[2] https://se.mathworks.com/help/stats/prob.negativebinomialdistribution.html
[3] http://search.r-project.org/R/library/stats/html/NegBinomial.html
Contributor guide
No contributing guide indexed for this repository
Assessment
This issue has not been assessed yet.