Skutecznie próbkuje progową dystrybucję Beta

10

Jak powinienem efektywnie próbkować z następującej dystrybucji?

xB(α,β), x>k

Jeśli k nie jest zbyt duże, próbowanie odrzucenia może być najlepszym podejściem, ale nie jestem pewien, jak postępować, gdy k jest duże. Być może istnieje jakieś asymptotyczne przybliżenie, które można zastosować?

użytkownik1502040
źródło
1
Nie jest jednoznacznie jasne, co tam zamierzasz przez „ ”. Czy masz na myśli skróconą dystrybucję beta (skróconą po lewej stronie o k )? xB(α,β), x>kk
Glen_b
@Glen_b dokładnie.
user1502040
5
Dla obu parametrów kształtu większych niż 1 rozkład beta jest wklęsły, więc obwiednie wykładnicze można wykorzystać do prób odrzucenia. Aby wygenerować nieskalowane warianty beta, próbujesz już ze skróconych rozkładów wykładniczych (co jest łatwe do zrobienia), dostosowanie tej metody powinno być proste.
Scortchi - Przywróć Monikę

Odpowiedzi:

14

Najprostszym i najogólniejszym sposobem, który stosuje się do dowolnego skróconego rozkładu (można go również uogólnić na obcinanie po obu stronach), jest zastosowanie odwrotnego próbkowania transformacji . Jeśli jest skumulowanym rozkładem odsetek, ustaw p 0 = F ( k ) i weźFp0=F(k)

UU(p0,1)X=F1(U)

gdzie jest próbką z F ściętą w lewo o k . Funkcja kwantyl F - 1 mapuje prawdopodobieństw próbek z F . Ponieważ pobieramy wartości U tylko z „obszaru”, który odpowiada wartościom rozkładu beta z regionu obciętego, spróbujesz tylko tych wartości.XFkF1FU

Metodę tę zilustrowano na poniższym obrazie, na którym obcięty obszar jest zaznaczony szarym prostokątem, punkty w kolorze czerwonym są rysowane z rozkładu a następnie przekształcane w próbki B ( 2 , 8 ) .U(p0,1)B(2,8)

Próbkowanie z transformacją odwrotną z rozkładu obciętego

Tim
źródło
5
(+1) Warto zauważyć, że funkcja kwantylu nie jest tak łatwa do oszacowania.
Scortchi - Przywróć Monikę
1
@Scortchi Jeśli a lub b mają wartość 1 lub przynajmniej liczbę całkowitą, istnieje niezła forma (patrz wikipedia ). A w Pythonie jest scipy.special.betaincodwrotność, aw R. jest pbeta.
Graipher,
3
@Graipher: Powinienem był powiedzieć „ogólnie tanio” - lepiej byłoby unikać Newtona-Raphsona lub innych iteracyjnych rozwiązań, jeśli to możliwe. (BTW to qbetafunkcja kwantyli w R.)
Scortchi - Przywróć Monikę
1
@Scortchi masz rację, ale w większości przypadków w przypadku współczesnych komputerów nie powinno to stanowić poważnego problemu. Polecam również to podejście, ponieważ jest ono bezpośrednio dostępne w większości programów i może być uogólnione na dowolną skróconą dystrybucję, tylko jeśli ktoś ma dostęp do funkcji kwantylowej.
Tim
1
Niewątpliwie dobrze jest mieć powszechnie stosowaną, łatwą do wdrożenia metodę, która ma rękę, której czas działania nie rośnie wraz z ; & dla dystrybucji z zamkniętymi funkcjami kwantylowymi, np. Weibull, musi być tak dobry, jak to tylko możliwe. Niemniej jednak podejrzewam, że k będzie musiał zostać ustawiony, aby odciąć dość dużą część rozkładu beta, zanim osiągnie efektywne algorytmy odrzucania i próbkowania, które są również dostępne w większości programów i które opierają się tylko na obliczeniu gęstości prawdopodobieństwa beta. kk
Scortchi - Przywróć Monikę
8

@ Odpowiedź Tima pokazuje, w jaki sposób można dopasować próbkowanie odwrotnej transformaty do rozkładów skróconych, uwalniając czas działania zależności od progu . Dalszą wydajność można uzyskać, unikając kosztownej oceny liczbowej funkcji kwantylu beta i stosując próbkowanie z transformacją odwrotną jako część prób odrzucania.k

