Algorytm Ziggurat - Ziggurat algorithm

Algorytm zikkurat jest algorytm dla próbek liczb pseudolosowych . Należący do klasy algorytmów próbkowania odrzucania , opiera się na źródle równomiernie rozłożonych liczb losowych, zwykle z generatora liczb pseudolosowych , a także z wstępnie obliczonych tabel. Algorytm służy do generowania wartości z monotonicznie malejącego rozkładu prawdopodobieństwa . Można ją również zastosować do symetrycznych rozkładów jednomodalnych , takich jak rozkład normalny , wybierając wartość z połowy rozkładu, a następnie losowo wybierając, z której połowy wartości uważa się, że została pobrana. Został opracowany przez George'a Marsaglia i innych w latach sześćdziesiątych.

Typowa wartość wygenerowana przez algorytm wymaga tylko wygenerowania jednej losowej wartości zmiennoprzecinkowej i jednego losowego indeksu tabeli, a następnie jednego wyszukiwania w tabeli, jednej operacji mnożenia i jednego porównania. Czasami (2,5% czasu, w przypadku rozkładu normalnego lub wykładniczego przy użyciu typowych rozmiarów tabel) wymaganych jest więcej obliczeń. Niemniej jednak algorytm jest obliczeniowo znacznie szybszy niż dwie najczęściej używane metody generowania liczb losowych o rozkładzie normalnym, metoda biegunowa Marsaglii i transformata Boxa-Mullera , które wymagają co najmniej jednego logarytmu i jednego pierwiastka kwadratowego dla każdej pary generowanych wartości . Jednak ponieważ algorytm ziggurat jest bardziej złożony w implementacji, najlepiej jest go używać, gdy wymagane są duże ilości liczb losowych.

Termin algorytm ziggurat pochodzi z artykułu Marsaglii z Wai Wan Tsang w 2000 roku; jest tak nazwany, ponieważ koncepcyjnie opiera się na pokryciu rozkładu prawdopodobieństwa prostokątnymi segmentami ułożonymi w malejącej kolejności, co daje figurę przypominającą ziggurat .

Image
Algorytm Ziggurat używany do generowania wartości próbek z rozkładem normalnym . (Dla uproszczenia pokazano tylko wartości dodatnie). Różowe kropki są początkowo liczbami losowymi o równomiernym rozkładzie. Żądana funkcja rozkładu jest najpierw dzielona na równe obszary „A”. Jedna warstwa i jest wybierana losowo przez jednolite źródło po lewej stronie. Następnie losowa wartość z górnego źródła jest mnożona przez szerokość wybranej warstwy, a wynik jest testowany x, aby zobaczyć, do którego obszaru wycinka należy z 3 możliwymi wynikami: 1) (lewy, jednolity czarny obszar) próbka wyraźnie pod krzywą i jest przekazywana bezpośrednio do wyjścia, 2) (prawy, pionowy obszar) wartość próbki może leżeć pod krzywą i musi być dalej testowana. W takim przypadku generowana jest losowa wartość y w wybranej warstwie i porównywana z f (x) . Jeśli mniej, punkt znajduje się pod krzywą i wyprowadzana jest wartość x . Jeśli nie (trzeci przypadek), wybrany punkt x jest odrzucany, a algorytm jest uruchamiany od początku.

Teoria operacji

Algorytm ziggurat to algorytm próbkowania odrzucania; losowo generuje punkt w rozkładzie nieco większym niż pożądany rozkład, a następnie sprawdza, czy wygenerowany punkt znajduje się w pożądanym rozkładzie. Jeśli nie, próbuje ponownie. Biorąc pod uwagę losowy punkt pod krzywą gęstości prawdopodobieństwa, jego współrzędna x jest liczbą losową o pożądanym rozkładzie.

Rozkład wybrany przez algorytm ziggurat składa się z n równych obszarów; n  - 1 prostokąty, które pokrywają większość pożądanego rozkładu, na górze nieprostokątnej podstawy, która zawiera koniec rozkładu.

Biorąc pod uwagę monotoniczną malejącą funkcję gęstości prawdopodobieństwa f ( x ), zdefiniowaną dla wszystkich x ≥ 0, podstawę zigguratu definiuje się jako wszystkie punkty wewnątrz rozkładu i poniżej y 1 = f ( x 1 ). Składa się on z prostokątnego obszaru od (0, 0) do ( x 1 ,  y 1 ) i (zazwyczaj nieskończonego) ogona rozkładu, gdzie x > x 1 (i y < y 1 ).

