Pytanie:
Jak pobrać próbkę z dyskretnej dystrybucji?
Barry
2013-08-21 01:40:40 UTC
view on stackexchange narkive permalink

Załóżmy, że mam rozkład rządzący możliwym wynikiem pojedynczej zmiennej losowej X. To jest coś w rodzaju [0,1, 0,4, 0,2, 0,3], gdzie X jest wartością 1, 2, 3, 4.

Czy możliwe jest próbkowanie z tego rozkładu, tj. generowanie liczb pseudolosowych dla każdego z możliwych wyników, biorąc pod uwagę prawdopodobieństwo tego wyniku. Więc jeśli chciałbym wiedzieć, jakie jest prawdopodobieństwo uzyskania 2, operacja próbkowania może zwrócić 0,34 lub coś podobnego.

Pytam o to, że próbuję zaimplementować politykę wyboru akcji dla metoda uczenia się ze wzmocnieniem oparta na pracy naukowej. Z tego, co zebrałem w artykule, autor jest w stanie próbkować rozkład przez „odwzorowanie rozkładu jednorodnego U [0,1] poprzez skumulowane funkcje gęstości prawdopodobieństwa otrzymane przez adaptacyjną całkowanie numeryczne”. Na tej podstawie sprawdza prawdopodobieństwa przejścia dla każdej próby ...

Byłbym wdzięczny za wszelkie informacje na ten temat ...

Z góry dziękuję

Istnieje wiele metod próbkowania dyskretnych rozkładów prawdopodobieństwa. Papier używa cdf (generuje uniform, $ U = u $ on (0,1), jeśli $ u <0,1 $ wyjście "1", jeśli to $ <0,1 + 0,4 $ wyjście "2" i tak dalej). Istnieją znacznie bardziej wydajne metody, jeśli problemem jest szybkość (np. Jeśli chcesz próbkować miliardy razy).
@Glen_b, czy mógłbyś wymienić bardziej wydajne metody próbkowania dyskretnego rv? To jest bardzo ciekawe.
@Riga zobacz moją odpowiedź poniżej
Tutaj jest fajny artykuł wyjaśniający „metodę aliasów”: http://www.keithschwarz.com/darts- designers-coins/
Pięć odpowiedzi:
#1
+29
jtobin
2013-08-21 02:19:57 UTC
view on stackexchange narkive permalink

Jasne. Oto funkcja R, która będzie próbkować z tej dystrybucji n razy, z zamiennikiem:

  sampleDist = function (n) {sample (x = c (1,2, 3,4), n, replace = T, prob = c (0,1, 0,4, 0,2, 0,3))} # > sampleDist (10) # [1] 4 2 2 2 2 2 4 1 2 2  

Jeśli chcesz zejść trochę niżej, możesz zobaczyć rzeczywisty algorytm używany przez sprawdzenie źródła R (napisanego w C):

  / * Próbkowanie nierównego prawdopodobieństwa ; przypadek z wymianą * n to długości p i dop. p zawiera prawdopodobieństwa, perm * zawiera rzeczywiste wyniki, a ans zawiera tablicę wartości *, które były próbkowane. * / static void ProbSampleReplace (int n, double * p, int * perm, int nans, int * ans) {double rU; int i, j; int nm1 = n - 1; / * tożsamości elementów rekordu * / for (i = 0; i < n; i ++) perm [i] = i + 1; / * sortuj prawdopodobieństwa w porządku malejącym * / revsort (p, perm, n); / * oblicz skumulowane prawdopodobieństwa * / for (i = 1; i < n; i ++) p [i] + = p [i - 1]; / * oblicz próbkę * / for (i = 0; i < nans; i ++) {rU = unif_rand (); for (j = 0; j < nm1; j ++) {if (rU < = p [j]) przerwa; } ans [i] = perm [j]; }}  
