Zikkuratový algoritmus

Aktuální verze stránky ještě nebyla zkontrolována zkušenými přispěvateli a může se výrazně lišit od verze recenzované 21. března 2018; kontroly vyžadují 4 úpravy .

Ziggurat Algorithm ( angl. Ziggurat Algorithm , Ziggurat Method ) je algoritmus pro vzorkování pseudonáhodných čísel . Jako zástupce třídy vzorkovacích algoritmů s odchylkou se ve své práci opírá o zdroj rovnoměrně rozložených náhodných čísel - obvykle generátor pseudonáhodných čísel nebo předem vypočítanou tabulku. Algoritmus se používá ke generování hodnot na základě monotónně klesajícího rozdělení pravděpodobnosti . Může být také aplikován na symetrické unimodální rozdělení, jako je normální, výběrem hodnot z jedné její poloviny a poté v případě potřeby přechodem na symetrickou hodnotu pomocí operace aritmetické negace. Jedním z autorů algoritmu vyvinutého v 60. letech 20. století je George Marsaglia .

V nejjednodušším případě vyžaduje výpočet hodnoty vrácené algoritmem pouze generování jednoho float a jednoho indexu náhodné tabulky, po kterém následuje jedno vyhledávání v tabulce, jedno násobení a jedno srovnání. Někdy (v mnohem menším počtu případů) jsou vyžadovány složitější výpočty. Tento algoritmus je však z výpočetního hlediska mnohem rychlejší než dvě nejběžněji používané metody pro generování normálně rozdělených náhodných čísel: Marsaglia polární metoda a Box-Mullerova transformace , které vyžadují výpočet alespoň jednoho logaritmu a jednoho čtverce . kořen pro každou dvojici generovaných hodnot. Protože je však algoritmus Zikkurat složitější na implementaci, nejčastěji se používá v případech, kdy je vyžadováno velké množství náhodných čísel.

Samotný termín „Zikkuratský algoritmus“ se objevuje ve společné práci Marsaglia a Wai Van Tsanga v roce 2000 a je tak pojmenován, protože je koncepčně založen na pokrytí rozdělení pravděpodobnosti pomocí pravoúhlých segmentů naskládaných na sebe v pořadí klesající velikosti (když při pohledu zdola nahoru), výsledkem je postava připomínající zikkurat .

Teoretický základ

Algoritmus zikkuratu je algoritmus vzorkování zkreslení. Náhodně vygeneruje bod, který se mírně odchyluje od požadovaného rozložení, a poté zkontroluje, zda vygenerovaný bod spadá přesně do něj. Pokud ne, algoritmus se pokusí znovu. Leží-li bod pod křivkou funkce hustoty pravděpodobnosti, pak jeho x -ová souřadnice bude požadované náhodné číslo s požadovaným rozdělením.

Distribuce, ze které algoritmus vzorky sestává z oblastí o stejné ploše; obdélník pokrývá hlavní část požadované distribuce a je "pyramida" na nepravoúhlé základně, která zahrnuje zbytek nebo "ocas" distribuce. $n$ $n-1$

Pro danou monotónně klesající funkci hustoty pravděpodobnosti definovanou pro všechny je základ zikkuratu definován jako všechny body v distribuci a pod některými . Skládá se z pravoúhlé části od do a (obvykle nekonečného) zbytku (ocasu) rozdělení, kde (a ). $f(x)$ $x\geqslant 0$ $y_{1}=f(x_{1})$ $(0, 0)$ $(x_{1},y_{1})$ $x>x_{1}$ $y<y_{1}$