Ta warstwa (nazywają to warstwy 0) ma powierzchnię A . Do tego należy dodać prostokątną warstwę szerokość x 1 i wysokość A / x 1 , a więc ma również obszar A . Wierzchołek tej warstwy znajduje się na wysokości y 2 = y 1 + A / x 1 i przecina funkcję gęstości w punkcie ( x 2 ,  y 2 ), gdzie y 2 = f ( x 2 ). Ta warstwa zawiera każdy punkt funkcji gęstości między y 1 i y 2 , ale (w przeciwieństwie do warstwy podstawowej) zawiera również punkty takie jak ( x 1 ,  y 2 ), które nie znajdują się w pożądanym rozkładzie.

Kolejne warstwy są następnie układane na wierzchu. Aby użyć wstępnie obliczonej tabeli o rozmiarze n (  typowe jest n = 256), wybiera się x 1 taką, że x n = 0, co oznacza, że ​​górna kostka, warstwa n - 1, osiąga szczyt rozkładu w (0, f (0) ) dokładnie.

Warstwa i rozciąga się w pionie od y i do y i +1 i można ją podzielić poziomo na dwa obszary: (ogólnie większą) część od 0 do x i +1, która jest w całości zawarta w pożądanym rozkładzie, oraz (małą) część od x i +1 do x i , który jest tylko częściowo zawarty.

Ignorując na chwilę problem warstwy 0 i mając dane jednorodne zmienne losowe U 0 i U 1  ∈ [0,1), algorytm zigguratu można opisać następująco:

  1. Wybierz losową warstwę 0 ≤ i < n .
  2. Niech x = U 0 x i .
  3. Jeśli x < x i +1 , zwróć x .
  4. Niech y = y i + U 1 ( y i +1 - y i ).
  5. Oblicz f ( x ). Jeśli y < f ( x ), zwraca x .
  6. W przeciwnym razie wybierz nowe liczby losowe i wróć do kroku 1.

Krok 1 sprowadza się do wyboru współrzędnej y o niskiej rozdzielczości . Krok 3 sprawdza, czy współrzędna x wyraźnie mieści się w żądanej funkcji gęstości, nie wiedząc więcej o współrzędnej y. Jeśli nie, krok 4 wybiera współrzędną y o wysokiej rozdzielczości, a krok 5 przeprowadza test odrzucenia.

W przypadku blisko rozmieszczonych warstw algorytm kończy się na kroku 3 przez bardzo dużą część czasu. Jednak dla górnej warstwy n - 1 ten test zawsze kończy się niepowodzeniem, ponieważ x n = 0.

Warstwę 0 można również podzielić na region centralny i krawędź, ale krawędź jest nieskończonym ogonem. Aby użyć tego samego algorytmu do sprawdzenia, czy punkt znajduje się w regionie centralnym, wygeneruj fikcyjne x 0 = A / y 1 . Spowoduje to wygenerowanie punktów z x < x 1 z prawidłową częstotliwością, aw rzadkich przypadkach, gdy zostanie wybrana warstwa 0, a x ≥ x 1 , użyj specjalnego algorytmu rezerwowego, aby wybrać losowo punkt z ogona. Ponieważ algorytm awaryjny jest używany rzadziej niż raz na tysiąc, szybkość nie jest niezbędna.

Tak więc pełny algorytm zigguratu dla jednostronnych rozkładów to:

  1. Wybierz losową warstwę 0 ≤ i < n .
  2. Niech x = U 0 x i
  3. Jeśli x < x i +1 , zwróć x .
  4. Jeśli i = 0, wygeneruj punkt z końca za pomocą algorytmu rezerwowego.
  5. Niech y = y i + U 1 ( y i +1 - y i ).
  6. Oblicz f ( x ). Jeśli y < f ( x ), zwraca x .
  7. W przeciwnym razie wybierz nowe liczby losowe i wróć do kroku 1.

Oczywiście w przypadku rozkładu dwustronnego wynik musi być zanegowany w 50% przypadków. Często można to zrobić wygodnie, wybierając U 0 ∈ (−1,1) i, w kroku 3, sprawdzając, czy | x | < x i +1 .

Algorytmy rezerwowe dla ogona

Ponieważ algorytm ziggurat generuje większość wyników tylko bardzo szybko i wymaga algorytmu rezerwowego zawsze, gdy x  >  x 1 , jest on zawsze bardziej złożony niż bardziej bezpośrednia implementacja. Algorytm awaryjny zależy oczywiście od dystrybucji.