Ok, teraz rozumiem co się dzieje, bardzo dziękuję za wszystkie odpowiedzi, naprawdę mam nadzieję, że to pomoże komuś innemu. Żałuję, że nie wybrałem wszystkich prawidłowych odpowiedzi. Dziękuję wszystkim
Dlaczego musisz sortować dystrybucję dyskretną?
@PavithranIyer Ja też tego nie widziałem (lub nie widziałem).Może to próba optymalizacji: najpierw testowanie pod kątem większych prawdopodobieństw oznacza, że częściej można zatrzymać się na początku pętli.Ale wątpię, czy jest to warte kosztu pierwszego sortowania, chyba że próbkujesz bardzo często (duże „nans”).Ale z drugiej strony, jeśli spróbujesz tylko kilka razy, nie zauważysz różnicy.
#2
+16
Glen_b
2013-08-22 08:17:49 UTC
view on stackexchange narkive permalink

W odpowiedzi na pytanie w komentarzach, oto zarys kilku potencjalnie * szybszych sposobów tworzenia dyskretnych dystrybucji niż metoda cdf.

* Mówię „potencjalnie”, ponieważ w niektórych dyskretnych przypadkach dobrze zaimplementowane podejście odwrotnego CDF może być bardzo szybkie. Ogólny przypadek jest trudniejszy do wykonania szybko bez wprowadzania dodatkowych sztuczek.

W przypadku czterech różnych wyników, jak w przykładzie w pytaniu, naiwna wersja podejścia odwrotnego CDF (lub efektywnie równoważnych podejść) to w porządku; ale jeśli istnieją setki (lub tysiące lub miliony) kategorii, może to stać się powolne, nie będąc trochę mądrzejszym (z pewnością nie chcesz przeszukiwać cdf sekwencyjnie , dopóki nie znajdziesz pierwszej kategorii w którym cdf przekracza losowy uniform). Istnieje kilka szybszych podejść niż to.

[Możesz zobaczyć kilka pierwszych rzeczy, o których wspomniałem poniżej, związanych z szybszymi niż sekwencyjnymi podejściami do lokalizowania wartości za pomocą indeksu i tak jest w pewnym sensie po prostu „mądrzejsza wersja korzystania z cdf”. Można oczywiście spojrzeć na „standardowe” podejścia do rozwiązywania pokrewnych problemów, takie jak „przeszukiwanie posortowanego pliku”, i skończyć z metodami o wiele szybszymi niż sekwencyjne; jeśli możesz wywołać odpowiednie funkcje, takie standardowe podejścia mogą często być wszystkim, czego potrzebujesz.]

W każdym razie do wydajnych metod generowania z dystrybucji dyskretnych.


1 ) „Metoda tabelowa”. Zamiast $ O (n) $ dla kategorii $ n $ , po skonfigurowaniu Wersja tego w (a) (jeśli dystrybucja jest odpowiednia) to $ O (1) $ .


a) Proste podejście - zakładanie racjonalnych prawdopodobieństw (wykonane na powyższym przykładzie danych):
- skonfiguruj tablicę z 10 komórkami, zawierającą „1”, cztery „2”, dwie „3” i trzy „4”. Wypróbuj to, używając dyskretnego munduru (łatwe do zrobienia z ciągłego munduru), a otrzymasz prosty, szybki kod.

b) Przypadek bardziej złożony - nie potrzebuje „ładnych” prawdopodobieństw. Użyj komórek 2 ^ k $ , a raczej zużyjesz kilka mniej. Rozważmy na przykład następujące kwestie:

  x 0 1 2 3 4 5 6P (X = x) 0,4581 0,0032 0,1985 0,3298 0,0022 0,0080 0,0002  