Tato úroveň (nazvěme ji úroveň 0) má plochu . K jeho vrcholu přidáme novou obdélníkovou úroveň šířky a výšky , takže jeho plocha bude také rovna . Vrchol této úrovně je ve výšce a protíná funkci hustoty v bodě , kde . Tato úroveň zahrnuje všechny body funkce hustoty mezi a , ale (na rozdíl od základní úrovně) zahrnuje i další body, jako například , které nepatří do požadované distribuce. $A$ $x_{1}$ $A/x_{1}$ $A$ ${\displaystyle y_{2}=y_{1}+A/x_{1))$ $(x_{2},y_{2})$ $y_{2}=f(x_{2})$ $y_1$ $y_2$ $(x_{1},y_{2})$

Všechny následující úrovně se překrývají stejným způsobem. Chcete-li použít předem vypočítanou tabulku velikostí ( používanou velmi často), měli byste zvolit takovou , aby horní obdélníková úroveň s číslem dosáhla vrcholu rozdělení přesně v bodě . $n$ $n=256$ $x_{1}$ $x_{n}=0$ $n-1$ $(0,f(0))$

Úroveň s číslem na výšku zaujímá místo od do a lze ji rozdělit na šířku do dvou oblastí: část od do (obvykle větší), která je celá obsažena v daném rozložení, a část od do (menší), která je uvnitř obsažena jen částečně. $i$ $y_{i}$ $y_{i+1}$ $0$ $x_{{i+1}}$ $x_{{i+1}}$ $x_{i}$

Když na chvíli zapomeneme na otázku speciálního případu s úrovní 0 a s rovnoměrně rozloženými čísly a , lze algoritmus popsat takto: $U_{0}$ $U_{1}$ $\in [0,1)$

Vyberte si náhodnou úroveň . $0\leqslant i<n$
Dejte . ${\displaystyle x=U_{0}x_{i))$
Pokud , vraťte se . $x<x_{i+1}$ $X$
Dejte . $y=y_{i}+U_{1}(y_{i+1}-y_{i})$
Vypočítejte . Pokud , vraťte se . $f(x)$ $y<f(x)$ $X$
V opačném případě vyberte nová náhodná čísla a vraťte se ke kroku 1.

Krok 1 je náhodné vzorkování úrovně. Krok 3 zkontroluje, zda souřadnice leží dobře v dané funkci hustoty i bez jakýchkoli informací o souřadnici . Pokud tomu tak není, krok 4 vypočítá souřadnici a krok 5 zkontroluje, zda je uvnitř požadované oblasti. $X$ $y$ $y$

Pokud je počet úrovní dostatečně velký a mají malou výšku, pak je stejná "riziková zóna", která se kontroluje po kroku 3, velmi malá a algoritmus se na značnou část času zastaví v kroku 3. Všimněte si, že horní úroveň však v tomto testu vždy selže, protože . $n$ $n-1$ $x_{n}=0$

Úroveň 0 lze také rozdělit na centrální a hraniční oblast, ale hraniční oblast bude obsahovat nekonečný zbytek funkce. Chcete-li použít stejný algoritmus ke kontrole, zda bod patří do centrální oblasti, vyplatí se vygenerovat figurínu . S body se souřadnicí se bude zacházet jednoduše a pro ten vzácný případ, kdy byla zvolena úroveň 0 a , budete muset použít speciální záložní algoritmus k náhodnému výběru bodu z "ocasu" funkce. Vzhledem k tomu, že takovýto záložní algoritmus bude používán extrémně zřídka (vzácnost je relativní a závisí na vrstvení), jeho rychlost nebude mít významný dopad na celkový výkon. $x_{0}=A/y_{1}$ ${\displaystyle x<x_{1))$ ${\displaystyle x\geqslant x_{1))$

Kompletní Zigguratův algoritmus pro nesymetrické rozdělení je tedy následující:

Vyberte si náhodnou úroveň . $0\leqslant i<n$
Dejte . ${\displaystyle x=U_{0}x_{i))$
Pokud , vraťte se . $x<x_{i+1}$ $X$
Pokud , vygenerujte bod z "ocasu" pomocí záložního algoritmu. $i=0$
Dejte . $y=y_{i}+U_{1}(y_{i+1}-y_{i})$
Vypočítejte . Pokud , vraťte se . $f(x)$ $y<f(x)$ $X$
V opačném případě vyberte nová náhodná čísla a vraťte se ke kroku 1.

