鉴于随着减少，我建议生成一个统一变量并计算累积和直到实现然后等于对应的。自从

$$\frac{\text{Beta}(\alpha+1, \beta+x)}{\text{Beta}(\alpha,\beta)}=\dfrac{\alpha}{\alpha+\beta+x}\dfrac{\beta+x-1}{\alpha+\beta+x-1}\cdots\dfrac{\beta}{\alpha+\beta}$$ $x$$u\sim\mathcal{U}(0,1)$ $$S_k=\sum_{x=0}^k \frac{\text{Beta}(\alpha+1, \beta+x)}{\text{Beta}(\alpha,\beta)}$$ $$S_k>u$$ $k$ $$\eqalign{R_x&=\frac{\text{Beta}(\alpha+1, \beta+x)}{\text{Beta}(\alpha,\beta)}\\&=\dfrac{\alpha}{\alpha+\beta+x}\dfrac{\beta+x-1}{\alpha+\beta+x-1}\cdots\dfrac{\beta}{\alpha+\beta}\\&=\frac{\alpha+\beta+x-1}{\alpha+\beta+x}\frac{\beta+x-1}{\alpha+\beta+x-1}R_{x-1}\\&=\frac{\beta+x-1}{\alpha+\beta+x}R_{x-1}}$$ 和计算可以完全避免使用 Gamma 函数。 $$S_k=S_k-1+R_k$$

这是一个Beta 负二项分布，在您的情况下参数，使用 Wikipedia 表示法。为整数时，它也称为Beta-Pascal分布。正如您在评论中指出的那样，这是贝叶斯负二项式模型中的预测分布，在成功概率上具有共轭 Beta。$r=1$$r$

因此，您可以通过对变量进行采样，然后对负二项式变量进行采样（，即说几何分布）。$\text{Beta}(\alpha,\beta)$$u$$\text{NB}(r,u)$$r=1$

此分发在 R 包中实现brr。采样器有 name rbeta_nbinom， pmf 有 namedbeta_nbinom等。符号是，，。查看：$a=r$$c=\alpha$$d=\beta$

> 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

查看代码，可以看到它实际上调用了包的ghyper（广义超几何）分布族SuppDists：

brr::rbeta_nbinom
function(n, a, c, d){
  rghyper(n, -d, -a, c-1)
}

事实上，BNB 分布被称为IV 型广义超几何分布。ghyper请参阅包中的帮助SuppDists。我相信这也可以在 Johnson & al 的《单变量离散分布》一书中找到。