( Moglibyśmy oczywiście mieć 10000 komórek i użyć poprzedniego dokładnego podejścia, ale co, jeśli te prawdopodobieństwa są, powiedzmy, irracjonalne? . Pomnóż prawdopodobieństwa przez $ 2 ^ k $ i obetnij, aby dowiedzieć się, ile komórek każdego typu potrzebujemy:

  x 0 1 2 3 4 5 6 TotP (X = x) 0,4581 0,0032 0,1985 0,3298 0,0022 0,0080 0,0002 1,0000 [256p (x)] 117 0 50 84 0 2 0 253  

Następnie ostatnie 3 komórki są w zasadzie „generuj zamiast tego z tej innej dystrybucji” (tj. p (x) - \ frac {\ lfloor 256 p (x) \ rfloor} {256} znormalizowane do pliku pmf):

  x * 0 1 2 3 4 5 6 TotP (X = x *) 0,091200 0,273067 0,272000 0,142933 0,187733 0,016000 0,017067 1,000000  

Tabelę „spillover” można wykonać dowolną rozsądną metodą (dostajesz się tutaj tylko w około 1% przypadków, nie musi to być tak szybkie). Tak więc $ \ frac {253} {256} $ czasu, w którym generujemy losowy uniform, używamy jego pierwszych 8 bitów do wybrania losowej komórki i wyprowadzamy wartość w komórka; po wstępnej konfiguracji wszystko to można zrobić bardzo szybko. W innym czasie $ \ frac {3} {256} $ trafiliśmy na komórkę z napisem „generuj z drugiej tabeli”. Prawie zawsze generujesz pojedynczy uniform na $ (0,1) $ i otrzymujesz dyskretną liczbę losową z mnożenia, obcięcia i kosztu dostępu do elementu tablicy.

2) Metoda „podniesienia histogramu do kwadratu”; jest to trochę powiązane z (1), ale każda komórka może w rzeczywistości wygenerować jedną z dwóch wartości, w zależności od (ciągłego) uniformu. Generujesz więc dyskretną wartość od 1 do n, a następnie w każdym z nich sprawdź, czy wygenerować jej główną wartość, czy drugą. Działa z ograniczonymi zmiennymi losowymi. Nie ma tabeli spillover i generalnie używa ona znacznie mniejszych tabel niż metoda (1). Zwykle jest skonfigurowany tak, że wybór 1: n wykorzystuje kilka pierwszych bitów jednolitej liczby losowej, a pozostała część określa, którą z dwóch wartości ma wyprowadzić ten bin.

Być może najprostszym sposobem opisania metody jest zrobienie tego na powyższym przykładzie:

Pomyśl o rozkładzie jako histogramie z 4 przedziałami:

original 'histogram'

Odcięliśmy szczyty najwyższych prętów i umieściliśmy je w krótszych, „prostując”. Średnia „wysokość” słupka będzie wynosić 0,25. Więc odciąliśmy 0,15 od drugiego taktu i umieściliśmy go w pierwszym i 0,05 poza czwartym i umieściliśmy go w trzecim:

'squaring off' the histogram

Zawsze można to zorganizować w w taki sposób, że żaden kosz nie kończy się na więcej niż 2 kolorach, chociaż jeden kolor może trafić do kilku pojemników.

Więc teraz wybierasz losowo jeden z 4 pojemników (wymaga 2 losowych bitów z wierzchu munduru). Następnie używasz pozostałych bitów, aby określić równomiernie rozłożone położenie w pionie i porównać z przerwą między kolorami, aby ustalić, którą z dwóch wartości wyprowadzić. Chociaż jest bardzo szybka, zwykle nie jest tak szybka jak metoda „tabelowa”.

-

Te metody można dostosować do obsługi zmiennych nieograniczonych, gdzie znowu jest „przeważnie szybka '.

Odniesienie: http://www.jstatsoft.org/v11/i03/paper

Stosunkowo powolna część z nich to tworzenie tabele wartości; są odpowiednie, gdy wiesz, co będziesz generować („musimy próbkować wartości z tej dystrybucji wiele razy w przyszłości”), zamiast próbować tworzyć to na bieżąco. „Musimy jak najszybciej pobrać próbkę miliona wartości z tego, ale nigdy nie będziemy musieli tego robić ponownie” stwarza inne priorytety; w wielu sytuacjach niektóre ze "standardowych podejść obliczeniowych" do wyszukiwania posortowanych wartości (np. szybsze wykonanie metody cdf) mogą być najlepszym wyborem.


Są jeszcze inne szybkie podejścia do generowania z dyskretnych dystrybucji. Starannie zakodowane, możesz zrobić bardzo szybkie generowanie. Na przykład:

