# Algoritmo zigurate

O algoritmo zigurate é um algoritmo para amostragem de números pseudoaleatórios . Pertencente à classe de algoritmos de amostragem de rejeição, ele depende de uma fonte subjacente de números aleatórios uniformemente distribuídos, normalmente de um gerador de números pseudoaleatórios, bem como de tabelas pré-calculadas. O algoritmo é usado para gerar valores de uma distribuição de probabilidade monotonicamente decrescente . Também pode ser aplicado a distribuições unimodais simétricas, como a distribuição normal, escolhendo um valor de uma metade da distribuição e, em seguida, escolhendo aleatoriamente de qual metade o valor é considerado como tendo sido extraído. Foi desenvolvido por George Marsaglia e outros na década de 1960.
Um valor típico produzido pelo algoritmo requer apenas a geração de um valor de ponto flutuante aleatório e um índice de tabela aleatório, seguido por uma consulta de tabela, uma operação de multiplicação e uma comparação. Às vezes (2,5% do tempo, no caso de uma distribuição normal ou exponencial ao usar tamanhos de tabela típicos) mais cálculos são necessários. No entanto, o algoritmo é computacionalmente muito mais rápidodo que os dois métodos mais comumente usados para gerar números aleatórios distribuídos normalmente, o método polar de Marsaglia e a transformada de Box-Muller, que requerem pelo menos um logaritmo e um cálculo de raiz quadrada para cada par de valores gerados. No entanto, como o algoritmo zigurate apresenta maior complexidade de implementação, ele é melhor usado quando grandes quantidades de números aleatórios são necessárias.
O termo algoritmo zigurate data do artigo de Marsaglia com Wai Wan Tsang em 2000; é assim chamado porque é conceitualmente baseado na cobertura da distribuição de probabilidade com segmentos retangulares empilhados em ordem decrescente de tamanho, resultando em uma figura que se assemelha a um zigurate .

## Teoria da operação

