This is a Beta negative binomial distribution, with parameter $r=1$ in your case, using the Wikipedia notation. It also named Beta-Pascal distribution when $r$ is an integer. As you noted in a comment, this is a predictive distribution in the Bayesian negative binomial model with a conjugate Beta prior on the success probability.
Thus you can sample it by sampling a $\text{Beta}(\alpha,\beta)$ variable $u$ and then sampling a negative binomial variable $\text{NB}(r,u)$ (with $r=1$ in your case, that is to say a geometric distribution).
This distribution is implemented in the R package brr
. The sampler has name rbeta_nbinom
, the pmf has name dbeta_nbinom
, etc. The notations are $a=r$, $c=\alpha$, $d=\beta$. Check:
> Alpha <- 2; Beta <- 3
> a <- 1
> all.equal(brr::dbeta_nbinom(0:10, a, Alpha, Beta), beta(Alpha+a, Beta+0:10)/beta(Alpha,Beta))
[1] TRUE
Looking at the code, one can see it actually calls the ghyper
(generalized hypergeometric) family of distributions of the SuppDists
package:
brr::rbeta_nbinom
function(n, a, c, d){
rghyper(n, -d, -a, c-1)
}
Ineed, the BNB distribution is known as a type IV generalized hypergeometric distribution. See the help of ghyper
in the SuppDists
package. I believe this can also be found in Johnson & al's book Univariate Discrete Distributions.