Jak pobierać próbki z

Chcę próbkować zgodnie z gęstością gdzie i są ściśle pozytywne. (Motywacja: może to być przydatne do próbkowania Gibbsa, gdy parametr kształtu gęstości gamma ma wcześniejszą jednolitość).

f (a) \propto \frac{c^{a} d^{a - 1}}{Γ (a)} 1_{(1, \infty)} (a)

$f(a) \propto \frac{c^a d^{a-1}}{\Gamma(a)} 1_{(1,\infty)}(a)$

c

$c$

d

$d$

Czy ktoś wie, jak łatwo pobierać próbki z tej gęstości? Może to standard i coś, o czym nie wiem?

Można myśleć głupie odrzucenia algorytmu sampliing, które bardziej lub mniej pracy (wybrać z trybu $a^*$ o $f$ , próbki $(a,u)$ od jednolite w dużym polu $[0,10a^*]\times [0,f(a^*)]$ i odrzuć, jeśli $u>f(a)$ ), ale (i) to wcale nie jest skuteczne i (ii) $f(a^*)$ będzie zbyt duży, aby komputer mógł go łatwo obsługiwać, nawet w przypadku umiarkowanie dużych $c$ i $d$ . (Zwróć uwagę, że tryb dla dużych $c$ i $d$ wynosi około $a=cd$ .)

Z góry dziękuję za wszelką pomoc!

distributions sampling gamma-distribution

— NF
źródło

+1 dobre pytanie. Nie jestem pewien, czy istnieje standardowe podejście.

— suncoolsu

Czy sprawdziłeś (pod kątem pomysłów) w „oczywistych” miejscach, takich jak np . Tekst Devroye'a ?

— kardynał

Tak, wypróbowałem już kilka pomysłów z tekstu Devroye.

stało się trudne dla mnie, aby dostać się gdziekolwiek z większość z nich, choć ... większość podejścia wydają się wymagać albo integrację (w celu odnalezienia CDF), rozkładowi na prostsze funkcji lub ograniczające przez prostszych funkcji ... ale

funkcja sprawia, że wszystkie te trudne. Jeśli ktoś ma pomysły na to, gdzie szukać podejść do tych podproblemów - np. Gdzie indziej funkcja

pojawia się w „niezbędny” sposób, jak tutaj (nie tylko jako stała normalizująca) w statystykach - to może być bardzo pomocne !

Γ (a)

$\Gamma(a)$

Γ

$\Gamma$

Γ

$\Gamma$

— NF

Istnieje ogromna różnica między przypadkiem

. Czy musisz uwzględnić oba te przypadki?

c d < 2

$c d \lt 2$

c d \geq 2

$c d \ge 2$

— whuber

To prawda - dzięki. Możemy założyć, że

c d \geq 2

$cd\geq 2$

— NF

Odpowiedzi:

Próbkowanie przy odrzuceniu będzie działało wyjątkowo dobrze, gdy i jest uzasadnione dla . $c d \ge \exp(5)$ $c d \ge \exp(2)$

Aby trochę uprościć matematykę, pozwól , napisz i zwróć na to uwagę $k = c d$ $x = a$

f (x) \propto \frac{k^{x}}{Γ (x)} d x

$f(x) \propto \frac{k^x}{\Gamma(x)} dx$

dla . Ustawienie podaje $x \ge 1$ $x = u^{3/2}$

f (u) \propto \frac{k^{u^{3 / 2}}}{Γ (u^{3 / 2})} u^{1 / 2} d u

$f(u) \propto \frac{k^{u^{3/2}}}{\Gamma(u^{3/2})} u^{1/2} du$

dla . Gdy , rozkład ten jest bardzo zbliżony do normalnego (i zbliża się, gdy staje się większy). W szczególności możesz $u \ge 1$ $k \ge \exp(5)$ $k$

Znajdź tryb numerycznie (używając np. Newtona-Raphsona). $f(u)$
Rozwiń do drugiego rzędu dotyczącego jego trybu. $\log{f(u)}$

Daje to parametry blisko przybliżonego rozkładu normalnego. Dla wysokiej dokładności ta przybliżona Normalna dominuje z wyjątkiem skrajnych ogonów. (Kiedy , może być konieczne trochę skalowanie normalnego pdf, aby zapewnić dominację.) $f(u)$ $k \lt \exp(5)$

Po wykonaniu tej wstępnej pracy dla dowolnej wartości i oszacowaniu stałej (jak opisano poniżej), uzyskanie losowej wariacji jest kwestią: $k$ $M \gt 1$