Pro symetrickou distribuci lze výsledek samozřejmě v 50 % případů jednoduše obrátit. Často může být vhodné vygenerovat a otestovat v kroku 3 . $U_{0}\in (-1,1)$ $|x|<x_{i+1}$

Záložní algoritmy pro konec funkce

Vzhledem k tomu, že algoritmus Zikkurat generuje většinu hodnot pouze velmi rychle a vyžaduje záložní algoritmus v případech , jsou věci složitější než přímá 6kroková implementace. Záložní algoritmus závisí na dané distribuci. $x>x_{1}$

V případě exponenciálního rozdělení je ocas ve formě distribučního těla. Jedním ze způsobů je vrátit se k nejzákladnějšímu algoritmu a vložit . Dalším způsobem je rekurzivně volat algoritmus Zikkurat a přidat k výsledku. $E=-\ln(U_{1})$ $x=x_{1}-\ln(U_{1})$ $x_{1}$

V případě normální distribuce Marsaglia navrhuje kompaktní algoritmus:

Dejte . ${\displaystyle x=-\ln(U_{1})/x_{1))$
Dejte . $y=-\ln(U_{2})$
Pokud , vraťte se . $2y>x^{2}$ ${\displaystyle x+x_{1))$
V opačném případě se vraťte ke kroku 1.

Protože tabulky mají víceméně typické velikosti, test v kroku 3 téměř vždy uspěje. $x_{1}\cca 3,5$

Optimalizace

Algoritmus lze efektivně provést pomocí předpočítaných tabulek a , ale existuje několik úprav, které jej ještě více urychlí: $x_{i}$ $y_{i}=f(x_{i})$

Nic v algoritmu nezávisí na tom, zda je funkce rozdělení pravděpodobnosti normalizována (hodnota integrálu je 1), takže odstranění normalizační konstanty může urychlit výpočet . $f(x)$
Většina jednotně distribuovaných generátorů náhodných čísel je založena na generátorech náhodných celých čísel, které vracejí celé číslo z rozsahu . Tabulka obsahující vám umožní používat taková čísla přímo jako . $[0,2^{32}-1]$ ${\displaystyle 2^{-32}x_{i))$ $U_{0}$
V případě práce se symetrickými distribucemi pomocí symetrického , jak je popsáno výše, lze náhodné celé číslo interpretovat jako číslo se znaménkem v rozsahu a lze použít škálovací faktor . $U_{0}$ $[-2^{31},2^{31}-1]$ $2^{-31}$
Namísto porovnávání s v kroku 3 je možné vypočítat předem a porovnat přímo s touto hodnotou. Pokud je generátor náhodných čísel celých čísel, lze hodnoty předem vynásobit (nebo případně ), aby bylo provedeno porovnání celých čísel. ${\displaystyle U_{0}x_{i))$ $x_{{i+1}}$ ${\displaystyle x_{i+1}/x_{i))$ $U_{0}$ $U_{0}$ $2^{32}$ $2^{{31}}$
Díky dvěma výše uvedeným změnám již tabulka nezpracovaných hodnot není potřeba a lze ji smazat. $x_{i}$
V případě IEEE 754 s jednoduchou přesností generování s plovoucí desetinnou čárkou , která používá 24bitovou mantisu (včetně implicitní 1), nejsou použity nejméně významné bity 32bitového celého náhodného čísla. Tyto bity lze použít při výběru úrovně. (zde [1] je podrobně popsána podstata problematiky).

Generování tabulky