3) metoda odrzucenia („akceptacja-odrzucenie”) może być wykonana z dyskretnymi dystrybucjami; jeśli masz dyskretną funkcję majoralizującą („obwiednię”), która jest skalowanym dyskretnym pmf, z którego możesz już wygenerować w szybki sposób, dostosowuje się ona bezpośrednio, aw niektórych przypadkach może być bardzo szybka. Bardziej ogólnie można skorzystać z możliwości generowania na podstawie ciągłych rozkładów (na przykład dyskretyzując wynik z powrotem do dyskretnej obwiedni).

Tutaj wyobraźmy sobie, że mamy jakąś dyskretną funkcję prawdopodobieństwa $ f $ , dla której nie mamy wygodnego cdf (lub inverse-cdf) - - rzeczywiście na tej ilustracji nie mieliśmy nawet stałej normalizującej, więc nasz wykres jest unormalizowany:

plot of an unnormalized unimodal discrete probability mass function on the natural numbers; it has its mode at 2 and eventually tails off approximately geometrically

Teraz musimy znaleźć wygodną do wygenerowania funkcję prawdopodobieństwa dyskretnego $ g $ , którą można pomnożyć przez stałą $ c $ i być wszędzie co najmniej tak duże jak $ f $ (musimy być pewni, że to obowiązuje dla wszystkich $ x $ wartości). To znaczy $ cg (x) \ geq f (x) $ dla wszystkich możliwych $ x $ wartości.

Czasami można łatwo zidentyfikować odpowiedni $ g $ , ale jedną z przydatnych opcji jest wybranie mieszanki oddzielnego munduru dla lewej części oraz dystrybucja z co najmniej tak ciężkim ogonem jak $ f $ po prawej stronie. Dwa całkiem dogodne opcje to rozkład geometryczny (kiedy ogon nie zmniejsza się wolniej niż wykładniczo) i coś w rodzaju dyskretyzowanego rozkładu Pareto lub dyskretyzowanego pół-Cauchy'ego, otrzymanego przez wzięcie $ \ lfloor X \ rfloor $ dla jakiejś losowej zmiennej Pareto lub pół-Cauchy'ego $ X $ (w obu przypadkach, gdy pmf spada wolniej niż wykładniczo).

(Sam geometrię można wygenerować poprzez dyskretyzację wykładniczą).

W tym przypadku dyskretny jednolity po lewej i geometryczny po prawej działają całkiem dobrze :

The previous discrete pmf with the aforementioned uniform&geometric envelope (majorizing function)

(Przypomnienie: to, co tu wykreślono, to nieznormalizowane pmf, więc oś Y nie przedstawia prawdopodobieństwa, ale coś proporcjonalnego prawdopodobieństwa)

Następnie procedura polega na zasymulowaniu proponowanej wartości $ x $ z $ g $ , symulowanie munduru, $ U $ na $ (0, cg (x)) $ i jeśli $ U<f $ , akceptując proponowany $ x $ (w przeciwnym razie odrzucając go i generując nowy proponowany $ x $ ).