Narysuj wartość z dominującego rozkładu normalnego . $u$ $g(u)$
Jeśli lub jeśli nowy jednolity wariant przekracza , wróć do kroku 1. $u \lt 1$ $X$ $f(u)/(M g(u))$
Zestaw . $x = u^{3/2}$

Oczekiwana liczba ocen powodu rozbieżności między i jest tylko nieznacznie większa niż 1. (Pewne dodatkowe oceny wystąpią z powodu odrzucenia wariantów mniejszych niż , ale nawet gdy jest tak niskie jak częstotliwość takich wystąpienia są niewielkie). $f$ $g$ $f$ $1$ $k$ $2$

Wykres f i g dla k = 5

Wykres ten pokazuje logarytm o g i F w funkcji u o . Ponieważ wykresy są tak blisko, musimy sprawdzić ich stosunek, aby zobaczyć, co się dzieje: $k=\exp(5)$

wykres współczynnika logarytmicznego

Wyświetla log logarytm ; uwzględniono współczynnik aby upewnić się, że logarytm jest dodatni w głównej części rozkładu; to znaczy aby zapewnić, z wyjątkiem ewentualnie w obszarach o niewielkim prawdopodobieństwem. Robiąc wystarczająco dużą, możesz zagwarantować, że $\log(\exp(0.004)g(u)/f(u))$ $M = \exp(0.004)$ $Mg(u) \ge f(u)$ $M$ $M \cdot g$ dominuje we wszystkich, z wyjątkiem najbardziej ekstremalnych ogonów (które i tak praktycznie nie mają szansy zostać wybranym w symulacji). Im większy , tym częściej będą występować odrzucenia. Ponieważ rośnie, można wybrać bardzo blisko , co praktycznie nie wiąże się z żadną karą. $f$ $M$ $k$ $M$ $1$

Podobne podejście działa nawet dla , ale dość duże wartości mogą być potrzebne, gdy , ponieważ jest zauważalnie asymetryczny. Na przykład, dla , aby uzyskać względnie dokładne , musimy ustawić : $k \gt \exp(2)$ $M$ $\exp(2) \lt k \lt \exp(5)$ $f(u)$ $k = \exp(2)$ $g$ $M=1$

Wykres dla k = 2

Górna czerwona krzywa jest wykresem podczas gdy dolna niebieska krzywa jest wykresem . Odrzucenie próbkowania względem spowoduje odrzucenie około 2/3 wszystkich losowań próbnych, potrojąc wysiłek: wciąż niezły. Prawa tylna ( i $\log(\exp(1)g(u))$ $\log(f(u))$ $f$ $\exp(1)g$ $u \gt 10$ $x \gt 10^{3/2} \sim 30$ ) będzie niedostatecznie reprezentowany w próbce odrzucenia (ponieważ nie dominuje już tam ), ale ten ogon zawiera mniej niż całkowitego prawdopodobieństwa. $\exp(1)g$ $f$ $\exp(-20) \sim 10^{-9}$

Podsumowując, po początkowej próbie obliczenia trybu i oszacowaniu kwadratycznego wyrażenia szeregu mocy wokół trybu - wysiłku, który wymaga co najwyżej kilkudziesięciu ocen funkcji - można użyć próbkowania odrzucenia na oczekiwany koszt od 1 do 3 (lub więcej) ocen na wariant. Mnożnik kosztów gwałtownie spada do 1, gdy wzrasta powyżej 5. $f(u)$ $k = c d$

Nawet jeśli potrzebne jest tylko jedno losowanie z , ta metoda jest rozsądna. Sprawdza się, kiedy potrzeba wielu niezależnych losowań dla tej samej wartości , ponieważ wówczas narzuty początkowych obliczeń są amortyzowane przez wiele losowań. $f$ $k$

Uzupełnienie

@Cardinal poprosił, dość rozsądnie, o poparcie niektórych analiz machania ręką w poprzednim. W szczególności, dlaczego transformacja sprawiają, że rozkład w przybliżeniu normalny? $x = u^{3/2}$