Je možné buď ponechat tabulku předem vypočítanou a úplnou, nebo pouze zahrnout hodnoty , , , a implementaci do zdrojového kódu a zbývající hodnoty vypočítat při inicializaci generátoru náhodných čísel (v závislosti na tom, co je pro nás dražší: výpočetní čas nebo paměť). $x_{i}$ $y_{i}$ $n$ $y_1$ $A$ $f^{{-1}}(y)$

Můžete najít a . Opakujte pro všechny úrovně zikkuratu. Mělo by to nakonec vyjít . $x_{i}=f^{-1}(y_{i})$ $y_{i+1}=y_{i}+A/x_{i}$ $n-1$ $y_{n}=f(0)$

Do konečného vyplňování tabulky je potřeba dát a , přičemž drobné nesrovnalosti (pokud opravdu vyšly malé) přijmete jako chyby zaokrouhlování . $x_{n}=0$ $y_{n}=f(0)$

Hledat a $x_{1}$ $A$

Pokud existuje počáteční hodnota (vypočtená, pokud ne přesně, pak přibližně), zbývá pouze vypočítat plochu ocasní části funkce, pro kterou . Můžete počítat pomocí numerických integračních metod . $x_{1}$ $t$ $x>x_{1}$

Dále je možné z oblasti ocasní části zjistit oblast základní úrovně: . $x_{1}$ $y_{1}=f(x_{1})$ $t$ $A=x_{1}y_{1}+t$

Potom se vypočítá řada a , jak je uvedeno výše. Pokud pro nějaký , pak počáteční hodnota byla příliš malá, což vedlo k velké ploše . Pokud , pak byla počáteční hodnota příliš velká. $y_{i}$ $x_{i}$ $y_{i}>f(0)$ $i<n$ $x_{1}$ $A$ $y_{n}<f(0)$ $x_{1}$

Vzhledem k výše uvedenému můžete použít numerické řešení rovnic (například metoda půlení ) k nalezení hodnoty , které se hodnota co nejvíce blíží . Alternativně lze zvážit a najít hodnoty pro oblast nejvyšší úrovně , co nejblíže požadované hodnotě . $x_{1}$ ${\displaystyle y_{n-1))$ $f(0)$ $x_{n-1}(f(0)-y_{n-1})$ $A$

Poznámky

↑ Jurgen A. Doornik. "Vylepšená metoda zikkuratu pro generování normálních náhodných vzorků" (anglicky) // Nuffield College, Oxford. - 2005. Archivováno 7. března 2016.

Literatura

George Marsaglia Metoda zikkuratu pro generování náhodných proměnných // Journal of Statistical Software . - 2000. - 7 s. - URL : webová stránka
Jurgen A. Doornik . Vylepšená metoda zikkuratu pro generování normálních náhodných vzorků. - Nuffield College, Oxford: 2005. - 21:00 - URL: práce
David B. Thomas, Philip HW Leong, Wayne Luk, John D. Villasenor . Generátory Gaussových náhodných čísel // ACM Computing Surveys. - 2007. - 38 s. - URL: práce
Boaz Nadler . Chyby v návrhu při implementaci metod Ziggurat a Monty Python (a některé poznámky k Matlab randn) // The Journal of Business. - 2006. - 16 s. - URL: práce
Edrees, Hassan M.; Cheung, Brian; Sandora, McCullen; Nummey, David; Stefan, Deian . Hardwarově optimalizovaný zikkuratový algoritmus pro vysokorychlostní generátory Gaussových náhodných čísel // Mezinárodní konference o inženýrství rekonfigurovatelných systémů a algoritmů 2009. Las Vegas. - URL: web
Marsaglia, George . Generování proměnné z konce normálního rozdělení // Technometrie. - 1964. - V. 6, č. 1. - C 101-102. - URL: web

Odkazy

Implementace C algoritmu pro normální a exponenciální hustotu funkcí je v podstatě kopie kódu z článku.
Implementace v C# a přehled samotného algoritmu.
Blogy Ziggurat Random Normal Generator společnosti MathWorks, zveřejněné Clevem Molerem, 18. května 2015.