Méthode de rejet

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

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

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

Autrement dit,

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

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

La méthode de rejet

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

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

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

1. Boucler:
   • Tirer ${\displaystyle Y}$ de densité ${\displaystyle g;}$
   • Tirer ${\displaystyle U}$ selon la loi uniforme U(0;1), indépendamment de ${\displaystyle Y;}$
2. Tant que ${\displaystyle U>{\tfrac {f(Y)}{c\,g(Y)}},}$ reprendre en 1;
3. Accepter ${\displaystyle Y}$ comme un tirage aléatoire de densité de probabilité ${\displaystyle 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 ${\displaystyle N,}$ est donc lui-même aléatoire. On peut montrer que ${\displaystyle N}$ suit la loi géométrique de paramètre ${\displaystyle {\frac {1}{c}},}$ c'est-à-dire

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

En effet,

${\displaystyle {\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)}}\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 ${\displaystyle N}$(c.-à-d. le nombre moyen d'itérations à effectuer avant d'obtenir une réalisation de la densité f ) vaut ${\displaystyle c}$.

${\displaystyle \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:

${\displaystyle 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 ${\displaystyle cg-f\geq 0,}$

${\displaystyle 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.

Le fait que ${\displaystyle f(x)\leq cg(x)}$ peut être écrit comme ${\displaystyle 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 ${\displaystyle 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:

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

tout en supposant qu'il existe une densité g telle que ${\displaystyle f(x)\leq cg(x)}$. Aucune autre hypothèse n'est nécessaire; en particulier, il ne faut pas imposer que ${\displaystyle 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(${\displaystyle Z\leq h_{1}(Y)}$, vrai, faux );
   • Si Suite alors
     • Si ${\displaystyle Z\leq h_{2}(Y)}$ alors Suite:= SI(${\displaystyle 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.