L' algorithme de la ziggourat est un algorithme d' échantillonnage de nombres pseudo-aléatoires . Appartenant à la classe des algorithmes d'échantillonnage par rejet , il repose...
La génération d'une valeur typique produite par l'algorithme ne nécessite que la génération d'une valeur aléatoire à virgule flottante et d'un index de table aléatoire, suivies d'une consultation de table, d'une multiplication et d'une comparaison. Dans 2,5 % des cas (pour une distribution normale ou exponentielle avec des tables de taille standard), des calculs supplémentaires sont parfois nécessaires. Néanmoins, l'algorithme est beaucoup plus rapide que les deux méthodes les plus couramment utilisées pour générer des nombres aléatoires suivant une loi normale : la méthode polaire de Marsaglia et la transformée de Box-Muller . Ces dernières requièrent au moins un calcul de logarithme et un calcul de racine carrée pour chaque paire de valeurs générées. Cependant, l'algorithme de ziggourat étant plus complexe à implémenter, il est préférable de l'utiliser lorsque de grandes quantités de nombres aléatoires sont requises.
Le terme « algorithme de ziggourat » date de l'article de Marsaglia avec Wai Wan Tsang en 2000 ; il est ainsi nommé car il est conceptuellement basé sur la couverture de la distribution de probabilité avec des segments rectangulaires empilés par ordre décroissant de taille, ce qui donne une figure qui ressemble à une ziggourat .
L'algorithme Ziggurat est utilisé pour générer des valeurs d'échantillon suivant une distribution normale . (Seules les valeurs positives sont affichées par souci de simplicité.) Les points roses représentent initialement des nombres aléatoires uniformément distribués. La fonction de distribution souhaitée est d'abord segmentée en zones égales « A ». Une couche i est sélectionnée aléatoirement par la source uniforme située à gauche. Ensuite, une valeur aléatoire issue de la source supérieure est multipliée par la largeur de la couche choisie, et le résultat x est testé afin de déterminer dans quelle région de la couche il se situe, avec trois résultats possibles : 1) (à gauche, zone noire continue) l'échantillon se situe clairement sous la courbe et peut être affiché immédiatement ; 2) (à droite, zone à rayures verticales) la valeur de l'échantillon peut se situer sous la courbe et doit être testée plus en détail. Dans ce cas, une valeur aléatoire y est générée dans la couche choisie et comparée à f(x) . Si elle est inférieure, le point se situe sous la courbe et la valeur x est affichée. Sinon (troisième cas), le point x choisi est rejeté et l'algorithme est redémarré.
Théorie de fonctionnement
L'algorithme de ziggourat est un algorithme d'échantillonnage par rejet ; il génère aléatoirement un point dans une distribution légèrement plus large que la distribution cible, puis vérifie si ce point appartient à la distribution cible. Si ce n'est pas le cas, il recommence. Étant donné un point aléatoire situé sous une courbe de densité de probabilité, son abscisse est un nombre aléatoire appartenant à la distribution cible.
La distribution choisie par l'algorithme de ziggourat est composée de n régions de surface égale ;
Étant donné une fonction de densité de probabilité monotone décroissante
Cette couche (appelons-la couche 0) a une aire A. Par-dessus, ajoutez une couche rectangulaire de largeur
D'autres couches sont ensuite empilées par-dessus. Pour utiliser un tableau précalculé de taille
Couche
En négligeant un instant le problème de la couche 0, et en considérant des variables aléatoires uniformes
Choisissez un calque aléatoire
Laisser
Si
Laisser
Calculer
Sinon, choisissez de nouveaux nombres aléatoires et retournez à l'étape 1.
L'étape 1 consiste à choisir une coordonnée y de faible résolution . L'étape 3 vérifie si la coordonnée x se situe clairement dans la fonction de densité souhaitée, sans connaître davantage la coordonnée y. Si ce n'est pas le cas, l'étape 4 choisit une coordonnée y de haute résolution, et l'étape 5 effectue le test de rejet.
Avec des couches très rapprochées, l'algorithme s'arrête à l'étape 3 dans une très grande partie des cas. Pour la couche supérieure
La couche 0 peut également être divisée en une région centrale et un bord, mais ce bord est une queue infinie. Pour utiliser le même algorithme afin de vérifier si le point se trouve dans la région centrale, générez un objet fictif.
L'algorithme complet de la ziggourat pour les distributions unilatérales est donc :
Choisissez un calque aléatoire
Laisser
Si
Si
Laisser
Calculer
Sinon, choisissez de nouveaux nombres aléatoires et retournez à l'étape 1.
Pour une distribution bilatérale, le résultat doit être inversé dans 50 % des cas. Cela peut souvent se faire facilement en choisissant
Algorithmes de repli pour la queue
L'algorithme de la ziggourat ne génère la plupart des résultats que très rapidement et nécessite un algorithme de secours en cas de besoin.x_{1 Elle est toujours plus complexe qu'une implémentation plus directe. L'algorithme de repli spécifique dépend de la distribution.
Pour une distribution exponentielle, la queue ressemble trait pour trait au corps de la distribution. Une solution consiste à utiliser l'algorithme le plus élémentaire E = −ln( U₁ ) et à poser x = x₁ − ln ( U₁ ) . Une autre consiste à appeler l'algorithme de la ziggourat de manière récursive et à ajouter x₁ au résultat.
Pour une distribution normale, Marsaglia propose un algorithme compact :
Soit x = −ln( U 1 )/ x 1 .
Soit y = −ln( U 2 ).
Si 2 y > x 2 , retournez x + x 1 .
Sinon, retournez à l'étape 1.
Étant donné que x₁ ≈ 3,5 pour des tables de taille courante, le test de l'étape 3 est presque toujours concluant. Comme −ln( U₁ ) suit une loi exponentielle, on peut utiliser une implémentation de cette loi.
Optimisations
L'algorithme peut être exécuté efficacement avec des tables précalculées de x i et y i = f ( x i ), mais il existe quelques modifications pour le rendre encore plus rapide :
Rien dans l'algorithme de la ziggourat ne dépend de la normalisation de la fonction de distribution de probabilité (intégrale sous la courbe égale à 1), la suppression des constantes de normalisation peut accélérer le calcul de f ( x ).
La plupart des générateurs de nombres aléatoires uniformes sont basés sur des générateurs de nombres aléatoires entiers qui renvoient un entier dans la plage [0, 2 32 − 1 ]. Une table de 2 −32 x i nous permet d'utiliser directement de tels nombres pour U 0 .
Lors du calcul de distributions bilatérales à l'aide d'un U 0 bilatéral comme décrit précédemment, l'entier aléatoire peut être interprété comme un nombre signé dans la plage [−2 31 , 2 31 − 1], et un facteur d'échelle de 2 −31 peut être utilisé.
Au lieu de comparer U₀xᵢ à xᵢ₊₁ à l' étape 3, il est possible de précalculer xᵢ₊₁ / xᵢ et de comparer directement U₀ à cette valeur . Si U₀ est un générateur de nombres aléatoires entiers, ces limites peuvent être prémultipliées par 2³² ( ou 2³¹ , selon le cas) afin d'effectuer une comparaison d'entiers .
Avec les deux modifications ci-dessus, le tableau des valeurs x i non modifiées n'est plus nécessaire et peut être supprimé.
Lors de la génération de valeurs à virgule flottante simple précision IEEE 754 , dont la mantisse est de 24 bits (incluant le 1 implicite en tête), les bits de poids faible d'un nombre aléatoire entier de 32 bits ne sont pas utilisés. Ces bits peuvent servir à sélectionner le numéro de couche. (Voir les références ci-dessous pour plus de détails.)
Les trois premières étapes peuvent être placées dans une fonction en ligne , qui peut appeler une implémentation hors ligne des étapes moins fréquemment nécessaires.
Générer les tables
Il est possible de stocker l'ensemble du tableau précalculé, ou simplement d'inclure les valeurs n , y1 , A et une implémentation de f −1 ( y ) dans le code source , et de calculer les valeurs restantes lors de l'initialisation du générateur de nombres aléatoires.
Lors du remplissage des valeurs du tableau, supposez simplement que x n = 0 et y n = f (0), et acceptez la légère différence dans la surface de la couche n − 1 comme une erreur d'arrondi.
Trouver x 1 et A
Étant donné une valeur initiale (estimée) de x₁ , nous avons besoin d'un moyen de calculer l'aire t de la queue pour laquelle x > x₁ . Pour la distribution exponentielle, il s'agit simplement de e⁻ˣ₁ , tandis que pour la distribution normale, en supposant que nous utilisons la fonction non normalisée f ( x ) = e⁻ˣ² / 2 , il s'agit de … erfc ( x /) . Pour des distributions plus complexes, une intégration numérique peut être nécessaire.
Avec cela en main, à partir de x 1 , nous pouvons trouver y 1 = f ( x 1 ), l'aire t dans la queue et l'aire de la couche de base A = x 1 y 1 + t .
Calculez ensuite les séries y <sub>i</sub> et x <sub>i </sub> comme précédemment. Si y <sub>i </sub> > f (0) pour tout i < n , alors l'estimation initiale x <sub>1</sub> était trop faible, ce qui conduit à une aire A trop grande . Si y <sub>n </sub> < f (0), alors l'estimation initiale x<sub> 1</sub> était trop élevée.
Dans ce cas , utilisez un algorithme de recherche de racines ( comme la méthode de dichotomie ) pour trouver la valeur x₁ qui rapproche au maximum yₙ₋₁ de f (0). Alternativement , cherchez la valeur qui rend l'aire de la couche supérieure, xₙ₋₁ ( f (0) − yₙ₋₁ ), aussi proche que possible de la valeur souhaitée A. Cela permet d'économiser une évaluation de fₙ₋₁ ( x ) et constitue en réalité la condition la plus pertinente.La variante de McFarland
Christopher D. McFarland a proposé une version encore plus optimisée. Celle-ci applique trois modifications algorithmiques, au prix de tables légèrement plus grandes.
Dans un premier temps, le cas général ne considère que les portions rectangulaires, de (0, y i −1 ) à ( x i , y i ). Les régions de forme irrégulière à droite de celles-ci (principalement presque triangulaires, plus la queue) sont traitées séparément. Cela simplifie et accélère le chemin rapide de l'algorithme .
Deuxièmement, on utilise la surface exacte des régions de forme irrégulière ; elle n’est pas arrondie pour inclure le rectangle entier de ( x <sub>i -1</sub> , y <sub>i</sub> ). Cela augmente la probabilité que le chemin le plus rapide soit emprunté.
Il en résulte notamment que le nombre de couches est légèrement inférieur à n . Bien que l'aire des portions de forme irrégulière soit calculée avec précision, le total dépasse l'aire d'une seule couche. L'aire par couche est ajustée afin que le nombre de couches rectangulaires soit un entier. Si la valeur initiale de i (0 ≤ i < n) est supérieure au nombre de couches rectangulaires, la phase 2 se poursuit.
Si la valeur recherchée se situe dans une zone de forme irrégulière, la méthode des alias est utilisée pour en sélectionner une, en fonction de sa surface réelle. Cela représente un léger surcroît de travail et nécessite des tables d'alias supplémentaires, mais permet de choisir l'une des faces droites des couches.
La région de forme irrégulière sélectionnée est soumise à un échantillonnage par rejet. Si un échantillon est rejeté, l'algorithme ne retourne pas au début. L'aire réelle de chaque région de forme irrégulière a servi à choisir une couche ; la boucle d'échantillonnage par rejet reste donc dans cette couche jusqu'à ce qu'un point soit sélectionné.
Troisièmement, la forme quasi triangulaire de la plupart des portions de forme irrégulière est exploitée, bien que cela doive être divisé en trois cas en fonction de la dérivée seconde de la fonction de distribution de probabilité dans la couche sélectionnée.
Si la fonction est convexe (comme la distribution exponentielle et la distribution normale pour | > 1), elle est strictement contenue dans le triangle inférieur. Deux variables aléatoires uniformes unitaires, U₁ et U₂ , sont choisies, et avant d'être mises à l'échelle sur le rectangle englobant la région de forme irrégulière, leur somme est testée. Si U₁ + U₂ > 1, le point se trouve dans le triangle supérieur et peut être symétrisé en (1 − U₁ , 1 − U₂ ) . Ensuite , si U₁ + U₂ < 1 − ε , pour une tolérance ε appropriée , le point est assurément en dessous de la courbe et peut être immédiatement accepté. Ce n'est que pour les points très proches de la diagonale qu'il est nécessaire de calculer la fonction de répartition f ( x ) pour effectuer un test de rejet exact. (La tolérance ε devrait en théorie dépendre de la couche, mais une seule valeur maximale peut être utilisée sur toutes les couches avec une perte négligeable.)
Si la fonction est concave (comme l'est la distribution normale pour | < 1), elle inclut une petite partie du triangle supérieur, la réflexion est donc impossible, mais les points dont les coordonnées normalisées satisfont peuvent être immédiatement acceptés, et les points pour lesquels peuvent être immédiatement rejetés.
Dans la couche qui chevauche | = 1, la distribution normale a un point d'inflexion et le test de rejet exact doit être appliqué si 1− ε < U 1 + U 2 < 1+ ε .
La queue est traitée comme dans l'algorithme Ziggurat original, et peut être considérée comme un quatrième cas pour la forme de la région de forme irrégulière à droite.
.Ce document numérote les couches à partir de 1 en partant du haut, et considère la couche 0 en bas comme un cas particulier, tandis que l'explication ci-dessus numérote les couches à partir de 0 en bas.
Implémentation en C de la méthode de la ziggourat pour la fonction de densité normale et la fonction de densité exponentielle , reprenant essentiellement le code présenté dans l'article. (Attention : ce code C suppose des entiers 32 bits.)
Implémentation AC# de l'algorithme de la ziggourat et aperçu de la méthode.
Jurgen A. Doornik (2005). « Une méthode de ziggourat améliorée pour générer des échantillons aléatoires normaux » (PDF) . Nuffield College, Oxford . Décrit les risques liés à l'utilisation des bits de poids faible du générateur de nombres aléatoires entiers pour choisir le numéro de couche.
Comportement normal par Cleve Moler, MathWorks, décrivant l'algorithme de ziggourat introduit dans la version 5 de MATLAB , 2001.
Blogs de MathWorks sur le générateur de nombres aléatoires normaux Ziggurat , publiés par Cleve Moler, le 18 mai 2015.
David B. Thomas, Philip HW Leong, Wayne Luk et John D. Villasenor (octobre 2007). « Générateurs de nombres aléatoires gaussiens » (PDF) . ACM Computing Surveys , vol. 39 , n° 4, p. 11, 1-38. doi : 10.1145/1287620.1287622 . ISSN 0360-0300 . S2CID 10948255. Lorsque le maintien d’une qualité statistique extrêmement élevée est la priorité absolue, et que, sous cette contrainte, la rapidité est également souhaitée, la méthode Ziggurat est souvent le choix le plus approprié.Comparaison de plusieurs algorithmes de génération de nombres aléatoires gaussiens .
Nadler, Boaz (2006). « Défauts de conception dans l’implémentation des méthodes Ziggurat et Monty Python (et quelques remarques sur Matlab randn) ». arXiv : math/0603058 .. Illustre les problèmes liés aux générateurs de nombres pseudo-aléatoires uniformes sous-jacents et comment ces problèmes affectent la sortie de l'algorithme ziggourat.
Edrees, Hassan M.; Cheung, Brian; Sandora, McCullen; Nummey, David; Stefan, Deian (13-16 juillet 2009). Algorithme de ziggourat optimisé matériellement pour les générateurs de nombres aléatoires gaussiens à haute vitesse (PDF) . Conférence internationale 2009 sur l'ingénierie des systèmes et algorithmes reconfigurables. Las Vegas. original ) le 24 janvier 2013 .
Marsaglia, George (septembre 1963). Génération d'une variable à partir de la queue de la distribution normale (Rapport technique). Laboratoires de recherche scientifique de Boeing. Note mathématique n° 322, numéro d'accès DTIC : AD0423993. Archivé de l'original le 10 septembre 2014 – via le Defense Technical Information Center .
Worldlex Wiki 