W świetle teorii transformacji Boxa-Coxa naturalne jest poszukiwanie transformacji mocy postaci (dla stałej , miejmy nadzieję, że nie różni się zbytnio od jedności), która sprawi, że rozkład „bardziej” będzie normalny. Przypomnijmy, że wszystkie rozkłady normalne są po prostu scharakteryzowane: logarytmy ich plików pdf są czysto kwadratowe, z zerowym składnikiem liniowym i bez wyrażeń wyższego rzędu. Dlatego możemy pobrać dowolny plik pdf i porównać go z rozkładem normalnym, rozszerzając jego logarytm jako szereg mocy wokół jego (najwyższego) szczytu. Szukamy wartości która czyni (przynajmniej) trzecią $x = u^\alpha$ $\alpha$ $\alpha$ moc zanika, przynajmniej w przybliżeniu: to jest najwięcej, co możemy mieć uzasadnioną nadzieję, że osiągnie jeden wolny współczynnik. Często działa to dobrze.

Ale jak uzyskać kontrolę nad tą konkretną dystrybucją? Po dokonaniu transformacji mocy jego pdf jest

f (u) = \frac{k^{u^{α}}}{Γ (u^{α})} u^{α - 1} .

$f(u) = \frac{k^{u^{\alpha}}}{\Gamma(u^{\alpha})} u^{\alpha-1}.$

Podejmuje logarytm i używać asymptotycznej ekspansję Stirlinga z : $\log(\Gamma)$

\log (fa (u)) \approx \log (k) u^{α} + (α - 1) \log (u) - α u^{α} \log (u) + u^{α} - \log (2) π u^{α}) / 2) + do u^{- α}

$\log(f(u)) \approx \log(k) u^\alpha + (\alpha - 1)\log(u) - \alpha u^\alpha \log(u) + u^\alpha - \log(2 \pi u^\alpha)/2 + c u^{-\alpha}$

(dla małych wartości , które nie są stałe). Działa to pod warunkiem, że jest dodatnie, co zakładamy, że tak jest (w przeciwnym razie nie możemy zaniedbać pozostałej części rozszerzenia). $c$ $\alpha$

Oblicz jego trzecią pochodną (która, podzielona przez Będzie współczynnikiem trzeciej mocy w szeregu mocy) i wykorzystaj fakt, że w szczycie pierwsza pochodna musi wynosić zero. Upraszcza to znacznie trzecią pochodną, dając (w przybliżeniu, ponieważ ignorujemy pochodną ) $3!$ $u$ $c$

- \frac{1}{2} u^{- (3 + α)} α (2 α (2 α - 3) u^{2 α} + (α^{2} - 5 α + 6) u^{α} + 12 c α) .

$-\frac{1}{2} u^{-(3+\alpha)} \alpha \left(2 \alpha(2 \alpha-3) u^{2 \alpha} + (\alpha^2 - 5\alpha +6)u^\alpha + 12 c \alpha \right).$

$k$ $u$ $\alpha$ $2\alpha$

2 α - 3 = 0.

$2 \alpha-3 = 0.$

That's why $\alpha = 3/2$ works so well: with this choice, the coefficient of the cubic term around the peak behaves like $u^{-3}$ , which is close to $\exp(-2 k)$ . Once $k$ exceeds 10 or so, you can practically forget about it, and it's reasonably small even for $k$ down to 2. The higher powers, from the fourth on, play less and less of a role as $k$ gets large, because their coefficients grow proportionately smaller, too. Incidentally, the same calculations (based on the second derivative of $log(f(u))$ at its peak) show the standard deviation of this Normal approximation is slightly less than $\frac{2}{3}\exp(k/6)$ , with the error proportional to $\exp(-k/2)$ .

— whuber
źródło

(+1) Great answer. Perhaps you could expand briefly on the motivation for your choice of transformation variable.

— cardinal

Nice addition. This makes a very, very complete answer!

— cardinal

Bardzo podoba mi się odpowiedź @ whubera; prawdopodobnie będzie bardzo wydajny i ma piękną analizę. Ale wymaga to głębokiego wglądu w odniesieniu do tej konkretnej dystrybucji. W sytuacjach, w których nie masz takiego wglądu (tak dla różnych dystrybucji), podoba mi się również następujące podejście, które działa dla wszystkich dystrybucji, w których plik PDF można dwukrotnie różnicować, a ta druga pochodna ma skończenie wiele korzeni. Konfiguracja wymaga sporo pracy, ale potem masz silnik, który działa dla większości dystrybucji, które możesz na niego rzucić.

