Algorytm zigguratu

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

Typowa wartość generowana przez algorytm wymaga jedynie wygenerowania jednej losowej wartości zmiennoprzecinkowej i jednego losowego indeksu tablicy, a następnie jednego przeszukiwania 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) [ potrzebne źródło ] wymaga więcej obliczeń. Niemniej jednak algorytm jest znacznie szybszy obliczeniowo [ potrzebne źródło ] niż dwie najczęściej stosowane metody generowania liczb losowych o rozkładzie normalnym, Metoda biegunowa Marsaglii i transformata Boxa-Mullera , które wymagają co najmniej jednego logarytmu i jednego obliczenia pierwiastka kwadratowego dla każdej pary wygenerowanych wartości. Ponieważ jednak algorytm zigguratu jest bardziej złożony w implementacji, najlepiej jest go stosować, gdy wymagane są duże ilości liczb losowych.

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

Algorytm Ziggurata używany do generowania przykładowych wartości o rozkładzie normalnym . (Dla uproszczenia pokazano tylko wartości dodatnie.) Różowe kropki to początkowo liczby losowe o jednolitym rozkładzie. Pożądana funkcja dystrybucji 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 wynikiem jest x przetestowany, aby zobaczyć, w który obszar wycinka się wpisuje z 3 możliwymi wynikami: 1) (lewy, pełny czarny obszar) próbka wyraźnie pod krzywą i jest przekazywana bezpośrednio do wyjścia, 2) (prawy, obszar w pionowe paski) wartość próbki może leżeć pod krzywą i musi być dalej testowany. W takim przypadku generowana jest losowa 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 i algorytm jest uruchamiany od początku.

Teoria operacji

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

Rozkład, z którego wybiera algorytm zigguratu, składa się z n regionów o równej powierzchni; n - 1 prostokątów, które pokrywają większość pożądanego rozkładu, na górze nieprostokątnej podstawy, która zawiera ogon 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ę 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 (nazwijmy ją warstwą 0) ma obszar A . Do tego dodaj prostokątną warstwę o szerokości x 1 i wysokości A / x 1 , aby miała również powierzchnię 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 ). Warstwa ta zawiera każdy punkt funkcji gęstości między y1 a y2 , ale (w przeciwieństwie do warstwy bazowej) zawiera również punkty takie jak ( x1 , y2 ) , 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 takie, że x n = 0, co oznacza, że ​​​​górna ramka, warstwa n - 1, osiąga szczyt rozkładu w (0, f (0) ) Dokładnie.

Warstwa i rozciąga się pionowo od yi do yi zawarta +1 i może być podzielona poziomo na dwa regiony: (zwykle większą) część od 0 do jest w całości xi +1 , która w pożądanym rozmieszczeniu, oraz (mała) część od x i +1 do x i , która jest zawarta tylko częściowo.

0 Pomijając na chwilę problem warstwy 0 i mając jednolite zmienne losowe U i U 1 ∈ [0,1), algorytm zigguratu można opisać jako:

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

y o niskiej rozdzielczości . Krok 3 sprawdza, czy x wyraźnie mieści się w żądanej funkcji gęstości, nie wiedząc więcej o współrzędnej y. Jeśli tak nie jest, 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ę w kroku 3 przez bardzo dużą część czasu. Jednak w przypadku warstwy wierzchniej n − 1 test ten zawsze kończy się niepowodzeniem, ponieważ x n = 0.

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 obszarze środkowym, wygeneruj fikcyjne x = A / y 1 . Spowoduje to wygenerowanie punktów o x < x 1 z prawidłową częstotliwością, aw rzadkich przypadkach, gdy wybrana zostanie warstwa 0 i x ≥ x 1 , użyj specjalnego algorytmu awaryjnego, aby wybrać losowo punkt z ogona. Ponieważ algorytm awaryjny jest używany rzadziej niż raz na tysiąc, szybkość nie jest istotna.

Zatem pełny algorytm zigguratu dla rozkładów jednostronnych to:

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

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

Algorytmy awaryjne dla ogona

Ponieważ algorytm zigguratu generuje większość danych wyjściowych bardzo szybko i wymaga algorytmu awaryjnego zawsze, gdy x > x 1 , jest on zawsze bardziej złożony niż bardziej bezpośrednia implementacja. Konkretny algorytm rezerwowy zależy od dystrybucji.

W przypadku rozkładu wykładniczego ogon wygląda tak samo jak ciało 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 wywołanie algorytmu zigguratu rekurencyjnie 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łów, test w kroku 3 prawie zawsze kończy się sukcesem.

optymalizacje

Algorytm można skutecznie wykonać z wcześniej obliczonymi tablicami x i oraz y i = f ( x i ), ale istnieją pewne modyfikacje, które czynią go jeszcze szybszym:

  • Nic w algorytmie zigguratu nie zależy od znormalizowanej funkcji rozkładu prawdopodobieństwa (całka pod krzywą równa 1), usunięcie stałych normalizujących może przyspieszyć obliczenie f ( x ).
  • 0 Większość jednolitych generatorów 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 użyć takich liczb bezpośrednio dla U .
  • 0 Podczas obliczania dwustronnych rozkładów przy użyciu dwustronnego U , jak opisano wcześniej, losową liczbę całkowitą można interpretować jako liczbę ze znakiem w zakresie [-2 31 , 2 31 - 1], a współczynnik skali 2 −31 można używany.
  • 000 Zamiast porównywać U 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 z tym. Jeśli U jest generatorem liczb losowych całkowitych, te granice mogą być wstępnie pomnożone przez 2 32 (lub 2 31 , odpowiednio), więc można 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 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 losowej nie są używane . Bity te mogą być użyte do wyboru numeru warstwy. (Zobacz odnośniki poniżej, aby uzyskać szczegółowe omówienie tego).
  • Pierwsze trzy kroki można umieścić w funkcji wbudowanej , która może wywoływać poza linią implementację rzadziej używanych 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 ja = f −1 ( y ja ) i y ja +1 = y ja + A / x ja . Powtórz n − 1 razy dla warstw zigguratu. Na końcu powinieneś mieć y n = f (0). Wystąpi błąd zaokrąglenia , ale warto sprawdzić , czy jest on akceptowalnie mały.

