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 .

O algoritmo zigurate é usado para gerar valores de amostra com uma distribuição normal . (Somente valores positivos são mostrados para simplificar.) Os pontos rosas são inicialmente números aleatórios uniformemente distribuídos. A função de distribuição desejada é primeiro segmentada em áreas iguais "A". Uma camada i é selecionada aleatoriamente pela fonte uniforme à esquerda. Em seguida, um valor aleatório da fonte superior é multiplicado pela largura da camada escolhida, e o resultado é testado em x para ver em qual região da camada ele se enquadra, com 3 resultados possíveis: 1) (esquerda, região preta sólida) a amostra está claramente abaixo da curva e pode ser emitida imediatamente, 2) (direita, região listrada verticalmente) o valor da amostra pode estar abaixo da curva e deve ser testado mais detalhadamente. Nesse caso, um valor y aleatório dentro da camada escolhida é gerado e comparado com f(x) . Se for menor, o ponto estará abaixo da curva e o valor x será gerado. Caso contrário (terceiro caso), o ponto escolhido x é rejeitado e o algoritmo é reiniciado do início.

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 ₁ = f(x₁). Consiste em uma região retangular de (0, 0) a (x₁ , y₁), e a cauda (tipicamente infinita) da distribuição, onde x > x₁ (e y < y₁ ).

Esta camada (vamos chamá-la camada 0) tem área A . Em cima dela, adicione uma camada retangular de largura x₁ e altura A / x₁, para que ela também tenha área A. O topo desta camada está na altura y ₂ = y₁ + A / x₁, e intercepta a função de densidade em um ponto ( x₂ , y₂ ), onde y₂ = f(x₂). Esta camada inclui todos os pontos na função de densidade entre y ₁ e y ₂, mas (ao contrário da camada base) também inclui pontos como (x₁ , y₂) 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₁ 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₀ e U₁ ∈ [0,1), o algoritmo do zigurate pode ser descrito como:

1. Escolha uma camada aleatória 0 ≤ i < n .
2. Seja x = U₀x_i .
3. Se x < x_i +1, retorna x .
4. Seja y = y_i + U₁ ( y_i +1 − y_i ).
5. Calcule f(x). Se y < f(x), retorne x .
6. 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₀ = A / y₁ fictício. Isso gerará pontos com x < x₁ com a frequência correta e, no caso raro de a camada 0 ser selecionada e x ≥ x₁, 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 é:

1. Escolha uma camada aleatória 0 ≤ i < n .
2. Seja x = U₀x_i
3. Se x < x_i +1, retorna x .
4. Se i = 0, gere um ponto a partir da cauda usando o algoritmo de fallback.
5. Seja y = y_i + U₁(y_i +1 − y_i ).
6. Calcule f (x). Se y < f (x), retorne x .
7. 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₀ ∈ (−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₁, é 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₁ ) e deixe x = x₁ − em( U₁ ). Outra é chamar o algoritmo do zigurate recursivamente e adicionar x₁ ao resultado.

Para uma distribuição normal, Marsaglia sugere um algoritmo compacto:

1. Seja x = −ln(U₁)/ x₁ .
2. Seja y = −ln(U₂).
3. Se 2y > x², retorne x + x₁ .
4. Caso contrário, volte para a etapa 1.

Como x ₁ ≈ 3,5 para tamanhos de tabela típicos, o teste na etapa 3 quase sempre é bem-sucedido. Como −ln(U₁) é 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³² − 1]. Uma tabela de 2⁻³² x _i permite que você use esses números diretamente para U₀ .
• Ao calcular distribuições de dois lados usando um U₀ de dois lados, conforme descrito anteriormente, o inteiro aleatório pode ser interpretado como um número assinado no intervalo [−2³¹ , 2³¹ − 1], e um fator de escala de 2⁻³¹ pode ser usado.
• Em vez de comparar U₀x_i com x_i +1 na etapa 3, é possível pré-calcular x_i +1/x _i e compare U₀ com isso diretamente. Se U₀ for um gerador de números aleatórios inteiros, esses limites podem ser pré-multiplicados por 2³² (ou 2³¹, 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₁, A e uma implementação de f⁻¹(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⁻¹(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₁ e A

Dado um valor inicial (estimado) x₁, você precisa de uma maneira de calcular a área t da cauda para a qual x > x₁. Para a distribuição exponencial, isto é apenas e^{− x₁}, 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 ₁, você pode encontrar y ₁ = f (x₁), a área t na cauda e a área da camada base A = x₁y₁ + 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₁ era muito baixa, levando a uma área A muito grande. Se y_n < f (0), então a estimativa inicial x₁ 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₁ 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⁻¹( 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₁ e U₂ 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₁ + U₂ > 1, o ponto está no triângulo superior e pode ser refletido para (1− U₁ , 1− U₂). Então, se U₁ + U₂ < 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₁ + U₂ ≤ 1 podem ser imediatamente aceitos, e pontos para os quais U₁ + U₂ > 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₁ + U₂ < 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.