Feature/issue 3121 beta neg binomial rng - #3126
Conversation
Sorry, something went wrong.
Jenkins Console Log Machine information No LSB modules are available. Distributor ID: Ubuntu Description: Ubuntu 20.04.3 LTS Release: 20.04 Codename: focal CPU: G++: Clang: |
Sorry, something went wrong.
SteveBronder
left a comment
There was a problem hiding this comment.
Two small changes then this is good to go!
Sorry, something went wrong.
| for (size_t n = 0; n < size_p; ++n) { | ||
| odds_ratio_p[n] = stan::math::exp(stan::math::log(p_vec.val(n)) | ||
| - stan::math::log((1 - p_vec.val(n)))); |
There was a problem hiding this comment.
You can add #include <stan/math/prim/fun/log.hpp> at the top to remove stan::math::log here
Sorry, something went wrong.
There was a problem hiding this comment.
This should be log1m as well
Sorry, something went wrong.
There was a problem hiding this comment.
Thanks! Will change it
Sorry, something went wrong.
| template <typename T_r, typename T_alpha, typename T_beta, class RNG> | ||
| inline typename VectorBuilder<true, int, T_r, T_alpha, T_beta>::type | ||
| beta_neg_binomial_rng(const T_r &r, const T_alpha &alpha, const T_beta &beta, | ||
| RNG &rng) { |
There was a problem hiding this comment.
| template <typename T_r, typename T_alpha, typename T_beta, class RNG> | |
| inline typename VectorBuilder<true, int, T_r, T_alpha, T_beta>::type | |
| beta_neg_binomial_rng(const T_r &r, const T_alpha &alpha, const T_beta &beta, | |
| RNG &rng) { | |
| template <typename T_r, typename T_alpha, typename T_beta, typename RNG> | |
| inline auto beta_neg_binomial_rng(const T_r& r, const T_alpha& alpha, const T_beta& beta, RNG& rng) { |
Sorry, something went wrong.
Jenkins Console Log Machine information No LSB modules are available. Distributor ID: Ubuntu Description: Ubuntu 20.04.3 LTS Release: 20.04 Codename: focal CPU: G++: Clang: |
Sorry, something went wrong.
SteveBronder
left a comment
There was a problem hiding this comment.
Looks good!
Sorry, something went wrong.
Summary
With this PR the rng of beta negative binomial distribution are added.
See issue #3121
Expressions involved :
for $$Y\sim \text{BNB}(r,\alpha,\beta)$$ we have
$$ \begin{align} Y &\sim \text{NB}(r,p) \\ p &\sim {\textrm {Beta}}(\alpha ,\beta ), \end{align} $$
where $\text{NB}(y~|~r,p) = \binom {y+r-1}{y} (1-p)^{y} p^{r}.$
In Stan, the negative binomial distribution is given by
$$ \text{NegBinomial}(y~|~\alpha,\beta) = \binom{y + \alpha - 1}{y} \left( \frac{\beta}{\beta+1} \right)^{\alpha} \left( \frac{1}{\beta + 1} \right)^{y} $$
so
$$p = \frac{\beta}{\beta+1} \Rightarrow \beta = \frac{p}{1-p}$$
Therefore, provided $r, \alpha, \beta$, we generate $$p \sim \text{Beta}(\alpha ,\beta )$$ and then sample $Y \sim \text{NB}(r, \frac{p}{1-p})$.
Tests
Test is written follow test/unit/math/prim/prob/beta_binomial_test.cpp
Side Effects
Under more extreme paremeter combination that makes the tail of the distribution heavy, neg_binomial_rng might throw domain_error because rng_from_gamma >= POISSON_MAX_RATE. See
math/stan/math/prim/prob/neg_binomial_rng.hpp
Line 60 in 4a812be
This is also a problem for neg_binomial_rng and poisson_rng. But due to their nature, this situation is relatively rare.
The reason behind this is that the Poisson random number generator in Boost will overflow when it exceeds the type range. Stan uses 32-bit int, so there is an upper limit. The current code includes POISSON_MAX_RATE to prevent overflow.
Release notes
beta_neg_binomial_rng is available if merged.
Checklist
Copyright holder: (fill in copyright holder information)
The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses:
- Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
- Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)
the basic tests are passing
the code is written in idiomatic C++ and changes are documented in the doxygen
the new changes are tested