O algoritmo zigurate é um algoritmo de amostragem de rejeição; ele gera aleatoriamente um ponto em uma distribuição ligeiramente maior que a distribuição desejada e, então, testa se o ponto gerado está dentro da distribuição desejada. Caso contrário, ele tenta novamente. Dado um ponto aleatório abaixo de uma curva de densidade de probabilidade, sua coordenada x é um número aleatório com a distribuição desejada.
A distribuição escolhida pelo algoritmo zigurate é composta por n regiões de área igual; n − 1 retângulos que cobrem a maior parte da distribuição desejada, sobre uma base não retangular que inclui a cauda da distribuição.
Dada uma função de densidade de probabilidade decrescente monótona f(x), definida para todo x ≥ 0, a base do zigurate é definida como todos os pontos dentro da distribuição e abaixo de y $_1$ = f(x$_1$). Consiste em uma região retangular de (0, 0) a (x$_1$ , y$_1$), e a cauda (tipicamente infinita) da distribuição, onde x > x$_1$ (e y < y$_1$ ).
Esta camada (vamos chamá-la camada 0) tem área A . Em cima dela, adicione uma camada retangular de largura x$_1$ e altura A / x$_1$, para que ela também tenha área A. O topo desta camada está na altura y $_2$ = y$_1$ + A / x$_1$, e intercepta a função de densidade em um ponto ( x$_2$ , y$_2$ ), onde y$_2$ = f(x$_2$). Esta camada inclui todos os pontos na função de densidade entre y $_1$ e y $_2$, mas (ao contrário da camada base) também inclui pontos como (x$_1$ , y$_2$) que não estão na distribuição desejada.
Outras camadas são então empilhadas em cima. Para usar uma tabela pré-computada de tamanho n ( n = 256 é típico), escolhe-se x$_1$ tal que x$_n$ = 0, o que significa que a caixa superior, camada n − 1, atinge o pico da distribuição em (0, f(0)) exatamente.
A camada i se estende verticalmente de y $_i$ para y$_{i+1}$, e pode ser dividido em duas regiões horizontalmente: a porção (geralmente maior) de 0 a x$_{i+1}$ que está inteiramente contido na distribuição desejada, e a (pequena) porção de x$_{i+1}$ para x$_i$, que está apenas parcialmente contido.
Ignorando por um momento o problema da camada 0, e dadas as variáveis aleatórias uniformes U$_0$ e U$_1$ ∈ [0,1), o algoritmo do zigurate pode ser descrito como:
* Escolha uma camada aleatória 0 ≤ i < n .
* Seja x = U$_0$x$_i$ .
* Se x < x$_{i+1}$, retorna x .
* Seja y = y$_i$ + U$_1$ ( y$_{i+1}$ − y$_i$ ).
* Calcule f(x). Se y < f(x), retorne x .
* Caso contrário, escolha novos números aleatórios e volte para a etapa 1.
O passo 1 consiste em escolher uma coordenada y de baixa resolução. A etapa 3 testa se a coordenada x está claramente dentro da função de densidade desejada sem saber mais sobre a coordenada y. Caso contrário, a etapa 4 escolhe uma coordenada y de alta resolução e a etapa 5 faz o teste de rejeição.
Com camadas muito próximas, o algoritmo termina na etapa 3 em uma fração muito grande do tempo. Para a camada superior n − 1, no entanto, esse teste sempre falha, porque x$_n$ = 0.
A camada 0 também pode ser dividida em uma região central e uma borda, mas a borda é uma cauda infinita. Para usar o mesmo algoritmo para verificar se o ponto está na região central, gere um x$_0$ = A / y$_1$ fictício. Isso gerará pontos com x < x$_1$ com a frequência correta e, no caso raro de a camada 0 ser selecionada e x ≥ x$_1$, use um algoritmo de fallback especial para selecionar um ponto aleatoriamente na cauda. Como o algoritmo de fallback é usado menos de uma vez em mil, a velocidade não é essencial.
Assim, o algoritmo zigurate completo para distribuições unilaterais é:
* Escolha uma camada aleatória 0 ≤ i < n .
* Seja x = U$_0$x$_i$
* Se x < x$_{i+1}$, retorna x .
* Se i = 0, gere um ponto a partir da cauda usando o algoritmo de fallback.
* Seja y = y$_i$ + U$_1$(y$_{i+1}$ − y$_i$ ).
* Calcule f (x). Se y < f (x), retorne x .
* Caso contrário, escolha novos números aleatórios e volte para a etapa 1.
Para uma distribuição bilateral, o resultado deve ser negado 50% das vezes. Isso pode ser feito convenientemente escolhendo U$_0$ ∈ (−1,1) e, na etapa 3, testando se | x | < x$_{i+1}$ .

### Algoritmos de fallback para a cauda

Porque o algoritmo zigurate só gera a maioria das saídas muito rapidamente e requer um algoritmo de fallback sempre que x > x$_1$, é sempre mais complexo do que uma implementação mais direta. O algoritmo de fallback específico depende da distribuição.
Para uma distribuição exponencial, a cauda se parece com o corpo da distribuição. Uma maneira é recorrer ao algoritmo mais elementar E = −ln( U$_1$ ) e deixe x = x$_1$ − em( U$_1$ ). Outra é chamar o algoritmo do zigurate recursivamente e adicionar x$_1$ ao resultado.
Para uma distribuição normal, Marsaglia sugere um algoritmo compacto:
* Seja x = −ln(U$_1$)/ x$_1$ .
* Seja y = −ln(U$_2$).
* Se 2y > x$^2$, retorne x + x$_1$ .
* Caso contrário, volte para a etapa 1.
Como x $_1$ ≈ 3,5 para tamanhos de tabela típicos, o teste na etapa 3 quase sempre é bem-sucedido. Como −ln(U$_1$) é uma variável distribuída exponencialmente, uma implementação da distribuição exponencial pode ser usada.