Dziękuję, Glen! Podejście 2 ^ k $ jest obiecujące. Czy mógłbyś wyjaśnić metodę „podniesienia histogramu do kwadratu”, podając przykład?
Jasne, mogę spróbować, chociaż przypomnij sobie, że twój pierwotny komentarz poprosił mnie o * nazwanie * podejścia ... i tak się nazywa. W międzyczasie istnieje dobre wyjaśnienie podstawowego procesu [tutaj] (http://www.robertowor.com/csci4151/lecture3.htm). Nazywa się to również [metodą Robin Hooda] (http://www.jstatsoft.org/v11/i03/paper).
@Riga Uzupełniłem wyjaśnienie o krótki zarys pomysłu na drugi przypadek na przykładzie w pytaniu.
Dziękuję Glen za poświęcony czas! Twoje referencje i wyjaśnienia to cenna wiedza.
Cześć @Glen_b,, czy mógłbyś wyjaśnić, w jaki sposób wymyśliłeś prawdopodobieństwa dotyczące drugiego / spillover tabeli?
@greendiod To $ p (x) - \ frac {\ lfloor 256 p (x) \ rfloor} {256} $ przeskalowane do sumy 1, więc w tym przykładzie zostanie pomnożone przez 256/3.
@Glen_b Ok, widzę, wypełniasz brakujące prawdopodobieństwo i renormalizujesz się, aby wygenerować na podstawie losowego munduru powyżej [0,1 (. Przypomina mi to metodę Ziggurata dla normalnych rv. Przy okazji, jak zdobyć pierwszy8 bitów losowego munduru?
Większość języków odpowiednich do implementacji szybkiego generowania zmiennych losowych oferuje arytmetykę bitową.Więc jeśli pracujesz w C, powiedzmy (lub piszesz w języku asemblera, aby jakiś układ scalony zszedł do niskiego poziomu, co robiłem wiele razy, gdy kompilator C lub jakikolwiek inny kompilator był nieodpowiedni), to jeststandard.Jeśli próbujesz pisać w języku wysokiego poziomu, któremu z jakiegoś powodu brakuje szybkiej manipulacji bitami, prawdopodobnie nie jesteś zainteresowany optymalizacją wydajności w tym zakresie - po prostu użyj innego munduru.... ctd
ctd ... Z drugiej strony, jeśli pytasz coś w rodzaju „jak mam manipulować bitami w * tym * języku”, jest to niewłaściwe forum dla tego pytania.
@Glen_b W rzeczywistości moje pytanie jest bardziej na poziomie konceptualnym (nawet jeśli jak zawsze, diabeł tkwi w szczegółach).Gdybym nie miał operatorów na poziomie bitowym, wybrałbym losową komórkę według klasycznej $ floor (U * 256) $, która następuje po $ U ([0, ..., 255]) $ (tak jak tyzaproponowałeś sobie generowanie z tabeli „spillover”, nawet jeśli może to być powolne).Teraz $ U $ byłoby, powiedzmy, podwójnym IEEE754 (ok, straszny szczegół implementacji), czy biorąc pierwsze 8 bitów $ U $ jest wystarczająco dobre?(dlaczego? jakikolwiek link?)
@green Jednolity RNG zazwyczaj generuje liczby całkowite w 0,1 $, ..., m-1 $ (dla jakichś $ m $) zamiast liczb zmiennoprzecinkowych, a następnie skaluje je do jakiegoś formatu zmiennoprzecinkowego (takiego jak double) w $[0,1) $ (dzieląc przez $ m $).Najłatwiej jest użyć górnych bitów liczby całkowitej, a nie górnych bitów części ułamkowej liczby podwójnej (lub innego używanego formatu zmiennoprzecinkowego).Unikasz dolnych bitów, ponieważ dla niektórych RNG nie są one wystarczająco losowe.
@Glen_b Zmyliła mnie twoja propozycja generacji z tabeli spillover.Rzeczywiście potrzebowałeś tego, aby przejść od Rnd (0, m -1) do U do Rnd (0, n - 1).dziękuję za cały wkład
Zwykle $ m $ będzie ogromne w porównaniu do $ n $
Masz na myśli $ U
@Foo tak.Jeśli uważasz, że jest to niejednoznaczne lub niewystarczająco jasne, mógłbym edytować.
@Glen_b Na wykresie musi być źle oznaczona oś pionowa, dyskretna zmienna losowa nie może mieć miary większej niż 1 w żadnej pojedynczej obserwacji.Suma wszystkich wartości musi sumować się do 1 i wszystkie są nieujemne.
@Lucas To nie jest źle oznakowane - moja odpowiedź wyraźnie stwierdza, że to, co jest wykreślane, jest * unormalizowane * $ f $, gdzie mówię: „* rzeczywiście na tej ilustracji nie mieliśmy nawet stałej normalizującej, więc nasz wykres jest unormalizowany *”(tj. zmienna px jest jedynie proporcjonalna do rzeczywistych prawdopodobieństw).Jeśli zinterpretowałbyś zmienną „px” jako prawdopodobieństwo, byłoby to błędne odczytanie mojej odpowiedzi.Mógłbym powtórzyć wyjaśnienie, że zmienna px reprezentuje nieznormalizowaną wersję pmf ($ f $) ponownie gdzieś później w odpowiedzi, jak sądzę.
@Glen_b Masz rację, odpowiedź mówi, że zagęszczenia są wyskalowane w górę - błędnie odczytałem - biorąc pod uwagę, że odpowiedź jest długa, może to ułatwić jej zrozumienie, jeśli powtórzysz ją gdzieś dalej.Przepraszamy za błędną interpretację.
Dzięki;Zrobiłem to teraz - wstawiłem zdanie pod ostatni wątek (ten z funkcją majoralizującą).Wybór nazwy zmiennej (ponieważ użyłem `px` zamiast czegoś w rodzaju` qx` say) mógł wpłynąć na wrażenie, że jest to rzeczywiste prawdopodobieństwo.
#3
+7
user1893354
2013-08-21 01:51:40 UTC
view on stackexchange narkive permalink

W Pythonie możesz zrobić coś takiego

  z scipy.stats import rv_discrete x = [1,2,3,4] px = [0.1,0.4,0.2,0.3] sample = rv_discrete (values ​​= (x, px)). rvs (size = 10)  

Co daje 10 próbek z dystrybucji. Możesz to powtórzyć, a następnie znaleźć proporcje próbek, które są 2.

#4
+6
Greg Snow
2013-08-21 02:22:14 UTC
view on stackexchange narkive permalink

Tak, jest to możliwe i dość łatwe, dokładnie w jaki sposób zależy od narzędzi, których używasz.

W R będzie to sample (1: 4, n, prob = c (0.1,0.4,0.2,0.3), replace = TRUE) gdzie n to liczba wartości, które chcesz próbkować.

W narzędziach bez równoważnej funkcji nadal możesz wygenerować jednolitą wartość, a wtedy twoja RV będzie równa 1, jeśli jest poniżej 0,1, 2, jeśli jest między 0,1 i 0,5, 3, jeśli między 0,5 a 0,7, i 4, jeśli jest większa niż 0,7 (to jest idea mapowania na kumulacja).

W swoim przykładzie możesz również pobierać próbki ze zbioru zawierającego 1, 4 2, 2 3 i 3 4, aby otrzymać te same prawdopodobieństwa.

+1. Jako działający przykład równoważnej funkcji w Excelu (żeby pokazać, jakie to proste), ustaw tablicę skumulowanych sum prawdopodobieństw (np. (0,0.1,0.5,0.7,1.0)), nadaj jej nazwę (np. „cumsum”) i użyj „= MATCH (RAND (), cumsum, 1)„ do wygenerowania wartości 1,2,3,4 z prawdopodobieństwami odpowiednio 0,1, 0,4, 0,2 i 0,3. To wyraźnie pokazuje, jak ważone próbkowanie jest powiązane z wyszukiwaniem w tablicy.
Bardzo dziękuję za tę odpowiedź, naprawdę bardzo mi pomogła, mam nadzieję, że pomoże to komuś innemu. Żałuję, że nie wybrałem wszystkich prawidłowych odpowiedzi. Dziękuję wszystkim
#5
+2
Nick Cox
2013-08-21 02:40:41 UTC
view on stackexchange narkive permalink

W Stata:

W Mata użyj rdiscrete () , jak opisano na http://www.stata.com/help.cgi?mf_runiform

W samym programie Stata jest wiele sposobów. Oto jeden:

 . gen rnd = runiform (). gen y = cond (rnd < = 0.1, 1, cond (rnd < = .5, 2, cond (rnd < = .7, 3, 4)))  


To pytanie i odpowiedź zostało automatycznie przetłumaczone z języka angielskiego.Oryginalna treść jest dostępna na stackexchange, za co dziękujemy za licencję cc by-sa 3.0, w ramach której jest rozpowszechniana.
Loading...