Jak szybko próbkować X, jeśli exp (X) ~ Gamma?

Mam prosty problem z próbkowaniem, w którym moja wewnętrzna pętla wygląda następująco:

v = sample_gamma(k, a)

gdzie sample_gammapróbki z rozkładu gamma tworzą próbkę Dirichleta.

Działa dobrze, ale w przypadku niektórych wartości k / a niektóre z niższych obliczeń są niedopełnione.

Dostosowałem go do używania zmiennych przestrzeni dziennika:

v = log(sample_gamma(k, a))

Po dostosowaniu całej reszty programu działa poprawnie (przynajmniej daje mi te same dokładne wyniki w przypadkach testowych). Jest jednak wolniejszy niż wcześniej.

Czy jest sposób na bezpośrednie próbkowanie bez użycia wolnych funkcji, takich jak ? Próbowałem google, ale nawet nie wiem, czy ta dystrybucja ma wspólną nazwę (log-gamma?). $X, \exp(X) \sim \text{Gamma}$ $\log()$

sampling gamma-distribution luispedro
źródło

Wszystko, co musisz zrobić, to podzielić każdą zmienną gamma przez ich sumę. Jak zatem występuje niedopełnienie? W jaki sposób przyjęcie logarytmu rozwiązuje ten problem (i tak nie można obliczyć sumy bez ponownego wykładnika)?

whuber

@ whuber W przestrzeni dziennika obliczasz sumę, a następnie odejmujesz ją od każdego elementu. Pozwala to uniknąć pierwszego punktu niedomiaru. Jest trochę dalszego przetwarzania, gdy te dirichlety służą jako składniki mieszaniny i ponownie zostają pomnożone przez małe liczby.

luispedro,

Dodawanie dzienników jest matematycznie niepoprawne: odpowiada raczej pomnożeniu gamma niż ich dodaniu. Tak, możesz uzyskać wyniki pracy, ale na pewno nie będą miały rozkładu Dirichleta! Znowu, jaka dokładnie jest natura pierwotnego niedomiaru i jakie obliczenia wykonujesz, kiedy to nastąpi? Jakie są rzeczywiste wartości, z którymi pracujesz?

whuber

@ whuber W moim opisie mogłem trochę za bardzo uprościć. Robię wszystko i {t = gamma (a, b); suma + = t; d [i] = log (t)}; logsum = log (suma); forall i {d [i] - = logsum; }. Wcześniej było to niedopełnione, jeśli a było bardzo małe.

luispedro,

Rozumiem: dla pobliżu 0 będziesz miał kłopoty bez względu na wszystko. Ciekawy problem!

α

$\alpha$

whuber

Odpowiedzi:

Rozważ mały parametr kształtu pobliżu 0, taki jak . W zakresie między 0 i , około , a więc jest w przybliżeniu gamma PDF . Można to zintegrować z przybliżonym CDF, . Odwracając go, widzimy potęgę : ogromny wykładnik. Dla powoduje to pewną szansę na niedopełnienie (wartość podwójnej precyzji mniejsza niż , więcej lub mniej). Oto wykres szansy na niedomiar w zależności od logarytmu dziesięciu podstaw $\alpha$ $\alpha = 1/100$ $\alpha$ $e^{-\alpha}$ $1$ $x^{\alpha-1}dx / \Gamma(\alpha)$ $F_\alpha(x) = \frac{x^\alpha}{\alpha \Gamma(\alpha)}$ $1/\alpha$ $\alpha = 1/100$ $10^{-300}$ $\alpha$ :

wprowadź opis zdjęcia tutaj

Jednym z rozwiązań jest wykorzystanie tego przybliżenia do generowania zmiennych log (Gamma): w efekcie spróbuj wygenerować zmienną Gamma, a jeśli jest ona zbyt mała, wygeneruj jej logarytm z przybliżonego rozkładu mocy (jak pokazano poniżej). (Wykonuj tę czynność wielokrotnie, aż dziennik znajdzie się w zakresie niedomiaru, aby był to prawidłowy zamiennik oryginalnej wariacji niedomiaru.) W obliczeniach Dirichleta odejmij maksimum wszystkich logarytmów od każdej wartości logu: domyślnie przeskaluje się wszystkie Gamma zmienia się, więc nie wpłynie to na wartości Dirichleta. Traktuj każdy wynikowy dziennik, który jest zbyt mały (powiedzmy, mniejszy niż -100), jako dziennik prawdziwego zera. Potęguj inne dzienniki. Teraz możesz kontynuować bez niedomiaru.

To potrwa jeszcze dłużej niż wcześniej, ale przynajmniej zadziała!

Aby wygenerować przybliżony log Gamma, należy zmienić za pomocą parametru kształtu , obliczenie wstępne . Jest to łatwe, ponieważ istnieją algorytmy do obliczania wartości log Gamma bezpośrednio . Wygeneruj jednolity losowy zmiennoprzecinkowy od 0 do 1, weź jego logarytm, podziel przez i dodaj do niego $\alpha$ $C = \log(\Gamma(\alpha)) + \log(\alpha)$ $\alpha$ $C$

Ponieważ parametr skali jedynie przeskalowuje zmienną, nie ma problemu z dostosowaniem jej w tych procedurach. Nie potrzebujesz go nawet, jeśli wszystkie parametry skali są takie same.

Edytować

