Probabilités de l'ingénieur : algorithmes stochastiques et simulation/Méthode de rejet

Début de la boite de navigation du chapitre

fin de la boite de navigation du chapitre

Principe

L'idée principale de la méthode de rejet est d’utiliser le tirage d'une variable aléatoire "simple" à simuler (par la méthode de la transformée inverse par exemple) pour simuler une fonction à densité plus complexe (avec de fortes irrégularités par exemple) ou dans des cas où la transformée inverse ne peut pas s'appliquer.

La méthode de rejet repose sur cette proposition :

Début d’un théorème

Théorème

Soit f une fonction de densité de probabilités. On suppose qu’il existe une densité de probabilités g telle que :

$\exists K>0\ ,\ \forall x\in \mathbb {R} \ ,\ f(x)\leq Kg(x)$

Soit alors Z suivant la loi de densité g, et $Y\sim {\mathcal {U}}([0;Kg(Z)])$.

Alors la variable aléatoire $X=\lbrace Z|Y\leq f(Z)\rbrace$ suit la loi de densité f.

Fin du théorème

Début d'une démonstration

Démonstration

On va montrer qu’il y a égalité des fonctions de répartiton : $\mathbb {P} (X\leq x_{0})=F(x_{0})$

$\mathbb {P} (X\leq x_{0})={\frac {\mathbb {P} \left[(X\leq x_{0})\cap (Y\leq f(X))\right]}{\mathbb {P} (Y\leq f(X))}}$.

Or, $\mathbb {P} \left[(X\leq x_{0})\cap (Y\leq f(X))\right]=\int _{-\infty }^{x_{0}}\mathbb {P} \left[Kg(x)U\leq f(x)\right]g(x)dx=\int _{-\infty }^{x_{0}}\mathbb {P} \left[U\leq {\frac {f(x)}{Kg(x)}}\right]g(x)dx=\int _{-\infty }^{x_{0}}{\frac {f(x)}{Kg(x)}}g(x)dx={\frac {F(x_{0})}{K}}$

et $\mathbb {P} (Y\leq f(X))=\int _{-\infty }^{+\infty }\mathbb {P} \left[Kg(x)U\leq f(x)\right]g(x)dx={\frac {1}{K}}$

D'où le résultat voulu.

Fin de la démonstration

Ainsi, connaissant la densité, on construit un algorithme qui permet de simuler la variable aléatoire de densité connue :

Algorithme de rejet

On suppose f, g et K connus.

On tire X suivant la loi de g et $Y\sim {\mathcal {U}}([0;Kg(X)])$

1. Si $Y\leq f(X)$, on garde X.
2. Sinon, on tire un nouveau couple $(X,Y)$ jusqu'à ce que la condition 1. soit remplie.

Évaluation du nombre de tirages

On peut évaluer le nombre de tirages nécessaires T afin d'obtenir un point acceptable. En effet, remarquons déjà :

$\int _{\mathbb {R} }(Kg-f)=K\int _{\mathbb {R} }g-\int _{\mathbb {R} }f=K-1\geq 0$

car f et g sont des fonctions densités, donc d'intégrales égales à 1.

Instinctivement, on voit que plus la constante K sera grande, plus le nombre de tirages nécessaires sera important (la zone de rejet sera d'autant plus grande). De plus, ce nombre de tirages est une variable aléatoire (on ne peut pas prévoir à chaque fois combien de tirages seront faits). Il suffit donc de trouver la loi de T pour pouvoir faire des prévisions.

Propriété

Dans un algorithme de rejet, le nombre de tirages nécessaires à chaque étape suit une loi géométrique de paramètre ${\frac {1}{K}}$

(Preuve à faire)

Ainsi, le nombre de tirages moyen est de l’ordre de l'espérance, soit K. Il convient donc de choisir la constante la plus proche de K possible.

Applications

Cas d'une densité à support compact

Pour le cas où la densité aurait une densité à support compact ${\textrm {Supp}}(f)=[a,b]$, le plus simple est de considérer la fonction densité majorante comme la densité de la loi uniforme sur [a,b], et de prendre $K=(b-a)\max _{x\in [a,b]}f(x)$

Sur l'illustration à droite, la zone d'acceptance $Y\leq \max _{x\in [a,b]}f(x)$ est sous la courbe bleue, et la zone de rejet est entre la courbe bleue et la droite rouge.

Cas général

Dans le cas où la fonction densité aurait un support infini (loi normale, loi log-normale), il faut chercher une densité de même support qu'on peut facilement simuler. Les exemples naturels sont les lois exponentielle (définie sur $\mathbb {R} _{+}$) et Laplace ou de Cauchy (définies sur $\mathbb {R}$), qu'on peut simuler par la méthode de la transformée inverse.