# Algoritmo ziggurat

L'algoritmo ziggurat è un algoritmo di rejection sampling utilizzato per il campionamento numerico pseudo-casuale. È stato sviluppato negli anni '60 da George Marsaglia ed è utilizzato per generare valori con una distribuzione di probabilità decrescente monotona partendo da una sorgente di numeri casuali uniformemente distribuiti (solitamente un generatore di numeri pseudo-casuali) e da una tabella precalcolata. Può anche essere applicato a distribuzioni unimodali simmetriche, come la distribuzione normale.
L'algoritmo ziggurat è più veloce rispetto ad altri metodi di generazione di numeri casuali normalmente utilizzati, come ad esempio il metodo polare di Marsaglia e la trasformazione di Box-Muller. Tuttavia, poiché questo algoritmo è anche più complesso da implementare, solitamente è utilizzato solo a fronte di una grande richiesta di numeri pseudo-casuali.
Il termine ziggurat è stato utilizzato per la prima volta in un articolo di George Marsaglia e Wai Wan Tsang nel 2000. Lo chiamarono così per via del modo in cui distribuisce la probabilità: la distribuzione ricorda la forma di una ziggurat a causa dei segmenti rettangolari impilati in ordine decrescente da cui è formata.

## Teoria del funzionamento

L'algoritmo ziggurat è basato sulla tecnica del rejection sampling: genera casualmente un punto in una distribuzione più grande di quella desiderata e controlla quindi se il punto generato è valido o meno.
Dato un punto casuale al di sotto di una curva di densità di probabilità, la sua coordinata x è un numero casuale con la distribuzione desiderata.
La distribuzione scelta dall'algoritmo ziggurat è composta da n regioni di uguale area.
Data una funzione di densità di probabilità decrescente monotona y = f (x), definita per tutti i valori x ≥ 0, la base della ziggurat è definita come tutti i punti interni alla distribuzione e sottostanti la retta y$_1$ = f (x$_1$). Questa consiste in una regione rettangolare individuata dai punti (0,0) a (x$_1$,y$_1$), più la coda (tipicamente infinita) della distribuzione, dove x > x$_1$ (e y < y$_1$). Questo livello è detto strato 0 ed ha area pari ad A.
Sopra di esso, viene costruito un altro strato rettangolare di larghezza x$_1$ e altezza A / x$_1$, così da mantenere anche per esso un'area pari ad A. La parte superiore di questo livello si trova a y$_2$ = y$_1$ + A / x$_1$ e interseca la funzione di densità di probabilità in un punto (x$_2$, y$_2$), dove y $_2$ = f (x$_2$). Questo strato include ogni punto della funzione di densità di probabilità compreso tra y$_1$ e y$_2$, ma, a differenza del livello base, include anche punti come (x$_1$, y$_2$) che non si trovano nella distribuzione desiderata.
Ulteriori strati vengono poi impilati mantenendo costante l'area. Per utilizzare una tabella precalcolata di dimensione n (256 solitamente), viene scelta x$_1$ tale che x$_n$ = 0, il che significa che il livello superiore, il livello n - 1, raggiunge esattamente il picco di distribuzione a (0, f (0)).
Generalizzando, lo strato i-esimo si estende verticalmente da y$_i$ a y$_{i+1}$ e può essere diviso orizzontalmente in due regioni: la porzione (generalmente più grande) da 0 a x$_{i+1}$ che è interamente contenuta nella distribuzione desiderata, e la porzione (generalmente più piccola) da x$_{i+1}$ a x$_i$, che è solo parzialmente contenuta nella distribuzione.
Ignorando per un momento il problema del livello 0, e date le variabili casuali uniformi U$_0$ e U$_1$ ∈ [0,1], l'algoritmo ziggurat opera in questi passi:
* Scegli un livello casuale 0 ≤ i < n .
* Sia x = U$_0$x$_i$ .
* Se x < x$_{i+1}$, restituisci x.
* Sia y = y$_i$ + U$_1$(y$_{i+1}$ - y$_i$ ).
* Calcola f (x). Se y < f (x), restituisci x.
* Altrimenti, scegli nuovi numeri casuali e torna al passaggio 1.
Il passaggio 1 equivale alla scelta di una coordinata y. Il passaggio 3 verifica se la coordinata x è all'interno della funzione di densità desiderata. In caso contrario, il passo 4 sceglie una coordinata y e il passo 5 esegue il test di accettazione.
Si noti che per il livello superiore n - 1 questo test fallisce sempre, perché x $_n$ = 0.

### Ottimizzazioni

L'algoritmo può essere eseguito in modo efficiente utilizzando tabelle precalcolate per x$_i$ e y$_i$ = f (x$_i$), ma si possono apportare alcuni miglioramenti per renderlo ancora più veloce:
* L'algoritmo ziggurat non ha niente che dipenda dalla normalizzazione della funzione di distribuzione di probabilità (integrale sotto la curva uguale a 1), quindi la rimozione delle costanti di normalizzazione può accelerare il calcolo di f (x).
* Invece di confrontare U$_0$x$_i$ con x$_{i+1}$ al punto 3, è possibile precalcolare x$_{i+1}$ / x$_i$ e confrontarlo con U$_0$ direttamente. Se U$_0$ è un generatore di numeri casuali intero, questi limiti possono essere moltiplicati per 2$^{32}$ (o 2$^{31}$, a seconda dei casi) in modo da poter utilizzare un confronto tra numeri interi.
* Quando si generano valori in virgola mobile a precisione singola (a standard IEEE 754), che hanno solo una mantissa a 24 bit (incluso l'iniziale implicita 1), i bit meno significativi di un numero casuale intero a 32 bit non vengono utilizzati. Questi bit possono essere usati per selezionare il numero di livelli.