Zasadniczo chodzi o to, aby użyć częściowego liniowego górnego ograniczenia do pliku PDF, który dostosowuje się podczas próbkowania odrzucania. Jednocześnie masz kawałek liniowo niższyzwiązany z plikiem PDF, co uniemożliwia zbyt częstą jego ocenę. Górne i dolne granice są określone przez akordy i styczne do wykresu PDF. Początkowy podział na przedziały jest taki, że w każdym przedziale PDF jest albo wklęsły, albo w całości wypukły; ilekroć musisz odrzucić punkt (x, y), dzielisz ten przedział na x. (Możesz także wykonać dodatkowy podział w x, jeśli musiałbyś obliczyć PDF, ponieważ dolna granica jest naprawdę zła.) To sprawia, że podziały występują szczególnie często, gdy górna (i dolna) granica jest zła, więc masz naprawdę dobrą zbliżenie pliku PDF zasadniczo za darmo. Szczegóły są trochę skomplikowane, aby uzyskać prawo, ale starałem się wyjaśnić, z czego większość w tej serii z bloga stanowisk - zwłaszczaostatni .

Te posty nie omawiają, co zrobić, jeśli plik PDF nie jest powiązany ani w domenie, ani w wartościach; Poleciłbym dość oczywiste rozwiązanie albo przekształcenia, które uczynią je skończonymi (co byłoby trudne do zautomatyzowania), albo zastosowania odcięcia. Wybrałbym granicę odcięcia w zależności od całkowitej liczby punktów, które spodziewasz się wygenerować, powiedzmy N , i wybrałem granicę tak, aby usunięta część miała mniej niż $1 / (10 N)$ prawdopodobieństwo. (Jest to dość łatwe, jeśli masz zamknięty formularz dla CDF; w przeciwnym razie może to być trudne.)

Ta metoda jest zaimplementowana w Maple jako domyślna metoda dla ciągłych dystrybucji zdefiniowanych przez użytkownika. (Pełne ujawnienie - Pracuję dla Maplesoft.)

Zrobiłem przykładowy przebieg, generując 10 ^ 4 punktów dla c = 2, d = 3, określając [1, 100] jako początkowy zakres wartości:

wykres

Odrzucono 23 (na czerwono), 51 punktów „w okresie próbnym”, które znajdowały się w międzyczasie między dolną granicą a faktycznym plikiem PDF, i 9949 punktów, które zostały zaakceptowane po sprawdzeniu tylko nierówności liniowych. To łącznie 74 oceny pliku PDF lub około jedna ocena PDF na 135 punktów. Współczynnik powinien się poprawiać, gdy generujesz więcej punktów, ponieważ przybliżenie staje się coraz lepsze (i odwrotnie, jeśli wygenerujesz tylko kilka punktów, stosunek będzie gorszy).

— Erik P.
źródło

A tak przy okazji - jeśli musisz oceniać plik PDF bardzo rzadko, ponieważ masz do niego dobrą dolną granicę, możesz pozwolić sobie na dłużej, więc możesz po prostu użyć biblioteki bignum (może nawet MPFR?) I ocenić funkcja gamma bez nadmiernego strachu przed przepełnieniem.

— Erik P.,

(+1) This is a nice approach. Thanks for sharing it.

— whuber

Problem przepełnienia rozwiązuje się poprzez wykorzystanie (prostych) relacji między gamma. Chodzi o to, że po normalizacji szczyt ma być w pobliżu

1

$1$ , jedyne ważne obliczenia mają formę

Γ (\exp (c d)) / Γ (x)

$\Gamma(\exp(c d))/\Gamma(x)$ gdzie

x

$x$ jest dość blisko

\exp (k)

$\exp(k)$ - cała reszta będzie tak bliska zeru, że możesz je zaniedbać. Ten stosunek można uprościć do znalezienia dwóch wartości

Γ

$\Gamma$ dla argumentów pomiędzy

1

$1$ i

2

$2$ plus suma małej liczby logarytmów: nie ma tam przelewu.

— whuber

@whuber re: Gammas: Ach tak - widzę, że zasugerowałeś to również powyżej. Dzięki!

— Erik P.,

Możesz to zrobić numerycznie, wykonując metodę inwersji, która mówi, że jeśli wstawisz jednolite (0,1) zmienne losowe do odwrotnego CDF, otrzymasz losowanie z rozkładu. Poniżej zamieściłem kod R, który to robi, i po kilku sprawdzeniach, które przeprowadziłem, działa dobrze, ale jest nieco niechlujny i jestem pewien, że możesz go zoptymalizować.

Jeśli nie znasz R, lgamma () jest log funkcji gamma; całka () oblicza określoną całkę 1-D; uniroot () oblicza pierwiastek funkcji za pomocą bisekcji 1-D.

# density. using the log-gamma gives a more numerically stable return for 
# the subsequent numerical integration (will not work without this trick)
f = function(x,c,d) exp( x*log(c) + (x-1)*log(d) - lgamma(x) )

