Ziggurat algoritmus

Az oldal jelenlegi verzióját még nem ellenőrizték tapasztalt közreműködők, és jelentősen eltérhet a 2018. március 21-én felülvizsgált verziótól ; az ellenőrzések 4 szerkesztést igényelnek .

A Ziggurat algoritmus ( angol. Ziggurat Algorithm , Ziggurat Method ) egy álvéletlen számok mintavételére szolgáló algoritmus . Mivel az eltéréssel rendelkező mintavételi algoritmusok osztályának képviselője, munkájában az egyenletes eloszlású véletlenszámok forrására támaszkodik - általában egy pszeudo-véletlenszám-generátorra vagy egy előre kiszámított táblázatra. Az algoritmus monoton csökkenő valószínűségi eloszláson alapuló értékek generálására szolgál . Alkalmazható egy szimmetrikus unimodális eloszlásra is, például a normál eloszlásra, ha az egyik feléből választunk értékeket, majd ha szükséges, az aritmetikai negációs művelettel szimmetrikus értékre váltunk. Az 1960-as években kidolgozott algoritmus egyik szerzője George Marsaglia .

A legegyszerűbb esetben az algoritmus által visszaadott érték kiszámításához csak egy lebegőpontos és egy véletlenszerű táblázatindex generálása szükséges, amelyet egy táblázatkeresés, egy szorzás és egy összehasonlítás követ. Néha (jóval kisebb számú esetben) bonyolultabb számításokra van szükség. Ez az algoritmus azonban számítási szempontból sokkal gyorsabb, mint a normál eloszlású véletlenszámok generálására leggyakrabban használt két módszer: a Marsaglia poláris módszer és a Box-Muller transzformáció , amelyekhez legalább egy logaritmus és egy négyzet kiszámítása szükséges. gyökér minden egyes generált értékpárhoz. Mivel azonban a Ziggurat algoritmus megvalósítása bonyolultabb, leggyakrabban olyan esetekben használják, amikor nagyszámú véletlenszámra van szükség.

Maga a "Ziggurat algoritmus" kifejezés Marsaglia és Wai Van Tsang 2000-es közös munkájában jelenik meg, és azért kapta ezt a nevet, mert elméletileg egy valószínűségi eloszláson alapul, ahol a téglalap alakú szegmensek egymásra vannak halmozva a méret csökkenésének sorrendjében (amikor alulról felfelé nézve), ami egy zikgurátra emlékeztető alakot eredményez .

Elméleti alap

A ziggurat algoritmus egy torzítási mintavételi algoritmus. Véletlenszerűen generál egy pontot, amely kissé eltér a kívánt eloszlástól, majd ellenőrzi, hogy a generált pont pontosan beleesik-e. Ha nem, az algoritmus újra próbálkozik. Ha a pont a valószínűségi sűrűségfüggvény görbéje alatt van, akkor annak x -koordinátája lesz a kívánt véletlenszám a kívánt eloszlással.

Az eloszlás, amelyből az algoritmus mintákat vesz , egyenlő területű régiókból áll; a téglalap lefedi a kívánt eloszlás fő részét, és "piramis" egy nem téglalap alakú alapon, amely magában foglalja az eloszlás maradékát vagy "végét". $n$ $n-1$

Egy adott monotonan csökkenő valószínűségi sűrűségfüggvény esetén , amely mindegyikre van definiálva, a zikgurát alapja az eloszláson belüli és néhány alatti összes pont . Ez áll egy téglalap alakú részből -tól -ig , és egy (általában végtelen) eloszlás maradékából (farokból), ahol (és ). $f(x)$ $x\geqslant 0$ $y_{1}=f(x_{1})$ $(0, 0)$ $(x_{1},y_{1})$ $x>x_{1}$ $y<y_{1}$

Ennek a szintnek (nevezzük 0-s szintnek) a területe a . Adjunk hozzá egy új téglalap alakú szélességi és magassági szintet a tetejéhez , így a területe is egyenlő lesz . Ennek a szintnek a teteje a magasságban van , és abban a pontban metszi a sűrűségfüggvényt, ahol . Ez a szint tartalmazza az és közötti összes sűrűségfüggvény-pontot , de (az alapszinttől eltérően) más pontokat is tartalmaz, például , amelyek nem tartoznak a kívánt eloszláshoz. $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})$