### Otimizações

O algoritmo pode ser executado eficientemente com tabelas pré-calculadas de x $_i$ e y $_i$ = f (x$_i$), mas há algumas modificações para torná-lo ainda mais rápido:
* Nada no algoritmo do zigurate depende da função de distribuição de probabilidade ser normalizada (integral sob a curva igual a 1). A remoção de constantes de normalização pode acelerar o cálculo de f(x).
* A maioria dos geradores uniformes de números aleatórios são baseados em geradores de números aleatórios inteiros que retornam um inteiro no intervalo [0, 2$^{32}$ − 1]. Uma tabela de 2$^{−32}$ x $_i$ permite que você use esses números diretamente para U$_0$ .
* Ao calcular distribuições de dois lados usando um U$_0$ de dois lados, conforme descrito anteriormente, o inteiro aleatório pode ser interpretado como um número assinado no intervalo [−2$^{31}$ , 2$^{31}$ − 1], e um fator de escala de 2$^{−31}$ pode ser usado.
* Em vez de comparar U$_0$x$_i$ com x$_{i+1}$ na etapa 3, é possível pré-calcular x$_{i+1}$/x $_i$ e compare U$_0$ com isso diretamente. Se U$_0$ for um gerador de números aleatórios inteiros, esses limites podem ser pré-multiplicados por 2$^{32}$ (ou 2$^{31}$, conforme apropriado) para que uma comparação inteira possa ser usada.
* Com as duas alterações acima, a tabela de valores x$_i$ não modificados não é mais necessária e pode ser excluída.
* Ao gerar valores de ponto flutuante de precisão simples IEEE 754, que têm apenas uma mantissa de 24 bits (incluindo o 1 inicial implícito), os bits menos significativos de um número aleatório inteiro de 32 bits não são usados. Esses bits podem ser usados para selecionar o número da camada. (Veja as referências abaixo para uma discussão detalhada sobre isso.)
* Os três primeiros passos podem ser colocados em uma função inline, que pode chamar uma implementação fora de linha dos passos menos frequentemente necessários.

### Gerando as tabelas

É possível armazenar a tabela inteira pré-calculada ou apenas incluir os valores n, y$_1$, A e uma implementação de f$^{−1}$(y) no código-fonte e calcule os valores restantes ao inicializar o gerador de números aleatórios.
Conforme descrito anteriormente, você pode encontrar x$_i$ = f$^{−1}$(y$_i$) e y$_{i+1}$ = y$_i$ + A/x$_i$ . Repita n − 1 vez para cada um das camadas do zigurate. No final, você deve ter y$_n$ = f(0). Haverá algum erro de arredondamento, mas é um teste de sanidade útil para verificar se ele é aceitavelmente pequeno.
Ao preencher os valores da tabela, basta assumir que x$_n$ = 0 e y$_n$ = f(0), e aceitar a ligeira diferença na área da camada n − 1 como erro de arredondamento.

### Encontrando x$_1$ e A

Dado um valor inicial (estimado) x$_1$, você precisa de uma maneira de calcular a área t da cauda para a qual x > x$_1$. Para a distribuição exponencial, isto é apenas e$^{−x1}$, enquanto para a distribuição normal, supondo que você esteja usando o não normalizado f (x) = e$^{−x2/2}$, isto é √π/2 erfc (x /√2). Para distribuições mais complicadas, pode ser necessária integração numérica.
Com isso em mãos, a partir de x $_1$, você pode encontrar y $_1$ = f (x$_1$), a área t na cauda e a área da camada base A = x$_1$y$_1$ + t .
Em seguida, calcule as séries y$_i$ e x$_i$ como acima. Se y$_i$ > f(0) para qualquer i < n, então a estimativa inicial x$_1$ era muito baixa, levando a uma área A muito grande. Se y$_n$ < f (0), então a estimativa inicial x$_1$ era muito alta.
Dado isso, use um algoritmo de busca de raízes (como o método da bissecção) para encontrar o valor x$_1$ que produz y$_{n−1}$ o mais próximo possível de f (0). Alternativamente, procure o valor que torna a área da camada superior, x$_{n−1}$(f(0) − y$_{n−1}$), o mais próximo possível do valor desejado A. Isso economiza uma avaliação de f$^{−1}$( x ) e é na verdade a condição de maior interesse.