Podczas wypełniania wartości w tabeli po prostu załóż, że x n = 0 i y n = f (0) i zaakceptuj niewielką różnicę w obszarze warstwy n − 1 jako błąd zaokrąglenia.

Znalezienie x 1 i A

  Biorąc pod uwagę początkową (zgadnij) x 1 , potrzebujesz sposobu na obliczenie obszaru t ogona, dla którego x > x 1 . Dla rozkładu wykładniczego jest to po prostu e − x 1 , podczas gdy dla rozkładu normalnego, zakładając, że używasz nieznormalizowanego f ( x ) = e − x 2 /2 , to jest √ π /2 erfc ( x / √ 2 ). W przypadku bardziej niewygodnych rozkładów może być wymagane całkowanie numeryczne .

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

Następnie oblicz serie y i oraz x i jak wyżej. Jeśli y i > f (0) dla dowolnego i < n , to początkowe oszacowanie x 1 było zbyt niskie, co prowadziło do zbyt dużego obszaru 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 tak blisko f (0), jak to możliwe. Alternatywnie, poszukaj wartości, która sprawia, że ​​obszar najwyższej warstwy, x n −1 ( f (0) − y n −1 ), jest jak najbardziej zbliżony do pożądanej wartości A. Oszczędza to jedną ocenę f -1 ( x ) i jest właściwie warunkiem największego zainteresowania.

  • Jerzego Marsaglia ; Wai Wan Tsang (2000). „Metoda Ziggurata do generowania zmiennych losowych” . Dziennik oprogramowania statystycznego . 5 (8) . Źródło 2007-06-20 . Ten artykuł numeruje warstwy od 1, zaczynając od góry, i czyni warstwę 0 na dole szczególnym przypadkiem, podczas gdy powyższe wyjaśnienie numeruje warstwy od 0 na dole.
  • C implementacja metody zigguratu dla normalnej funkcji gęstości i wykładniczej funkcji gęstości , czyli zasadniczo kopia kodu w pracy. (Potencjalni użytkownicy powinni mieć świadomość, że ten kod C zakłada 32-bitowe liczby całkowite).
  • Implementacja AC# algorytmu zigguratu i omówienie metody.
  • Jurgen A. Doornik (2005). „Ulepszona metoda zigguratu do generowania normalnych losowych próbek” (PDF) . Nuffield College w Oksfordzie . Źródło 2007-06-20 . {{ cite journal }} : Cite journalwymaga |journal= ( pomoc ) Opisuje zagrożenia związane z użyciem najmniej znaczących bitów generatora liczb całkowitych do wybrania numeru warstwy.
  • Normalne zachowanie Cleve Moler, MathWorks, opisujące algorytm zigguratu wprowadzony w MATLAB wersja 5, 2001.
  • Ziggurat Random Normal Generator Blogs of MathWorks, opublikowane przez Cleve Moler, 18 maja 2015 r.
  •    Davida B. Thomasa; Philipa HW Leonga; Wayne Luk; John D. Villasenor (październik 2007). „Gaussowskie generatory liczb losowych” (PDF) . Ankiety komputerowe ACM . 39 (4): 11:1-38. doi : 10.1145/1287620.1287622 . ISSN 0360-0300 . S2CID 10948255 . Źródło 2009-07-27 . [Gdy] utrzymanie niezwykle wysokiej jakości statystycznej jest priorytetem, a z zastrzeżeniem tego ograniczenia pożądana jest również szybkość, metoda zigguratu będzie często najwłaściwszym wyborem. Porównanie kilku algorytmów generowania liczb losowych Gaussa .
  • Nadler, Boaz (2006). „Wady projektowe we wdrażaniu metod Ziggurat i Monty Python (i kilka uwag na temat Matlab Randn)” . arXiv : matematyka/0603058 . . Ilustruje problemy z podstawowymi jednolitymi generatorami liczb pseudolosowych oraz ich wpływ na dane wyjściowe algorytmu zigguratu.
  • Edrees, Hassan M.; Cheung, Brian; Sandora, McCullen; Nummey, Dawid; Stefan, Deian (13-16 lipca 2009). Zoptymalizowany sprzętowo algorytm Ziggurat dla szybkich generatorów liczb losowych Gaussa (PDF) . Międzynarodowa konferencja 2009 na temat inżynierii systemów i algorytmów rekonfigurowalnych. Las Vegas.
  • Marsaglia, George (wrzesień 1963). Generowanie zmiennej z ogona rozkładu normalnego (raport techniczny). Laboratoria Badań Naukowych Boeinga. Nota matematyczna nr 322, numer dostępu DTIC AD0423993. Zarchiwizowane od oryginału w dniu 10 września 2014 r. - za pośrednictwem Centrum Informacji Technicznej Obrony .
Pobrano z „ https://en.wikipedia.org/w/index.php?title=Ziggurat_algorithm&oldid=1102322020 ”