Funkcja gęstości rozkładu beta o parametrach kształtu i β podwójnie obciętych przy k 1 < k 2 (dla nieco większej ogólności) wynosiαβk1<k2

f(x)=x(α1)(1x)(β1)B(k2,α,β)B(k1,α,β)

Weź dowolną monotonicznie rosnącą część gęstości między i x U : dla α , β > 1 jest to log-wklęsły, abyś mógł ją objąć funkcją wykładniczą narysowaną stycznie do dowolnego punktu wzdłuż niej:xLxUα,β>1

g(x)=cλeλ(xxL)

Znajdź , ustawiając gradienty gęstości logarytmicznej równeλ

i znajdźc, obliczając, ile gęstości wykładniczej należy przeskalować, aby osiągnąć gęstość w tym punkcie c=f(x)

λ=a1xb11x
c
c=f(x)λeλ(xxL)

wprowadź opis zdjęcia tutaj

Najlepsza obwiednia do prób odrzucenia to ta z najmniejszym obszarem poniżej. Obszar to Podstawiając wyrażenia w x na λ

A=c(1eλ(xUxL))
xλc , a porzucenie stałych współczynników daje

Q(x)=xza(1-x)b(za+b-2))x-za+1[exp((b-1)(x-xL.)1-x+xL.(za-1)x-(za-1))-exp((b-1)(x-xU)1-x+xU(za-1)x-(za-1))]

reQrexxreQrex=0

k1k2)U-log(1-U)λλ

wprowadź opis zdjęcia tutaj

Piękno tego podejścia polega na tym, że cała ciężka praca jest w konfiguracji. Po zdefiniowaniu funkcji obwiedni, obliczonej stałej normalizującej dla obciętej gęstości beta, pozostaje wygenerować jednolite losowe zmienne i wykonać na nich kilka prostych operacji arytmetycznych, logów i mocy oraz porównań. Zaostrzenie funkcji obwiedni - za pomocą linii poziomych lub więcej krzywych wykładniczych - może oczywiście zmniejszyć liczbę odrzuceń.

Scortchi - Przywróć Monikę
źródło
1
+1 Fajny pomysł. Ponieważ Beta jest w przybliżeniu normalna w przypadku niewielkich lub dużych wartości jej parametrów, w zależności od tego, jak blisko są, użycie obwiedni Gaussa może być jeszcze bardziej wydajne.
whuber
@ whuber: Chciałem rozkładu z kopertą kwantową o zamkniętej formie dla obwiedni, ale przypuszczam, że można efektywnie wygenerować skrócone zmienne Gaussa przy użyciu jednego z tych dobrych przybliżeń funkcji kwantylu Gaussa. Nadal jestem zainteresowany tym, co byś zrobił, kiedyα<1 lub β<1.
Scortchi - Przywróć Monikę
1
Dla takich małych α i β, chcesz przejść na wykładniczy ogon. Nie jestem jednak pewien, co rozumiesz przez „formę zamkniętą”: kiedy przyjrzysz się komputerowym implementacjom wykładniczym, gaussowskim, funkcjom hipergeometrycznym itp. - to znaczy wszystkim funkcjom niealgebraicznym - odkryjesz, że żadna z nich mają one ogólną „formę zamkniętą”: są obliczane przez kolejne aproksymacje, takie jak szereg Taylora, ułamki częściowe lub asymptotyczne rozszerzenia. Nie ma dużej różnicy między obliczeniem wykładniczym a obliczeniem kwantyla Gaussa.
whuber
@whuber: (1) Podejście do konstruowania kopert nie zadziałałoby, ponieważ gęstości nie są wklęsłe. (2) (a) Miałem na pewno funkcje algebraiczne + logi i potęgi, trig. funkcje, gdybym został poproszony, a może nawet funkcje gamma - przyznaję, że nie miałem dokładnego pojęcia. (b) Zdobyte punkty - szybkie oceny funkcji nie ograniczają się do tych z zamkniętymi formularzami.
Scortchi - Przywróć Monikę
1
Dobra uwaga na temat wklęsłości kłody. Podejrzewam, że podział prawa energetycznego powinien stanowić dobrą kopertęα<1 lub β<1.
whuber