## Variação de McFarland

Christopher D. McFarland propôs uma versão ainda mais otimizada. Isso aplica três alterações algorítmicas, às custas de tabelas um pouco maiores.
Primeiro, o caso comum considera apenas as porções retangulares, de (0, y$_{i−1}$ ) para (x$_i$ , y$_i$) As regiões de formato estranho à direita destas (a maioria quase triangular, mais a cauda) são tratadas separadamente. Isso simplifica e acelera o caminho rápido do algoritmo.
Em segundo lugar, a área exata das regiões de formato ímpar é usada; elas não são arredondadas para incluir o retângulo inteiro para (x$_{i-1}$ , y$_i$ ). Isso aumenta a probabilidade de que o caminho rápido seja usado.
Uma consequência importante disso é que o número de camadas é ligeiramente menor que n . Mesmo que a área das porções de formato ímpar seja tomada exatamente, o total soma mais do que o valor de uma camada. A área por camada é ajustada para que o número de camadas retangulares seja um número inteiro. Se o 0 inicial ≤ i < n excede o número de camadas retangulares, a fase 2 prossegue.
Se o valor procurado estiver em qualquer uma das regiões de formato ímpar, o método de alias será usado para escolher uma, com base em sua área real. Isso representa uma pequena quantidade de trabalho adicional e requer tabelas de alias adicionais, mas escolhe um dos lados direitos das camadas.
A região de formato ímpar escolhida é rejeitada, mas se uma amostra for rejeitada, o algoritmo não retorna ao início. A área real de cada região de formato estranho foi usada para escolher uma camada, de modo que o loop de rejeição-amostragem permanece nessa camada até que um ponto seja escolhido.
Terceiro, o formato quase triangular da maioria das porções de formatos estranhos é explorado, embora isso deva ser dividido em três casos, dependendo da segunda derivada da função de distribuição de probabilidade na camada selecionada.
Se a função for convexa (como a distribuição exponencial está em toda parte, e a distribuição normal é para |x| > 1), então a função está estritamente contida dentro do triângulo inferior. Dois desvios uniformes unitários U$_1$ e U$_2$ são escolhidos e, antes de serem dimensionados para o retângulo que envolve a região de formato ímpar, sua soma é testada. Se U$_1$ + U$_2$ > 1, o ponto está no triângulo superior e pode ser refletido para (1− U$_1$ , 1− U$_2$). Então, se U$_1$ + U$_2$ < 1− ε, para alguma tolerância adequada ε, o ponto está definitivamente abaixo da curva e pode ser imediatamente aceito. Somente para pontos muito próximos da diagonal é necessário calcular a função de distribuição f(x) para executar um teste de rejeição exato. (A tolerância ε deve, em teoria, depender da camada, mas um único valor máximo pode ser usado em todas as camadas com pouca perda.)
Se a função for côncava (como a distribuição normal é para |x| < 1), inclui uma pequena porção do triângulo superior, de modo que a reflexão é impossível, mas pontos cujas coordenadas normalizadas satisfazem U$_1$ + U$_2$ ≤ 1 podem ser imediatamente aceitos, e pontos para os quais U$_1$ + U$_2$ > 1+ε podem ser imediatamente rejeitados.
Na camada que atravessa |x| = 1, a distribuição normal tem um ponto de inflexão e o teste de rejeição exato deve ser aplicado se 1− ε < U$_1$ + U$_2$ < 1+ ε.
A cauda é tratada como no algoritmo zigurate original e pode ser considerada um quarto caso para o formato da região de formato ímpar à direita.