Az összes következő szint ugyanúgy van egymásra rakva. Előre kiszámított ( nagyon gyakran használt) mérettáblázat használatához olyat kell választani , hogy , így a felső négyszögletes szint a számmal pontosan a pontban érje el az eloszlás csúcsát . $n$ $n=256$ $x_{1}$ $x_{n}=0$ $n-1$ $(0,f(0))$

Egy szám magasságú szint től ig terjedő helyet foglal el , és szélességében két részre osztható: egy től ig (általában nagyobb), amely teljes egészében egy adott eloszláson belül van, és egy től ig (kisebbre), amelyet csak részben tartalmaz. $én$ $y_{i}$ $y_{i+1}$ $0$ $x_{{i+1}}$ $x_{{i+1}}$ $x_{i}$

Egy pillanatra megfeledkezve a 0 szintű speciális eset kérdéséről, egyenletes eloszlású számokkal és az algoritmus a következőképpen írható le: $U_{0}$ $U_{1}$ $\in [0,1)$

Válasszon egy véletlenszerű szintet . $0\leqslant i<n$
Tedd . ${\displaystyle x=U_{0}x_{i))$
Ha vissza . $x<x_{i+1}$ $x$
Tedd . $y=y_{i}+U_{1}(y_{i+1}-y_{i})$
Számolja ki . Ha vissza . $f(x)$ $y<f(x)$ $x$
Ellenkező esetben válasszon új véletlen számokat, és térjen vissza az 1. lépéshez.

Az 1. lépés a szint véletlenszerű mintavétele. A 3. lépés azt ellenőrzi, hogy a koordináta jól esik-e az adott sűrűségfüggvényen belül, még akkor is, ha nincs információ a koordinátáról . Ha nem, a 4. lépés kiszámítja a koordinátát , az 5. lépés pedig ellenőrzi, hogy a kívánt területen belül van-e. $x$ $y$ $y$

Ha elég nagy a szintek száma, és kicsi a magasságuk, akkor ugyanaz a "kockázati zóna", amelyet a 3. lépés után ellenőrizünk, nagyon kicsi, és az algoritmus az idő jelentős részében a 3. lépésnél megáll. Vegye figyelembe, hogy a felső szint azonban mindig megbukik ebben a tesztben, mert . $n$ $n-1$ $x_{n}=0$

A 0. szint felosztható központi és határterületre is, de a határterület a függvény végtelen maradékát tartalmazza. Ha ugyanazt az algoritmust szeretnénk használni annak ellenőrzésére, hogy egy pont a központi területhez tartozik-e, érdemes egy dummyt generálni . A koordinátákkal rendelkező pontok kezelése egyszerűen történik, és abban a ritka esetben, amikor a 0 és a szintet választották , egy speciális tartalék algoritmussal kell véletlenszerűen kiválasztani egy pontot a függvény "farkából". Mivel egy ilyen tartalék algoritmust rendkívül ritkán fognak használni (a ritkaság relatív és a rétegződéstől függ), sebessége nem lesz jelentős hatással az általános teljesítményre. $x_{0}=A/y_{1}$ ${\displaystyle x<x_{1))$ ${\displaystyle x\geqslant x_{1))$

Így a nem szimmetrikus eloszlás teljes Ziggurat algoritmusa a következő:

Válasszon egy véletlenszerű szintet . $0\leqslant i<n$
Tedd . ${\displaystyle x=U_{0}x_{i))$
Ha vissza . $x<x_{i+1}$ $x$
Ha , állítson elő egy pontot a "farokból" a tartalék algoritmus segítségével. $i=0$
Tedd . $y=y_{i}+U_{1}(y_{i+1}-y_{i})$
Számolja ki . Ha vissza . $f(x)$ $y<f(x)$ $x$
Ellenkező esetben válasszon új véletlen számokat, és térjen vissza az 1. lépéshez.

