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ść).c

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!

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

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

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

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

1. Znajdź tryb numerycznie (używając np. Newtona-Raphsona).$f(u)$

2. 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 k < exp ( 5 ) , 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 M > 1 (jak opisano poniżej), uzyskanie losowej wariacji jest kwestią:$k$$M \gt 1$

1. Narysuj wartość z dominującego rozkładu normalnego g ( u ) .$u$$g(u)$

2. Jeśli lub jeśli nowy jednolity wariant X przekracza f ( u ) / ( M g ( u ) ) , wróć do kroku 1.$u \lt 1$$X$$f(u)/(M g(u))$

3. Zestaw .$x = u^{3/2}$

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

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

Wyświetla log logarytm ; uwzględniono współczynnik M = exp ( 0,004 ), aby upewnić się, że logarytm jest dodatni w głównej części rozkładu; to znaczy aby zapewnić, M g ( U ) ≥ F ( u ) z wyjątkiem ewentualnie w obszarach o niewielkim prawdopodobieństwem. Robiąc M wystarczająco dużą, możesz zagwarantować, że M ⋅ $\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 M , tym częściej będą występować odrzucenia. Ponieważ k rośnie, M można wybrać bardzo blisko 1 , 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 M mogą być potrzebne, gdy exp ( 2 ) < k < exp ( 5 ) , ponieważ f ( u ) jest zauważalnie asymetryczny. Na przykład, dla k = exp ( 2 ) , aby uzyskać względnie dokładne g , musimy ustawić M = 1 :$k \gt \exp(2)$$M$$\exp(2) \lt k \lt \exp(5)$$f(u)$$k = \exp(2)$$g$$M=1$

Górna czerwona krzywa jest wykresem podczas gdy dolna niebieska krzywa jest wykresem log ( f ( u ) ) . Odrzucenie próbkowania f względem exp ( 1 ) g spowoduje odrzucenie około 2/3 wszystkich losowań próbnych, potrojąc wysiłek: wciąż niezły. Prawa tylna ( U > 10 i x > 10 3 / 2 ~ $\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 f ), ale ten ogon zawiera mniej niż exp ( - 20 ) ∼ 10 - 9 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 k = c d 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 k , 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

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

(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 u 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ą c )$3!$$u$$c$

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

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

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:

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

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.