W przypadku rozkładu wykładniczego ogon wygląda tak samo, jak korpus rozkładu. Jednym ze sposobów jest powrót do najbardziej elementarnego algorytmu E  = −ln ( U 1 ) i niech x  =  x 1  - ln ( U 1 ). Innym jest rekurencyjne wywołanie algorytmu ziggurat i dodanie x 1 do wyniku.

Dla rozkładu normalnego Marsaglia sugeruje zwarty algorytm:

  1. Niech x = −ln ( U 1 ) / x 1 .
  2. Niech y = −ln ( U 2 ).
  3. Jeśli 2 y > x 2 , zwróć x  +  x 1 .
  4. W przeciwnym razie wróć do kroku 1.

Ponieważ x 1 ≈ 3,5 dla typowych rozmiarów stołu, test w kroku 3 prawie zawsze kończy się pomyślnie.

Optymalizacje

Algorytm można wykonać wydajnie z wstępnie obliczonymi tabelami x i oraz y i = f ( x i ), ale są pewne modyfikacje, aby uczynić go jeszcze szybszym:

  • Nic w algorytmie ziggurat nie zależy od znormalizowania funkcji rozkładu prawdopodobieństwa (całka pod krzywą równą 1), usunięcie stałych normalizujących może przyspieszyć obliczenia f ( x ).
  • Większość generatorów jednorodnych liczb losowych opiera się na generatorach liczb całkowitych, które zwracają liczbę całkowitą z zakresu [0, 2 32 - 1]. Tablica 2 −32 x i pozwala na użycie takich liczb bezpośrednio dla U 0 .
  • Podczas obliczania rozkładów dwustronnych przy użyciu dwustronnego U 0, jak opisano wcześniej, losową liczbę całkowitą można zinterpretować jako liczbę ze znakiem z zakresu [-2 31 , 2 31 - 1], a współczynnik skali 2 - 31 można być użytym.
  • Zamiast porównywać U 0 x i do x i +1 w kroku 3, możliwe jest wstępne obliczenie x i +1 / x i i bezpośrednie porównanie U 0 z tym. Jeśli U 0 jest generatorem liczb losowych całkowitoliczbowych, granice te mogą być przemnożone wstępnie przez 2 32 (lub odpowiednio 2 31 ), tak aby można było zastosować porównanie liczb całkowitych.
  • Dzięki powyższym dwóm zmianom tabela niezmodyfikowanych wartości x i nie jest już potrzebna i może zostać usunięta.
  • Podczas generowania wartości zmiennoprzecinkowych o pojedynczej precyzji IEEE 754 , które mają tylko 24-bitową mantysę (w tym niejawną wiodącą 1), najmniej znaczące bity 32-bitowej liczby całkowitej nie są używane. Te bity mogą być użyte do wybrania numeru warstwy. (Zobacz odniesienia poniżej, aby uzyskać szczegółowe omówienie tego).
  • Pierwsze trzy kroki można umieścić w funkcji inline , która może wywołać zewnętrzną implementację rzadziej potrzebnych kroków.

Generowanie tabel

Możliwe jest przechowywanie całej wstępnie obliczonej tabeli lub po prostu uwzględnienie wartości n , y 1 , A i implementacji f  −1 ( y ) w kodzie źródłowym i obliczenie pozostałych wartości podczas inicjalizacji generatora liczb losowych.

Jak opisano wcześniej, można znaleźć x i = f  −1 ( y i ) i y i +1 =  y i  +  A / x i . Powtórz n  - 1 razy dla warstw zigguratu. Na koniec powinieneś mieć y n  =  f (0). Oczywiście wystąpi jakiś błąd zaokrąglenia , ale użytecznym testem poczytalności jest sprawdzenie, czy jest on akceptowalnie mały.

Wypełniając wartości tabeli, po prostu załóżmy, że x n  = 0 i y n  =  f (0) i zaakceptuj niewielką różnicę w powierzchni warstwy n  - 1 jako błąd zaokrąglenia.

Znalezienie x 1 i A

Biorąc pod uwagę początkowe (przypuszczenie) x 1 , potrzebujesz sposobu na obliczenie powierzchni t ogona, dla której x > x 1 . Dla rozkładu wykładniczego, to tylko e - x 1 , zaś dla rozkładu normalnego, przy założeniu, że są za pomocą nieznormalizowanego f ( x ) = e - x 2 /2 , to jest √ gatunku / 2   ERFC ( x / √ 2 ) . W przypadku bardziej niewygodnych rozkładów może być wymagana integracja numeryczna .