Szimmetrikus eloszlás esetén az eredmény természetesen az esetek 50%-ában egyszerűen megfordítható. Gyakran kényelmes lehet a 3. lépésben generálni és tesztelni . $U_{0}\in (-1,1)$ $|x|<x_{i+1}$

Tartalékalgoritmusok egy függvény végéhez

Mivel a Ziggurat algoritmus csak nagyon gyorsan generálja az értékek nagy részét , és esetenként tartalék algoritmust igényel , a dolgok bonyolultabbak, mint a közvetlen 6 lépéses megvalósítás. A tartalék algoritmus az adott eloszlástól függ. $x>x_{1}$

Exponenciális eloszlás esetén a farok eloszlási test formájában van. Az egyik módja az, hogy visszatérünk a legelemibb algoritmushoz , és feltesszük . Egy másik módszer a Ziggurat algoritmus rekurzív meghívása, és hozzáadásával az eredményhez. $E=-\ln(U_{1})$ $x=x_{1}-\ln(U_{1})$ $x_{1}$

Normál eloszlás esetén a Marsaglia egy kompakt algoritmust javasol:

Tedd . ${\displaystyle x=-\ln(U_{1})/x_{1))$
Tedd . $y=-\ln(U_{2})$
Ha vissza . $2y>x^{2}$ ${\displaystyle x+x_{1))$
Ellenkező esetben térjen vissza az 1. lépéshez.

Mivel a táblázatok többé-kevésbé jellemző méretűek, a 3. lépésben végzett teszt szinte mindig sikeres. $x_{1}\kb. 3,5$

Optimalizálások

Az algoritmus hatékonyan elvégezhető előre kiszámított és előre kiszámított táblák használatával , de van néhány módosítás, amely még jobban felgyorsítja: $x_{i}$ $y_{i}=f(x_{i})$

Az algoritmusban semmi nem függ attól, hogy a valószínűségi eloszlásfüggvény normalizálva van-e (az integrál értéke 1), így a normalizáló állandó eltávolítása felgyorsíthatja a kiszámítását . $f(x)$
A legtöbb egyenletes eloszlású véletlenszám-generátor véletlen egész szám generátorokon alapul, amelyek egy egész számot adnak vissza egy tartományból . A táblázatot tartalmazó táblázat lehetővé teszi az olyan számok közvetlen használatát, mint a . $[0.2^{32}-1]$ ${\displaystyle 2^{-32}x_{i))$ $U_{0}$
Abban az esetben, ha szimmetrikus eloszlásokkal dolgozunk a fent leírtak szerint, egy véletlenszerű egész szám értelmezhető előjeles számként a tartományban , és használható egy skálázási tényező . $U_{0}$ $[-2^{31},2^{31}-1]$ $2^{-31}$
A 3. lépésben szereplővel való összehasonlítás helyett lehetőség van előre kiszámítani és közvetlenül ezzel az értékkel összehasonlítani. Ha egy egész véletlenszám-generátor, akkor az értékeket előre meg lehet szorozni -val (vagy adott esetben), így egész számokat kell összehasonlítani. ${\displaystyle U_{0}x_{i))$ $x_{{i+1}}$ ${\displaystyle x_{i+1}/x_{i))$ $U_{0}$ $U_{0}$ $2^{32}$ $2^{{31}}$
A fenti két változtatással a nyersértéktáblázatra már nincs szükség, és törölhető. $x_{i}$
Az IEEE 754 egypontos lebegőpontos generálása esetén, amely 24 bites mantisszát használ (beleértve az implicit 1-et is), a 32 bites véletlenszám legkisebb jelentőségű bitjei nem kerülnek felhasználásra. Ezek a bitek szintkiválasztásnál használhatók. (itt [1] a kérdés lényegét ismertetjük részletesen).

Táblázat generálás

