Как сделать выборку из дискретного распределения по неотрицательным целым числам?

У меня есть следующее дискретное распределение, где - известные константы:$\alpha,\beta$

$$p(x;\alpha,\beta) = \frac{\text{Beta}(\alpha+1, \beta+x)}{\text{Beta}(\alpha,\beta)} \;\;\;\;\text{for } x = 0,1,2,\dots$$

Каковы некоторые подходы для эффективной выборки из этого распределения?

Это бета-отрицательное биномиальное распределение с параметром $r=1$ в вашем случае с использованием записи Википедии. Он также называется бета-Паскаль распределения, когда $r$ является целым числом. Как вы отметили в комментарии, это прогнозное распределение в байесовской отрицательной биномиальной модели с сопряженной бета-оценкой до вероятности успеха.

Таким образом, вы можете сэмплировать его, выбрав переменную u и затем выбрав отрицательную биномиальную переменную NB ( r , u ) (с r = 1 в вашем случае, то есть с геометрическим распределением).$\text{Beta}(\alpha,\beta)$$u$$\text{NB}(r,u)$$r=1$

Это распределение осуществляется в пакете R brr. У сэмплера есть имя rbeta_nbinom, у pmf есть имя dbeta_nbinomи т. Д. Обозначения , c = α , d = β . Проверьте:$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пакете. Я полагаю, что это также можно найти в книге Джонсона и др. Univariate Discrete Distribution .

Учитывая, что уменьшается с ростомx, я предлагаю сгенерировать равномерную переменнуюu∼U(0,1)и вычислить накопленные суммыSk=k ∑ x=0Beta(α+1,β+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}$$ $x$$u\sim\mathcal{U}(0,1)$ покаSk>u.Тогда реализация будет равна соответствующемуk. Поскольку R x = бета ( α + 1 , β + x $$S_k=\sum_{x=0}^k \frac{\text{Beta}(\alpha+1, \beta+x)}{\text{Beta}(\alpha,\beta)}$$ $$S_k>u$$ $k$иSk=Sk-1+Rkв вычислениях можно полностью избежать использования гамма-функций. $$\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}}$$ $$S_k=S_k-1+R_k$$