Mając to w ręku, z x 1 można znaleźć y 1 = f ( x 1 ), pole t na ogonie i pole warstwy podstawowej A = x 1 y 1  +  t .

Następnie oblicz serie y i oraz x i jak powyżej. Jeżeli Y i > C (0) dla każdego i < n , to szacunkowe x 1 była zbyt mała, co prowadzi do zbyt duży obszar A . Jeśli y n < f (0), to początkowe oszacowanie x 1 było zbyt wysokie.

Biorąc to pod uwagę, użyj algorytmu znajdowania pierwiastków (takiego jak metoda bisekcji ), aby znaleźć wartość x 1, która daje y n −1 możliwie najbliżej f (0). Alternatywnie, należy zwrócić uwagę na wartość, która powoduje, że obszar w warstwie wierzchniej, x n -1 ( f (0), -  y n -1 ) jak najbliżej do żądanej wartości A , jak to możliwe. To oszczędza jednej oceny f  −1 ( x ) i jest w rzeczywistości warunkiem najbardziej interesującym.

Bibliografia

  • George Marsaglia ; Wai Wan Tsang (2000). „Metoda Ziggurata do generowania zmiennych losowych” . Journal of Statistical Software . 5 (8) . Źródło 2007-06-20 . Ten artykuł numeruje warstwy od 1, zaczynając od góry, i sprawia, że ​​warstwa 0 na dole jest przypadkiem specjalnym, podczas gdy wyjaśnienie powyżej numeruje warstwy od 0 na dole.
  • C implementacja metody ziggurat dla funkcji gęstości normalnej i funkcji gęstości wykładniczej , czyli w zasadzie kopia kodu w artykule. (Potencjalni użytkownicy powinni być świadomi, że ten kod C zakłada 32-bitowe liczby całkowite).
  • Implementacja AC # algorytmu ziggurat i omówienie metody.
  • Jurgen A. Doornik (2005). „Ulepszona metoda zigguratu do generowania normalnych próbek losowych” (PDF) . Nuffield College w Oksfordzie . Źródło 2007-06-20 . Cite Journal wymaga |journal= ( pomoc ) Opisuje zagrożenia związane z użyciem najmniej znaczących bitów generatora liczb całkowitych losowych w celu wybrania numeru warstwy.
  • Normal Behavior By Cleve Moler, MathWorks, opisujący algorytm ziggurat wprowadzony w MATLAB w wersji 5, 2001.
  • Blogi Ziggurat Random Normal Generator firmy MathWorks, opublikowane przez Cleve Moler, 18 maja 2015 r.
  • David B. Thomas; Philip HW Leong; Wayne Luk; John D. Villasenor (październik 2007). „Generatory liczb losowych Gaussa” (PDF) . Ankiety ACM Computing . 39 (4): 11: 1–38. doi : 10,1145 / 1287620,1287622 . ISSN   0360-0300 . S2CID   10948255 . Źródło 2009-07-27 . [Jeśli] utrzymanie ekstremalnie wysokiej jakości statystycznej jest priorytetem, a przy tych ograniczeniach pożądana jest również szybkość, metoda Ziggurata będzie często najwłaściwszym wyborem. Porównanie kilku algorytmów generowania liczb losowych Gaussa .
  • Nadler, Boaz (2006). „Błędy projektowe w implementacji metod Ziggurat i Monty Python (i kilka uwag na temat Matlab randn)”. arXiv : matematyka / 0603058 . . Ilustruje problemy związane z jednolitymi generatorami liczb pseudolosowych oraz ich wpływ na wyniki algorytmu ziggurat.
  • Edrees, Hassan M .; Cheung, Brian; Sandora, McCullen; Nummey, David; Stefan, Deian (13–16 lipca 2009). Zoptymalizowany sprzętowo algorytm Ziggurat dla szybkich generatorów liczb losowych Gaussa (PDF) . 2009 Międzynarodowa konferencja nt. Inżynierii systemów rekonfigurowalnych i algorytmów. Las Vegas.
  • Marsaglia, George (wrzesień 1963). Generowanie zmiennej na podstawie ogona rozkładu normalnego (raport techniczny). Boeing Scientific Research Labs. Nota matematyczna nr 322, numer dostępu DTIC AD0423993 - za pośrednictwem Obronnego Centrum Informacji Technicznej .