El algoritmo ziggurat es un algoritmo de muestreo de números pseudoaleatorios . Perteneciente a la clase de algoritmos de muestreo de rechazo , se basa en una fuente subyacente de números aleatorios distribuidos uniformemente, normalmente de un generador de números pseudoaleatorios , así como de tablas precalculadas. El algoritmo se utiliza para generar valores a partir de una distribución de probabilidad monótonamente decreciente . También se puede aplicar a distribuciones unimodales simétricas , como la distribución normal , eligiendo un valor de la mitad de la distribución y luego eligiendo aleatoriamente de qué mitad se considera que se ha extraído el valor. Fue desarrollado por George Marsaglia y otros en la década de 1960.
Un valor típico producido por el algoritmo solo requiere la generación de un valor de punto flotante aleatorio y un índice de tabla aleatorio, seguido de una búsqueda en la tabla, una operación de multiplicación y una comparación. A veces (el 2,5 % del tiempo, en el caso de una distribución normal o exponencial cuando se utilizan tamaños de tabla típicos) [ cita requerida ] se requieren más cálculos. Sin embargo, el algoritmo es computacionalmente mucho más rápido [ cita requerida ] que los dos métodos más comúnmente utilizados para generar números aleatorios distribuidos normalmente, el método polar de Marsaglia y la transformada de Box-Muller , que requieren al menos un cálculo de logaritmo y una raíz cuadrada para cada par de valores generados. Sin embargo, dado que el algoritmo ziggurat es más complejo de implementar, se utiliza mejor cuando se requieren grandes cantidades de números aleatorios.
El término algoritmo zigurat data del artículo de Marsaglia con Wai Wan Tsang en 2000; se llama así porque se basa conceptualmente en cubrir la distribución de probabilidad con segmentos rectangulares apilados en orden decreciente de tamaño, dando como resultado una figura que se asemeja a un zigurat .
El algoritmo ziggurat es un algoritmo de muestreo por rechazo; genera aleatoriamente un punto en una distribución ligeramente mayor que la distribución deseada y luego prueba si el punto generado está dentro de la distribución deseada. Si no es así, vuelve a intentarlo. Dado un punto aleatorio debajo de una curva de densidad de probabilidad, su coordenada x es un número aleatorio con la distribución deseada.
La distribución que elige el algoritmo zigurat se compone de n regiones de área igual; n − 1 rectángulos que cubren la mayor parte de la distribución deseada, sobre una base no rectangular que incluye la cola de la distribución.
Dada una función de densidad de probabilidad decreciente y monótona f ( x ), definida para todo x ≥ 0, la base del zigurat se define como todos los puntos dentro de la distribución y por debajo de y 1 = f ( x 1 ). Esta consiste en una región rectangular desde (0, 0) hasta ( x 1 , y 1 ), y la cola (normalmente infinita) de la distribución, donde x > x 1 (e y < y 1 ).
Esta capa (llamémosla capa 0) tiene un área A. Sobre esta, añada una capa rectangular de ancho x 1 y altura A / x 1 , por lo que también tiene un área A. La parte superior de esta capa está a una altura y 2 = y 1 + A / x 1 , e interseca la función de densidad en un punto ( x 2 , y 2 ), donde y 2 = f ( x 2 ). Esta capa incluye todos los puntos de la función de densidad entre y 1 e y 2 , pero (a diferencia de la capa base) también incluye puntos como ( x 1 , y 2 ) que no están en la distribución deseada.
Luego se apilan más capas sobre la parte superior. Para utilizar una tabla precalculada de tamaño n ( n = 256 es lo típico), se elige x 1 de modo que x n = 0, lo que significa que el cuadro superior, capa n − 1, alcanza el pico de la distribución exactamente en (0, f (0)).
La capa i se extiende verticalmente desde y i hasta y i +1 , y se puede dividir en dos regiones horizontalmente: la porción (generalmente más grande) de 0 a x i +1 que está completamente contenida dentro de la distribución deseada, y la porción (pequeña) de x i +1 a x i , que está solo parcialmente contenida.
Ignorando por un momento el problema de la capa 0, y dadas las variables aleatorias uniformes U 0 y U 1 ∈ [0,1), el algoritmo del zigurat se puede describir como:
El paso 1 consiste en elegir una coordenada y de baja resolución . El paso 3 comprueba si la coordenada x está claramente dentro de la función de densidad deseada sin saber más sobre la coordenada y. Si no lo está, el paso 4 elige una coordenada y de alta resolución y el paso 5 realiza la prueba de rechazo.
Con capas muy espaciadas, el algoritmo termina en el paso 3 una fracción muy grande de las veces. Sin embargo, para la capa superior n − 1, esta prueba siempre falla, porque x n = 0.
La capa 0 también se puede dividir en una región central y un borde, pero el borde es una cola infinita. Para utilizar el mismo algoritmo para comprobar si el punto está en la región central, genere un x 0 = A / y 1 ficticio . Esto generará puntos con x < x 1 con la frecuencia correcta y, en el caso poco frecuente de que se seleccione la capa 0 y x ≥ x 1 , utilice un algoritmo de reserva especial para seleccionar un punto al azar de la cola. Debido a que el algoritmo de reserva se utiliza menos de una vez en mil, la velocidad no es esencial.
Por lo tanto, el algoritmo zigurat completo para distribuciones unilaterales es:
En el caso de una distribución bilateral, el resultado debe ser negado el 50 % de las veces. Esto se puede hacer de manera conveniente eligiendo U 0 ∈ (−1,1) y, en el paso 3, comprobando si | x | < x i +1 .
Dado que el algoritmo zigurat solo genera la mayoría de los resultados muy rápidamente y requiere un algoritmo de respaldo siempre que x > x 1 , siempre es más complejo que una implementación más directa. El algoritmo de respaldo específico depende de la distribución.
En una distribución exponencial, la cola se parece al cuerpo de la distribución. Una forma de hacerlo es recurrir al algoritmo más elemental E = −ln( U 1 ) y dejar x = x 1 − ln( U 1 ). Otra forma de hacerlo es llamar al algoritmo del zigurat de forma recursiva y sumar x 1 al resultado.
Para una distribución normal, Marsaglia sugiere un algoritmo compacto:
Como x 1 ≈ 3,5 para tamaños de tabla típicos, la prueba del paso 3 casi siempre es exitosa. Como −ln( U 1 ) es una variable distribuida exponencialmente, se puede utilizar una implementación de la distribución exponencial.
El algoritmo se puede ejecutar de manera eficiente con tablas precalculadas de x i e y i = f ( x i ) , pero hay algunas modificaciones para hacerlo aún más rápido:
Es posible almacenar la tabla completa precalculada, o solo incluir los valores n , y 1 , A y una implementación de f −1 ( y ) en el código fuente, y calcular los valores restantes al inicializar el generador de números aleatorios.
Como se describió anteriormente, puedes encontrar x i = f −1 ( y i ) y y i +1 = y i + A / x i . Repite n − 1 veces para las capas del zigurat. Al final, deberías tener y n = f (0). Habrá algún error de redondeo , pero es una prueba de cordura útil para ver que es aceptablemente pequeño.
Al completar los valores de la tabla, simplemente suponga que x n = 0 e y n = f (0), y acepte la ligera diferencia en el área de la capa n − 1 como error de redondeo.
Dado un valor inicial (aproximado) x 1 , se necesita una forma de calcular el área t de la cola para la cual x > x 1 . Para la distribución exponencial, esto es simplemente e − x 1 , mientras que para la distribución normal, suponiendo que se está utilizando la f ( x ) = e − x 2 /2 no normalizada , esto es √ π /2 erfc ( x / √ 2 ). Para distribuciones más complicadas, puede ser necesaria la integración numérica .
Con esto en la mano, a partir de x 1 , puedes encontrar y 1 = f ( x 1 ), el área t en la cola y el área de la capa base A = x 1 y 1 + t .
Luego, calcule las series y i y x i como se indicó anteriormente. Si y i > f (0) para cualquier i < n , entonces la estimación inicial x 1 fue demasiado baja, lo que llevó a un área A demasiado grande . Si y n < f (0), entonces la estimación inicial x 1 fue demasiado alta.
Teniendo esto en cuenta, utilice un algoritmo de búsqueda de raíces (como el método de bisección ) para encontrar el valor x 1 que produce y n −1 lo más cerca posible de f (0). Alternativamente, busque el valor que hace que el área de la capa superior, x n −1 ( f (0) − y n −1 ), sea lo más cercana posible al valor deseado A. Esto ahorra una evaluación de f −1 ( x ) y es, en realidad, la condición de mayor interés.
Christopher D. McFarland propuso una versión aún más optimizada. [1] Esta aplica tres cambios algorítmicos, a expensas de tablas ligeramente más grandes.
En primer lugar, el caso común considera solo las porciones rectangulares, desde (0, y i −1 ) hasta ( x i , y i ). Las regiones de forma irregular a la derecha de estas (en su mayoría casi triangulares, más la cola) se manejan por separado. Esto simplifica y acelera la ruta rápida del algoritmo .
En segundo lugar, se utiliza el área exacta de las regiones de forma irregular; no se redondean para incluir todo el rectángulo a ( x i −1 , y i ). Esto aumenta la probabilidad de que se utilice la ruta rápida.
Una consecuencia importante de esto es que el número de capas es ligeramente menor que n . Aunque el área de las partes con formas irregulares se toma exactamente, el total suma más del valor de una capa. El área por capa se ajusta de modo que el número de capas rectangulares sea un número entero. Si el valor inicial 0 ≤ i < n excede el número de capas rectangulares, se continúa con la fase 2.
Si el valor buscado se encuentra en alguna de las regiones de forma irregular, se utiliza el método de alias para elegir una, en función de su área real. Esto supone una pequeña cantidad de trabajo adicional y requiere tablas de alias adicionales, pero elige uno de los lados derechos de las capas.
La región de forma irregular elegida se somete a un muestreo de rechazo, pero si se rechaza una muestra, el algoritmo no vuelve al principio. Se utilizó el área real de cada región de forma irregular para elegir una capa, por lo que el bucle de muestreo de rechazo permanece en esa capa hasta que se elige un punto.
En tercer lugar, se aprovecha la forma casi triangular de la mayoría de las porciones de forma irregular, aunque esto debe dividirse en tres casos dependiendo de la segunda derivada de la función de distribución de probabilidad en la capa seleccionada.
Si la función es convexa (ya que la distribución exponencial está en todas partes y la distribución normal es para | x | > 1), entonces la función está estrictamente contenida dentro del triángulo inferior. Se eligen dos desviaciones uniformes unitarias U 1 y U 2 y, antes de escalarlas al rectángulo que encierra la región de forma irregular, se prueba su suma. Si U 1 + U 2 > 1, el punto está en el triángulo superior y se puede reflejar en (1− U 1 , 1− U 2 ). Entonces, si U 1 + U 2 < 1− ε , para alguna tolerancia adecuada ε , el punto está definitivamente debajo de la curva y se puede aceptar de inmediato. Solo para puntos muy cercanos a la diagonal es necesario calcular la función de distribución f ( x ) para realizar una prueba de rechazo exacta. (La tolerancia ε debería depender en teoría de la capa, pero se puede usar un único valor máximo en todas las capas con poca pérdida).
Si la función es cóncava (como lo es la distribución normal para | x | < 1), incluye una pequeña porción del triángulo superior, por lo que la reflexión es imposible, pero los puntos cuyas coordenadas normalizadas satisfacen U 1 + U 2 ≤ 1 pueden aceptarse inmediatamente, y los puntos para los que U 1 + U 2 > 1+ ε pueden rechazarse inmediatamente.
En la capa que se extiende a ambos lados de | x | = 1, la distribución normal tiene un punto de inflexión y se debe aplicar la prueba de rechazo exacta si 1− ε < U 1 + U 2 < 1+ ε .
La cola se maneja como en el algoritmo Ziggurat original y puede considerarse como un cuarto caso para la forma de la región de forma irregular a la derecha.
[C]uando la prioridad número uno es mantener una calidad estadística extremadamente alta y, sujeto a esa restricción, también se desea velocidad, el método Ziggurat será a menudo la opción más adecuada.Comparación de varios algoritmos para generar números aleatorios gaussianos .