Méthode de rejet

La méthode du rejet est une méthode utilisée dans le domaine des probabilités.

But

La méthode de rejet est utilisée pour engendrer indirectement une variable aléatoire $X$, de densité de probabilité $f,$ lorsqu'on ne sait pas simuler directement la loi de densité de probabilité $f$ (c'est le cas par exemple si $f$ n'est pas une densité classique, mais aussi pour la loi de Gauss^{[réf. nécessaire]}).

Soit $(Y,U)$ un couple de variables aléatoires indépendantes tirées selon une loi uniforme, i.e. $(Y,U)$ est un point tiré uniformément dans le carré unité. On peut alors montrer que la distribution de $X$ est la loi conditionnelle de $Y$ sachant l'événement

$M=\{U\leq f_{X}(Y)\}.$

Autrement dit,

${f_{X}(x)=f_{Y}(x|M)}.$

Pour simuler une suite de variables aléatoires réelles $\left(X_{n}\right)_{n\geq 1}$ de distribution identique à celle de $X,$ il suffit donc, dans une suite de tirages de couples $(Y_{i},U_{i})$ uniformes indépendants, de sélectionner les $Y_{i}$ correspondant aux tirages $(Y_{i},U_{i})$ vérifiant $M,$ et de rejeter les autres.

Algorithme

La méthode de rejet

On voudrait simuler une variable aléatoire réelle $X$ de densité de probabilité $f$. On suppose

• qu'il existe une autre densité de probabilité $g$ telle que le ratio ${\frac {f}{g}}$ soit borné, disons par $c$ (i.e. $f\leq cg$),
• qu'on sache simuler $Y$ de densité $g.$

La version basique de la méthode de rejet prend la forme suivante:

1. Boucler:
   • Tirer $Y$ de densité $g;$
   • Tirer $U$ selon la loi uniforme U(0;1), indépendamment de $Y;$
2. Tant que $U>{\tfrac {f(Y)}{c\,g(Y)}},$ reprendre en 1;
3. Accepter $Y$ comme un tirage aléatoire de densité de probabilité $f.$

On remarque que l'algorithme comporte une boucle dont la condition porte sur des variables aléatoires. Le nombre d'itérations, notons-le $N,$ est donc lui-même aléatoire. On peut montrer que $N$ suit la loi géométrique de paramètre ${\frac {1}{c}},$ c'est-à-dire

$\mathbb {P} \left(N=k\right)\ =\ \left(1-{\tfrac {1}{c}}\right)^{k-1}\,{\tfrac {1}{c}}\ 1_{k\geq 1}.$

En effet,

${\frac {1}{c}}\ =\ \int _{\mathbb {R} }\,{\tfrac {f(y)}{c\,g(y)}}\ g(y)\,dy\ =\ \int _{\mathbb {R} }\,\mathbb {P} \left(U\leq {\tfrac {f(y)}{c\,g(y)}}\mid Y=y\right)\ g(y)\,dy\ =\ \mathbb {P} \left(U\leq {\tfrac {f(Y)}{c\,g(Y)}}\right)$

est la probabilité, lors d'une itération, de terminer la boucle, et, par conséquent, d'accepter Y. Par suite, l'espérance de $N$(c.-à-d. le nombre moyen d'itérations à effectuer avant d'obtenir une réalisation de la densité f ) vaut $c$.

$\mathbb {E} \left[N\right]\ =\ c.$

On a donc tout intérêt à choisir c le plus petit possible. En pratique, une fois la fonction g choisie, le meilleur choix de c est donc la plus petite constante qui majore le ratio f/g, c'est-à-dire:

$c=\sup {\frac {f(x)}{g(x)}}.$

Notons que, soit c est supérieur strict à 1, soit f=g, la deuxième solution étant assez théorique : en effet, comme $cg-f\geq 0,$

$0\ \leq \ \int _{\mathbb {R} }\left(cg-f\right)\ =\ c\ \int _{\mathbb {R} }g\ -\ \int _{\mathbb {R} }f\ =\ c-1.$

On a donc intérêt à choisir c le plus proche de 1 possible, pour que le nombre d'itérations moyen soit proche de 1 lui aussi. Bref, le choix de l'enveloppe g est primordial:

• le tirage de la loi g doit être facile ;
• l'évaluation de f(x)/g(x) doit être aisée ;
• la constante c doit être la plus petite possible ;
• la fonction cg doit majorer la densité f.

Les deux derniers points conduisent à rechercher une fonction g dont le graphe "épouse" étroitement celui de f.

Généralisations

Le fait que $f(x)\leq cg(x)$ peut être écrit comme $f(x)=cg(x)h(x)$ où h est une fonction à valeurs dans [0;1]. On remplace l'étape 2 de l'algorithme initial par la condition:

Tant que $U/h(X)>1$, reprendre en 1

Une autre généralisation peut être considérée lorsque l'évaluation du ratio f/g est délicate. On cherche alors à encadrer la fonction f par deux fonctions facilement évaluables:

$h_{1}(x)\leq f(x)\leq h_{2}(x)$

tout en supposant qu'il existe une densité g telle que $f(x)\leq cg(x)$. Aucune autre hypothèse n'est nécessaire; en particulier, il ne faut pas imposer que $h_{2}(x)\leq cg(x)$. L'algorithme prend alors la forme suivante:

1. Suite := vrai
2. Tant que Suite
   • tirer Y selon g;
   • tirer U selon la loi uniforme U(0;1), indépendamment de Y;
   • Z := U c g(Y);
   • Suite := SI($Z\leq h_{1}(Y)$, vrai, faux );
   • Si Suite alors
     • Si $Z\leq h_{2}(Y)$ alors Suite:= SI($Z\leq f(Y)$,vrai,faux) fin si
   • Fin si
3. fin tant que
4. retourne Y comme un tirage de f.

Dans cet algorithme, les fonctions h permettent de ne recourir à une comparaison à f (et donc à son évaluation) que très rarement.

Voir aussi

Bibliographie

• (en) Robert, C.P. et George Casella, Monte Carlo Statistical Methods, New York, Springer-Verlag, 2004.
• (en) J. von Neumann, « Various techniques used in connection with random digits. Monte Carlo methods », Nat. Bureau Standards, vol. 12,‎ 1951, p. 36–38.
• (en) Luc Devroye, Non-Uniform Random Variate Generation, New York, Springer-Verlag, 1986 (lire en ligne), « 2, section 3, p. 40 »
• (en) Sheldon Ross, A First Course in Probability, 8, 2002, 442-7 p. (lire en ligne)

Articles connexes

• Méthode de la transformée inverse et son graphe
• Algorithme de Metropolis-Hastings