W innej odpowiedzi OP opisuje metodę, w której moc wariantu jednolitego (zmienna ) jest mnożona przez . Działa to, ponieważ pdf wspólnej dystrybucji tych dwóch wariantów jest równy . Aby znaleźć pdf , podstawiamy , dzielimy przez jakobian i integrujemy . Całka musi zawierać się w przedziale od do ponieważ , skąd $1/\alpha$ $B(\alpha)$ $\Gamma(\alpha+1)$ $\left(\alpha x^{\alpha-1}\right) \left(y^{\alpha}e^{-y}dy/\Gamma(\alpha+1)\right)$ $z = xy$ $y \to z/x$ $x$ $x$ $z$ $\infty$ $0 \le y \le 1$

pdf (z) = \frac{α}{Γ (α + 1)} \int_{z}^{\infty} (x^{α} / x) e^{- x} (z / x)^{α - 1} d x d z = \frac{1}{Γ (α)} z^{α - 1} e^{- z} d z,

$\text{pdf}(z)=\frac{\alpha}{\Gamma(\alpha+1)}\int_z^{\infty}{\left(x^\alpha / x\right) e^{-x} (z/x)^{\alpha-1} dx} dz = \frac{1}{\Gamma(\alpha)}z^{\alpha-1}e^{-z}dz,$

który jest pdf dystrybucji . $\Gamma(\alpha)$

Chodzi o to, że gdy , wartość narysowana z jest mało prawdopodobna, a sumując log i razy log niezależnego munduru zmieniamy będzie miał zmienną . Dziennik będzie prawdopodobnie bardzo negatywny, ale ominiemy konstrukcję jego antilogu, który spadnie w reprezentacji zmiennoprzecinkowej. $0 \lt \alpha \lt 1$ $\Gamma(\alpha+1)$ $1/\alpha$ $\Gamma(\alpha)$

Whuber
źródło

Po prostu argument, aby twoja edycja była bardziej elegancka, tak naprawdę nie musisz odwoływać się do integracji. Po prostu skorzystaj z faktu, że oraz że . Są to standardowe właściwości rozkładów beta i gamma. Ponadto, gdy mamy mniej więcej , co może być szybsze do symulacji ( ) niż ogólna zmienna losowa .

\frac{Γ (α)}{Γ (α) + Γ (1)} \sim B e t a (α, 1)

$\frac{\Gamma(\alpha)}{\Gamma(\alpha)+\Gamma(1)}\sim Beta(\alpha,1)$

Γ (α) + Γ (1) \sim Γ (α + 1)

$\Gamma(\alpha)+\Gamma(1)\sim \Gamma(\alpha+1)$

α \approx 0

$\alpha\approx 0$

y \sim e x p o (1)

$y\sim expo(1)$

- \log (u)

$-\log(u)$

Γ (α + 1)

$\Gamma(\alpha+1)$

probabilityislogic

Odpowiadam na własne pytanie, ale znalazłem całkiem dobre rozwiązanie, nawet jeśli nie do końca je rozumiem. Patrząc na kod z Biblioteki Naukowej GNU, oto, jak pobiera próbki zmiennych gamma ( rjest to generator liczb losowych, ato i to ): $\alpha$ b $\beta$

  if (a < 1)
    {
      double u = gsl_rng_uniform_pos (r);
      return gsl_ran_gamma (r, 1.0 + a, b) * pow (u, 1.0 / a);
   }

gsl_ran_gammajest funkcją, która zwraca losową próbkę gamma (więc powyższe jest wywołaniem rekurencyjnym), podczas gdy gsl_rng_uniform_poszwraca równomiernie rozłożoną liczbę w ( jest to wartość ściśle dodatnia, ponieważ gwarantuje się, że nie zwróci 0.0). $(0,1)$ _pos

Dlatego mogę pobrać dziennik ostatniego wyrażenia i użyć go

return log(gsl_ran_gamma(r, 1.0 + a, b)) + log(u)/a;

Aby dostać to, czego chciałem. Mam teraz dwa log()połączenia (ale jedno mniej pow()), ale wynik jest prawdopodobnie lepszy. Wcześniej, jak zauważył whuber, miałem coś podniesionego do potęgi , potencjalnie ogromnej liczby. Teraz w przestrzeni logicznej mnożę przez . Tak więc jest mniej prawdopodobne, że będzie niedostateczne. $1/a$ $1/a$

luispedro
źródło

Czy możesz wyjaśnić, co robią gsl_rng_uniform_pos i gsl_ran_gamma? Domyślam się, że pierwszy zwraca jednolitą losową wartość od 0 do r, a drugi jest związany z wartością Gamma (1 + a, b) - może to jest niepełna Gamma? Ogólnie rzecz biorąc, wygląda to bardzo zbliżone do przybliżenia, które zasugerowałem (z wyjątkiem tego, że podczas przeglądu jest oczywiste, że zapomniałem określić podział przez część , co jest istotne!)

α

$\alpha$

whuber

Zredagowałem swoją odpowiedź, aby teraz podać więcej szczegółów.

luispedro 16.03.11

Dziękuję: ale co to jest „r”? (Uwaga: rekursja jest ograniczona: zostanie wykonane co najwyżej jedno wywołanie rekurencyjne, ponieważ a> 0 oznacza 1,0 + a> 1.)

whuber

r jest generatorem liczb losowych (skąd otrzymujesz liczby losowe).

luispedro 16.03.11

Ach, to jest sprytne: iloczyn i niezależnej odmiany okazuje się być . Zredagowałem swoją odpowiedź, aby wskazywała na twoje rozwiązanie i wyjaśnia, dlaczego to działa.

Γ (α + 1)

$\Gamma(\alpha+1)$

B (α, 1)

$B(\alpha,1)$

Γ (α)

$\Gamma(\alpha)$

whuber