MĂ©thode ziggourat
La mĂ©thode ziggourat est un algorithme pour engendrer des nombres alĂ©atoires de loi non uniforme. Il s'agit d'une mĂ©thode de rejet et peut ĂȘtre choisie pour simuler une variable alĂ©atoire ayant une densitĂ© strictement monotone. Cette mĂ©thode peut aussi ĂȘtre appliquĂ©e Ă des lois symĂ©triques unimodales telles que la loi normale en choisissant un point sur l'une des moitiĂ©s et en choisissant le cĂŽtĂ© alĂ©atoirement. Cette mĂ©thode a Ă©tĂ© dĂ©veloppĂ©e par George Marsaglia et Wai Wan Tang.
Comme la plupart des algorithmes de ce type, il suppose que l'on dispose d'un gĂ©nĂ©rateur de nombres alĂ©atoires de loi uniforme, en gĂ©nĂ©ral un gĂ©nĂ©rateur pseudo-alĂ©atoire. Dans la plupart des cas, comme la loi normale ou la loi exponentielle, l'algorithme consiste Ă engendrer un nombre flottant, un index alĂ©atoire de table, faire une recherche dans une table, une multiplication et une comparaison. Cet algorithme est considĂ©rablement plus rapide que les autres mĂ©thodes pour simuler des variables alĂ©atoires de loi normale, comme la mĂ©thode de Box-Muller qui requiĂšrent de calculer une racine carrĂ©e ou un logarithme. D'un autre cĂŽtĂ©, cette mĂ©thode est plus complexe Ă mettre en Ćuvre et nĂ©cessite des tables prĂ©calculĂ©es, de sorte qu'il vaut mieux l'utiliser lorsque l'on a besoin de nombres alĂ©atoires en grande quantitĂ©.
Le nom de cette méthode provient du fait qu'elle couvre la densité de la loi avec des segments rectangulaires empilés par ordre de taille décroissant, ce qui ressemble donc à une ziggourat.
Théorie
La méthode ziggourat est une méthode de rejet. La position d'un point est engendrée aléatoirement à l'intérieur d'une courbe délimitant une surface légÚrement plus grande que celle délimitée par la densité de la loi considérée. On teste alors si ce point se trouve dans cette derniÚre surface. Si c'est le cas, on retourne l'abscisse du point. Sinon, on rejette ce point et on en tire un nouveau.
Pour construire la surface qui recouvre la densité, on choisit une surface composée de n régions de tailles égales, dont n-1 sont des rectangles empilés au-dessus d'une région non rectangulaire qui couvre la queue de la densité de la loi.
Si f est la densitĂ© de la loi Ă simuler, la base de la ziggourat se compose donc d'un rectangle allant dont le coin infĂ©rieur gauche a pour coordonnĂ©es (0 ; 0) et le coin supĂ©rieur a pour coordonnĂ©es (x0 , f(x0)), et de l'ensemble des points situĂ©s sous la courbe (x , f(x)) pour x â„ x0. Cette couche dĂ©limite une rĂ©gion d'aire A. On place au-dessus une rĂ©gion rectangulaire dont le coin infĂ©rieur gauche est (0 , f(x0)) et le coin supĂ©rieur droit (x1 , f(x1)), oĂč x1 est pour que l'aire de ce rectangle soit Ă©gale Ă A. On construit de mĂȘme un rectangle de coordonnĂ©es dĂ©finies par (0 , f(x1)) et (x2 , f(x2)) d'aire A. On construit ainsi une suite de points x0, x1,..., xn pour un nombre n de couches donnĂ© Ă l'avance (typiquement, n-256). Les coordonnĂ©es des points x0, x1,..., xn dĂ©pendent de n.
En ignorant pour l'instant le problĂšme de la premiĂšre couche qui n'est pas rectangulaire, l'algorithme est le suivant :
- On choisit aléatoirement de façon une couche i avec une probabilité 1/n.
- Si i = 0, on utilise alors un algorithme spécifique (voir plus bas)
- On se donne une réalisation U0 d'une variable aléatoire uniforme sur [0 ; 1[
- On pose x' = U0xi-1
- Si , alors on retourne x'
- Sinon on se donne une réalisation U1 d'une variable aléatoire uniforme sur [0 ; 1[
- On calcule avec
- Si , on retourne x'
- Sinon, on recommence l'algorithme depuis le début.
Pour une loi dont la densité est symétrique, il suffit juste de tester si en utilisant pour U0 une variable aléatoire uniforme sur ]-1 ; 1[.
![](https://img.franco.wiki/i/Ziggurat_method.gif)
Algorithmes pour la queue de la densité
Lorsque la densité à un support non compact (comme c'est le cas pour la loi normale, la loi exponentielle...), il est alors nécessaire d'utiliser un algorithme spécifique lorsque la premiÚre tranche est sélectionnée. Cet algorithme dépend de la loi.
La premiĂšre couche peut se diviser en une rĂ©gion centrale rectangulaire, et une rĂ©gion avec une queue infinie. On peut utiliser le mĂȘme algorithme en ajoutant un point fictif et en tirant un point oĂč U0 une rĂ©alisation d'une variable alĂ©atoire uniforme sur [0 ; 1 [. Si , on retourne x. Sinon, on a besoin de tirer une rĂ©alisation de cette variable alĂ©atoire sachant qu'elle est plus grande que la valeur x0. Pour la loi exponentielle de paramĂštre λ, cela se fait trĂšs facilement par oĂč U est la rĂ©alisation d'une variable alĂ©atoire uniforme sur [0 ; 1 [.
Une autre possibilité consiste à appeler l'algorithme ziggourat de façon récursive sur la queue de la densité.
Pour la loi normale, G. Marsaglia[1] propose la méthode suivante :
- On calcule x = -ln(U1)/x0 oĂč U1 est la rĂ©alisation d'une variable alĂ©atoire uniforme sur [0 ; 1 [.
- On calcule y = -ln(U2) oĂč U2 est la rĂ©alisation d'une variable alĂ©atoire uniforme sur [0 ; 1 [.
- Si 2y > x2 alors on retourne x + x0,
- Sinon, on revient au premier pas.
Comme pour une taille de table typique (n = 128), x0 â 3,5, le test du pas no 3 est presque toujours vĂ©rifiĂ©. Puisque -ln(Ui) renvoie une rĂ©alisation d'une variable alĂ©atoire exponentielle, si l'on dispose d'une mĂ©thode ziggourat pour l'exponentielle, on peut l'utiliser Ă cette fin.
Optimisations
L'algorithme peut tourner de façon efficace en utilisant des tables prĂ©calculĂ©es pour les xi et les yi, mais il existe des astuces pour le rendre encore plus rapide. La principale idĂ©e est que rien ne requiert que l'on utilise pour f une fonction d'aire 1. On peut donc utiliser cet algorithme sur des densitĂ©s non normalisĂ©es, ce qui accĂ©lĂšre le calcul de f(x) (par exemple, pour la loi normale, on peut utiliser f(x) = exp(-x2/2) au lieu de f(x) = exp(-x2/2)/â2Ï ).
Comment engendrer les tables ?
Il est possible de stocker dans l'algorithme la table entiĂšre des points xi et , ou bien juste les valeurs n, A et x0 et de calculer les autres valeurs juste durant la phase d'initialisation des calculs.
Comment trouver x0 et A ?
Ătant donnĂ© un nombre de niveau n, on a besoin de trouver la valeur de x0 et de A telles que et . Pour la loi exponentielle de paramĂštre λ, on a l'aire . La fonction de rĂ©partition de la loi normale est trĂšs souvent incluse dans de nombreuses bibliothĂšques numĂ©riques. Pour d'autres lois, il peut ĂȘtre nĂ©cessaire d'utiliser des procĂ©dures d'intĂ©gration numĂ©rique. L'idĂ©e est alors d'utiliser un algorithme de recherche de zĂ©ro en partant d'une estimation initiale de x0 afin de rĂ©soudre
en calculant les successivement jusqu'à de sorte que l'aire de chaque tranche soit égale à pour une estimation donnée de x0.
Notes et références
- (en) George Marsaglia, « Generating a Variable from the Tail of the Normal Distribution », ASQC and the American Statistical Association,â (lire en ligne [PDF])
Voir aussi
Articles connexes
Bibliographie
- (en) George Marsaglia et Wai Wan Tsang, « The Ziggurat Method for Generating Random Variables », dans Journal of Statistical Software, volume 5, no 8, 2000. DOI 10.18637/jss.v005.i08
- (en) Jurgen A. Doornik, « An Improved Ziggurat Method to Generate Normal Random Samples », Nuffield College, Oxford.
- (en) David B. Thomas, Philip G.W. Leong, Wayne Luk et John D. Villasenor, « Gaussian Random Number Generators », dans ACM Computing Surveys, volume 39, 2007. DOI 10.1145/1287620.1287622