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()$

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α = 1 / 100 α e - α 1 x α - 1 d x / Γ ( α ) F a ( x ) = x $\alpha$$\alpha = 1/100$$\alpha$$e^{-\alpha}$$1$$x^{\alpha-1}dx / \Gamma(\alpha)$ 1/αα=1/10010-300$F_\alpha(x) = \frac{x^\alpha}{\alpha \Gamma(\alpha)}$$1/\alpha$$\alpha = 1/100$$10^{-300}$$\alpha$ :

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 niegoC = log ( Γ ( α ) ) + log ( α ) α $\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$

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)$

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 /$1/a$$1/a$