# brute force calculation of the CDF, calculating the normalizing constant numerically
F = function(x,c,d) 
{
   g = function(x) f(x,c,d)
   return( integrate(g,1,x)$val/integrate(g,1,Inf)$val )
}

# Using bisection to find where the CDF equals p, to give the inverse CDF. This works 
# since the density given in the problem corresponds to a continuous CDF. 
F_1 = function(p,c,d) 
{
   Q = function(x) F(x,c,d)-p
   return( uniroot(Q, c(1+1e-10, 1e4))$root )
}

# plug uniform(0,1)'s into the inverse CDF. Testing for c=3, d=4. 
G = function(x) F_1(x,3,4)
z = sapply(runif(1000),G)

# simulated mean
mean(z)
[1] 13.10915

# exact mean
g = function(x) f(x,3,4)
nc = integrate(g,1,Inf)$val
h = function(x) f(x,3,4)*x/nc
integrate(h,1,Inf)$val
[1] 13.00002 

# simulated second moment
mean(z^2)
[1] 183.0266

# exact second moment
g = function(x) f(x,3,4)
nc = integrate(g,1,Inf)$val
h = function(x) f(x,3,4)*(x^2)/nc
integrate(h,1,Inf)$val
[1] 181.0003

# estimated density from the sample
plot(density(z))

# true density 
s = seq(1,25,length=1000)
plot(s, f(s,3,4), type="l", lwd=3)

Główną arbitralną rzeczą, którą tu robię, jest założenie tego $(1,10000)$ jest wystarczającym przedziałem dla podziału - byłem leniwy co do tego i może być bardziej skuteczny sposób na wybranie tego przedziału. W przypadku bardzo dużych wartości obliczenia numeryczne CDF (powiedzmy $> 100000$ ) kończy się niepowodzeniem, więc przedział musi być poniżej tego. CDF jest faktycznie równy 1 w tych punktach (chyba że $c, d$ są bardzo duże), więc prawdopodobnie można by dołączyć coś, co zapobiegłoby błędnemu obliczeniu CDF dla bardzo dużych wartości wejściowych.

Edycja: kiedy $cd$ jest bardzo duży, przy tej metodzie występuje problem numeryczny. Jak zauważył Whuber w komentarzach, kiedy to się stanie, rozkład jest zasadniczo zdegenerowany w swoim trybie, co czyni go trywialnym problemem próbkowania.

— Makro
źródło

Metoda jest poprawna, ale strasznie bolesna! Jak oceniasz, ile ocen funkcji jest potrzebnych dla jednej zmiennej losowej? Tysiące? Kilkadziesiąt tysięcy?

— whuber

Jest dużo obliczeń, ale tak naprawdę nie zajmuje to dużo czasu - na pewno znacznie szybciej niż próbkowanie odrzucania. Symulacja, którą pokazałem powyżej, zajęła mniej niż minutę. Problem polega na tym, że kiedy

c d

$cd$ jest duży, wciąż się psuje. Dzieje się tak, ponieważ musi obliczyć ekwiwalent

(c d)^{x}

$(cd)^{x}$ dla dużych

x

$x$ . Jednak każde zaproponowane rozwiązanie będzie miało ten problem - staram się dowiedzieć, czy jest jakiś sposób, aby to zrobić w skali dziennika i dokonać ponownej transformacji.

— Makro

Minuta na 1000 wariantów nie jest zbyt dobra: będziesz czekał godzinami na jedną dobrą symulację Monte-Carlo. Możesz przejść cztery rzędy wielkości szybciej, stosując próbkowanie odrzucania. Sztuką jest odrzucić z dokładnym przybliżeniem

f

$f$ zamiast w odniesieniu do jednolitego rozkładu. W odniesieniu do obliczeń: oblicz

a \log (c d) - \log (Γ (a))

$a \log(c d) - \log(\Gamma(a))$ (oczywiście obliczając log Gamma bezpośrednio), a następnie potęgować. To pozwala uniknąć przepełnienia.

— whuber

That is what I do for the computation - it still doesn't avoid overflow. You can't exponentiate a number greater than around 500 on a computer. That quantity gets much larger than that. I mean "pretty good" comparing it with the rejection sampling the OP mentioned.

— Makro

I did notice that the "standard deviation rule" that normals follow (68% within 1, 95% within 2, 99.7% within 3) did apply. So basically for large

c d

$cd$ it's a point mass at the mode. From what you say, the threshold where this occurs before the numerical problems, so this still works. Thanks for the insight

— Macro