Algoritmo ziggurat
# 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_0x_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_0x_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.
