Jak próbkować z dyskretnego rozkładu na liczbach całkowitych nieujemnych?

10

Mam następujący dyskretny rozkład, w którym są znanymi stałymi: $\alpha,\beta$

p (x; α, β) = \frac{Beta (α + 1, β + x)}{Beta (α, β)} dla x = 0, 1, 2), \dots

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

Jakie są podejścia do skutecznego próbkowania z tej dystrybucji?

— jII
źródło

9

Jest to rozkład dwumianowy Beta-ujemny z parametrem w twoim przypadku, z wykorzystaniem notacji Wikipedia. Nazwano także rozkładem Beta-Pascal, gdy jest liczbą całkowitą. Jak zauważyłeś w komentarzu, jest to rozkład predykcyjny w bayesowskim ujemnym modelu dwumianowym ze sprzężoną Beta przed prawdopodobieństwem sukcesu. $r=1$ $r$

W ten sposób możesz próbkować go, próbkując zmienną a następnie próbkując ujemną zmienną dwumianową (przy twoim przypadku , to znaczy powiedzieć rozkład geometryczny). $\text{Beta}(\alpha,\beta)$ $u$ $\text{NB}(r,u)$ $r=1$

Rozkład ten jest realizowany w pakiecie R brr. Próbnik ma nazwę rbeta_nbinom, pmf ma nazwę dbeta_nbinomitp. Notacje to , , . Czek: $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

Patrząc na kod, widać, że faktycznie wywołuje on ghyper(uogólnioną hipergeometryczną) rodzinę dystrybucji SuppDistspakietu:

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

Innymi słowy, rozkład BNB jest znany jako uogólniony rozkład hipergeometryczny typu IV . Zobacz pomoc ghyperw SuppDistspakiecie. Wierzę, że można to również znaleźć w książce Johnson & al Univariate Discrete Distribution .

— Stéphane Laurent
źródło

Ta odpowiedź jest świetna, ale byłoby jeszcze lepiej, gdybyś udowodnił, że zaksięgowana gęstość OP jest taka sama jak ujemna gęstość dwumianowa.

— Sycorax mówi Przywróć Monikę

1

@ user777 Myślę, że autor OP sam to udowodnił, biorąc pod uwagę jego komentarz do odpowiedzi Xiana (tylna predykcyjna dystrybucja w ujemnym modelu dwumianowym z koniugatem Beta przed).

— Stéphane Laurent,

10

Jeśli się uwzględni

\frac{Beta (α + 1, β + x)}{Beta (α, β)} = \frac{α}{α + β + x} \frac{β + x - 1}{α + β + x - 1} \dots \frac{β}{α + β}

$\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}$ zmniejsza się z

x

$x$ , Proponuję wygenerować zmienną jednolitą

u \sim U (0, 1)

$u\sim\mathcal{U}(0,1)$ i obliczanie skumulowanych kwot

{S.}_{k} = \sum_{x = 0}^{k} \frac{Beta (α + 1, β + x)}{Beta (α, β)}

$S_k=\sum_{x=0}^k \frac{\text{Beta}(\alpha+1, \beta+x)}{\text{Beta}(\alpha,\beta)}$ aż do

{S.}_{k} > u

$S_k>u$ Realizacja jest wtedy równa odpowiadającej

k

$k$ . Od

\begin{aligned} R_{x} & = \frac{Beta (α + 1, β + x)}{Beta (α, β)} \\ = \frac{α}{α + β + x} \frac{β + x - 1}{α + β + x - 1} \dots \frac{β}{α + β} \\ = \frac{α + β + x - 1}{α + β + x} \frac{β + x - 1}{α + β + x - 1} R_{x - 1} \\ = \frac{β + x - 1}{α + β + x} R_{x - 1} \end{aligned}

$\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}}$ i

{S.}_{k} = {S.}_{k} - 1 + R_{k}

$S_k=S_k-1+R_k$ w obliczeniach można całkowicie uniknąć korzystania z funkcji gamma.

— Xi'an
źródło

1

(+1) Korzystanie

S_{k} = 1 - \frac{Γ (a + b) Γ (b + k + 1)}{Γ (b) Γ (a + b + k + 1)}

$S_k = 1-\frac{\Gamma (a+b) \Gamma (b+k+1)}{\Gamma (b) \Gamma (a+b+k+1)}$ znacznie przyspieszy pracę.

— whuber

1

Re edycja: Podejrzewam, że wykorzystanie funkcji gamma będzie jednak pomocne w rozwiązywaniu problemu

k

$k$ pod względem

u

$u$ ,

α

$\alpha$ , i

β

$\beta$ . Na przykład można znaleźć wstępne przybliżenie do

u

$u$ za pomocą wzoru Stirlinga w ocenie

Γ (b + k + 1)

$\Gamma(b+k+1)$ i

Γ (a + b + k + 1)

$\Gamma(a+b+k+1)$ a następnie dopracowując to za pomocą kilku kroków Newtona-Raphsona. Potrzebują one oceny log Gamma i jej pochodnej. Oczywiście jeśli

α

$\alpha$ i

β

$\beta$ oba są integralne, to rozwiązanie jest źródłem wielomianu - ale nawet wtedy używanie Gammy może być nadal dobrym rozwiązaniem.

— whuber

1

Świetna odpowiedź! Zaakceptowałem odpowiedź udzieloną przez SL, ponieważ zwróciło mi ona uwagę na kluczową kwestię (nie jest częścią pierwotnego pytania), że próbkowanie z predykcji tylnej jest równoważne próbkowaniu parametru z tylnej, a następnie próbkowanie danych z prawdopodobieństwa. W szczególności powyższa funkcja rozkładu jest tylną predykcją danych geometrycznych z Beta przed parametrem

p

$p$ .

— jII