Lehetőség van arra, hogy a táblázat előre kiszámított és teljes legyen, vagy csak az értékeket , , , és a megvalósítást belefoglalja a forráskódba , és kiszámítja a fennmaradó értékeket a véletlenszám-generátor inicializálása során (attól függően, hogy mi az drágább számunkra: számítási idő vagy memória). $x_{i}$ $y_{i}$ $n$ $y_1$ $A$ $f^{{-1}}(y)$

Megtalálhatja és . Ismételje meg a zikgurat minden szintjén. A végén sikerülnie kell . $x_{i}=f^{-1}(y_{i})$ $y_{i+1}=y_{i}+A/x_{i}$ $n-1$ $y_{n}=f(0)$

A táblázat végső kitöltésekor a és -t kell beírni, elfogadva a kis inkonzisztenciákat (ha valóban kicsinek tűntek) kerekítési hibaként . $x_{n}=0$ $y_{n}=f(0)$

Keresés és $x_{1}$ $A$

Ha van kezdeti érték (kiszámítva, ha nem pontosan, akkor hozzávetőlegesen), akkor csak a függvény farok részének területét kell kiszámítani, amelyre . Számíthat numerikus integrációs módszerekkel . $x_{1}$ $t$ $x>x_{1}$

Továbbá a farokszakasz területéből kiolvasható az alapszint területe: . $x_{1}$ $y_{1}=f(x_{1})$ $t$ $A=x_{1}y_{1}+t$

Ezután a sorozat és a fentiek szerint kerül kiszámításra. Ha valamelyik esetén, akkor a kezdeti érték túl kicsi volt, ami nagy területet eredményezett . Ha , akkor a kezdeti érték túl nagy volt. $y_{i}$ $x_{i}$ $y_{i}>f(0)$ $i<n$ $x_{1}$ $A$ $y_{n}<f(0)$ $x_{1}$

A fentiek alapján az egyenletek numerikus megoldásával (például a felezési módszerrel ) kereshet olyan értéket , amelynél az érték a lehető legközelebb van . Alternatív megoldásként megfontolhatja és megtalálhatja a legfelső szintű terület értékeit, a lehető legközelebb a kívánt értékhez . $x_{1}$ ${\displaystyle y_{n-1))$ $f(0)$ $x_{n-1}(f(0)-y_{n-1})$ $A$

Jegyzetek

↑ Jurgen A. Doornik. "An Improved Ziggurat Method Generate Normal Random Samples" (angol) // Nuffield College, Oxford. - 2005. Archiválva : 2016. március 7.

Irodalom

George Marsaglia A Ziggurat módszer véletlen változók generálására // Journal of Statistical Software . - 2000. - 7 p. - URL : weboldal
Jürgen A. Doornik . Továbbfejlesztett Ziggurat-módszer normál véletlenszerű minták generálására. - Nuffield College, Oxford: 2005. - 9 p. - URL: állás
David B. Thomas, Philip HW Leong, Wayne Luk, John D. Villasenor . Gauss-féle véletlenszám-generátorok // ACM Computing Surveys. - 2007. - 38 p. - URL: állás
Boaz Nadler . Tervezési hibák a Ziggurat és a Monty Python módszerek megvalósításában (és néhány megjegyzés a Matlab randn-ról) // The Journal of Business. - 2006. - 16 p. - URL: állás
Edrees, Hassan M.; Cheung, Brian; Sandora, McCullen; Nummey, David; Stefan, Deian . Hardverre optimalizált Ziggurat algoritmus nagysebességű Gauss-féle véletlenszám-generátorokhoz // 2009-es Nemzetközi Konferencia az Újrakonfigurálható Rendszerek és Algoritmusok Tervezéséről. Las Vegas. - URL: webhely
Marsaglia, George . Változó generálása a normál eloszlás farkából // Technometria. - 1964. - V. 6., 1. sz. - C 101-102. - URL: webhely

Linkek

A normál és exponenciális sűrűségfüggvények algoritmusának C megvalósítása lényegében a cikk kódjának másolata.
Megvalósítás C# nyelven és magának az algoritmusnak az áttekintése.
A MathWorks Ziggurat Random Normal Generator Blogjai, közzétette: Cleve Moler, 2015. május 18.