1.5. Méthode du rejet

Soit \((E, \mathcal{E})\) un espace mesurable et \(\mu\) une mesure \(\sigma\)–finie positive de référence sur \(E\). Typiquement \(E = \mathbf{R}^d\), la tribu \(\mathcal{E}\) est la tribu borélienne et \(\mu\) est la mesure de Lebesgue sur \(\mathbf{R}^d\). Dans ce cadre abstrait on peut aussi traiter les variables aléatoires discrètes en prenant la mesure de comptage pour mesure de référence \(\mu\). On veut simuler selon une loi \(\nu\) mesure de probabilité absolument continue par rapport à \(\mu\) de densité de Radon-Nykodim \(f\) i.e.

\[\begin{equation*} \exists f:(E, \mathcal{E}) \to \mathbf{R}_+ \quad \text{t.q.} \quad \nu = f \mu, \end{equation*}\]

et \(\int f \operatorname{d}\! \mu = 1\) (car \(\nu\) est une probabilité). On rappelle pour \(X\) une variable aléatoire de loi \(\nu\) on a

\[\begin{equation*} \forall \phi \in \mathcal{C}_b, \quad \mathbf{E}[\phi(X)] = \int \phi(x) f(x) \mu(\operatorname{d}\! x) \end{equation*}\]

On suppose qu’on sait simuler selon une loi auxiliaire \(\bar \nu\) qui est absolument continue par rapport à \(\mu\) et de densité \(g\) vérifiant \(g > 0\) \(\mu-p.p.\) et

(1.5)\[ \exists c \ge 1, \, \forall x \in \mathbf{R}^d, \quad f(x) \le c g(x). \]

Le résultat suivant permet d’obtenir des réalisations de la loi \(\nu\) à partir de réalisations de la loi auxiliaire \(\bar \nu\).

Theorem 1.2

Soit \((U_n, Y_n)_{n \ge 1}\) une suite de v.a. i.i.d. de loi \(\mathcal{U}(]0,1[) \otimes \bar \nu\) (où \(\bar \nu\) est de densité \(g > 0\) \(\mu-p.p.\) par rapport à \(\mu\) et vérifiant (1.5)). Soit

()\[\begin{equation} \tau = \min\big\{k \ge 1, c U_k g(Y_k) < f(Y_k) \big\}. \end{equation}\]

Alors \(\tau\) suit une loi géométrique de paramètre \(p = 1/c\) et \(Y_\tau\) a pour loi \(\nu\) (de densité \(f\) par rapport à \(\mu\)).

Proof. Le fait que \(\tau\) suit une loi géométrique est vrai par définition de la loi géométrique en considérant l’expérience de Bernoulli \((B_n)_{n \ge 1}\) avec \(B_n = \mathbf{1}_{\{c U_n g(Y_n) < f(Y_n)\}}\). Les \((B_n)_{n \ge 1}\) sont bien i.i.d. de loi de Bernoulli de paramètre \(p = \mathbf{P}\big(c U_1 g(Y_n) < f(Y_n) \big)\). On calcule \(p\) dans la suite.

Tout d’abord on prouve que \(Y_\tau\) a même loi que la loi conditionnelle \(Y_1 | \{c U_1 g(Y_1) < f(Y_1)\}\) c’est à dire que pour tout \(B \in \mathcal{E}\) on a

\[\begin{equation*} \mathbf{P} \big(Y_\tau \in B\big) = \mathbf{P} \big( Y_1 \in B | c U_1 g(Y_1) < f(Y_1) \big). \end{equation*}\]

En utilisant la partition \((\{\tau = n\}, n \ge 1)\) de \(\Omega\) on a

\[\begin{align*} \mathbf{P} \big(Y_\tau \in B \big) &= \mathbf{P} \Big(\cup_{n \ge 1} \big(\{\tau = n\} \cap \{Y_n \in B\} \big) \Big), \\ & = \sum_{n \ge 1} \mathbf{E} \big[ \mathbf{1}_{\{\tau = n\}} \mathbf{1}_{\{Y_n \in B\}} \big]. \end{align*}\]

Par définition \(\{\tau = n\} = \big\{c U_1 g(Y_1) \ge f(Y_1),\dots, c U_{n-1} g(Y_{n-1}) \ge f(Y_{n-1}), c U_n g(Y_n) < f(Y_n)\big\}\) et comme la suite \((U_n, Y_n)_{n \ge 1}\) est i.i.d. on a

\[\begin{align*} \mathbf{P} \big(Y_\tau \in B \big) = \sum_{n \ge 1} \mathbf{P} \big(c U_1 g(Y_1) \ge f(Y_1) \big)^{n-1} \mathbf{P} \big(c U_1 g(Y_1) < f(Y_1), Y_1 \in B \big). \end{align*}\]

En posant \(p = \mathbf{P} \big(c U_1 g(Y_1) < f(Y_1) \big)\) et en utilisant \(\sum_{p \ge 1} (1-p)^{n-1} p = 1\) on a le résultat

\[\begin{equation*} \mathbf{P} \big(Y_\tau \in B \big) = \frac{1}{p} \mathbf{P} \big(c U_1 g(Y_1) < f(Y_1), Y_1 \in B \big) = \mathbf{P} \big(Y_1 \in B | c U_1 g(Y_1) < f(Y_1)\big). \end{equation*}\]

On détermine maintenant la loi conditionnelle \(Y_1 | \{c U_1 g(Y_1) < f(Y_1)\}\). Soit \(\phi:\mathbf{R}^d \to \mathbf{R}\) une fonction continue bornée, alors

\[\begin{align*} \mathbf{E} \big[ \phi(Y_1) \mathbf{1}_{\{c U_1 g(Y_1) < f(Y_1)\}} \big] = \int_{\mathbf{R}^d} \phi(y) \Big( \int_0^1 \mathbf{1}_{c u g(y) < f(y)} \operatorname{d}\! u \Big) g(y) \mu(\operatorname{d}\! y). \end{align*}\]

Le point clé est de remarquer que \(\int_0^1 \mathbf{1}_{c u g(y) < f(y)} \operatorname{d}\! u = \frac{f(y)}{c g(y)}\) donc

\[\begin{align*} \mathbf{E} \big[ \phi(Y_1) \mathbf{1}_{\{c U_1 g(Y_1) < f(Y_1)\}} \big] &= \int_{\mathbf{R}^d} \phi(y) \frac{f(y)}{c g(y)} g(y) \mu(\operatorname{d}\! y), \\ & = \frac{1}{c} \int_{\mathbf{R}^d} \phi(y) f(y) \mu(\operatorname{d}\! y). \end{align*}\]

En considérant la fonction \(\phi = 1\) on en déduit notamment que \(\mathbf{P} \big(c U_1 g(Y_1) < f(Y_1) \big) = p = \frac{1}{c}\) et donc

\[\begin{equation*} \frac{\mathbf{E} \big[ \phi(Y_1) \mathbf{1}_{\{c U_1 g(Y_1) < f(Y_1)\}} \big]} {\mathbf{P} \big(c U_1 g(Y_1) < f(Y_1) \big)} = \frac{1}{c} \int_{\mathbf{R}^d} \phi(y) f(y) \mu(\operatorname{d}\! y) = \mathbf{E} [\phi(X) ], \end{equation*}\]

c’est à dire que la loi conditionnelle \(Y_1 | \{c U_1 g(Y_1) < f(Y_1) \}\) est la loi \(\nu\).

Remark 1.1

Si on pose par récurrence, pour tout \(n \ge 1\),

\[\begin{equation*} \begin{cases} \tau_1 = \min\big\{k \ge 1, \, c U_k g(Y_k) < f(Y_k) \big\} \\ \tau_{n + 1} = \min\big\{k \ge \tau_n + 1, \, c U_k g(Y_k) < f(Y_k) \big\} \\ \end{cases} \end{equation*}\]

alors la suite \(X_n = Y_{\tau_n}\) est une suite i.i.d. de loi \(\nu\).

1.5.1. Efficacité de la méthode du rejet

En pratique il est important de s’intéresser à l’efficacité de la méthode du rejet. La variable aléatoire \(\tau\) représente le nombre nécessaire de variables aléatoires uniformes pour simuler une variable aléatoire de loi \(\nu\). La variable \(\tau\) suit une loi géométrique de paramètre \(p = \frac{1}{c}\) où \(c > 1\) est le paramètre apparaissant dans l’hypothèse de domination (1.5). En moyenne le nombre de variables uniformes utilisées pour 1 réalisation de \(X\) est

\[\begin{equation*} \mathbf{E}[\tau] = \frac{1}{p} = c. \end{equation*}\]

On peut aussi s’intéresser au rendement moyen qui est défini par

\[\begin{equation*} \mathbf{E}\big[\frac{1}{\tau} \big] = \sum_{n \ge 1} \frac{1}{n} p(1-p)^{n-1} = - \frac{p}{1-p} \log(p). \end{equation*}\]

Plus \(c > 1\) est proche de 1 et plus la méthode est efficace. Ce paramètre est important et en pratique on cherche à construire une loi auxiliaire \(\bar \nu\) proche de la loi à simuler \(\nu\) (densité \(g\) proche de \(f\)).

1.5.2. Interprétation géométrique de la méthode du rejet

Il est possible de faire réinterpréter le théorème précédent et de voir le temps aléatoire \(\tau\) comme le premier instant où l’on tombe dans un domaine particulier: sous la densité \(f\) selon laquelle on souhaite simuler. On rappelle que l’hypographe d’une fonction \(h:\mathbf{R}^d \to \mathbf{R}_+\) est le domaine défini par

\[\begin{equation*} \mathcal{H}(f) = \big\{ (x, y) \in \mathbf{R}^d \times \mathbf{R}_+, \, 0 \le y \le f(x) \big\}. \end{equation*}\]

Proposition 1.4

Soit \(f: \mathbf{R}^d \to \mathbf{R}_+\) une densité de probabilité i.e. \(\int f(x) \operatorname{d}\! x = 1\).

1. Soit \(X\) de densité \(f\) et \(U \sim \mathcal{U}\big(]0,1[\big)\) indépendante de \(X\). Alors pour tout \(c >0\), le couple \((X, c U f(X))\) est uniforme sur \(\mathcal{H}(c f)\).

2. Si \((X, Y)\) est uniforme sur \(\mathcal{H}(f)\) alors \(X\) suit la loi de densité \(f\). \end{enumerate}

Proof. 1. Soit \(X\) de densité \(f\), \(U \sim \mathcal{U}\big(]0,1[\big)\) indépendante de \(X\), \(c > 0\) et une fonction test continue bornée \(g:\mathbf{R}^d \times \mathbf{R}_+ \to \mathbf{R}\). Alors

\[\begin{equation*} \mathbf{E}\big[ g(X, c U f(X)) \big] = \int_{\mathbf{R}^d} \int_0^1 g(x, c u f(x)) \operatorname{d}\! u f(x) \operatorname{d}\! x, \end{equation*}\]

et par le changement de variable \(y = u c f(x)\) on a

\[\begin{equation*} \mathbf{E} \big[ g(X, c U f(X)) \big] = \int_{\mathcal{H}(cf)} g(x, y) \frac{1}{c} \operatorname{d}\! x \operatorname{d}\! y. \end{equation*}\]

On conclut en remarquant que \(\lambda_{d+1}(\mathcal{H}(c f)) = c\). 2. Soit \((X, Y) \sim \mathcal{U}\big(\mathcal{H}(f) \big)\) et \(g:\mathbf{R}^d \to \mathbf{R}\) une fonction test continue bornée. Alors

\[\begin{align*} \mathbf{E}[g(X)] &= \int_{\mathcal{H}(f)} g(x) \operatorname{d}\! x \operatorname{d}\! y, \\ & = \int_{\mathbf{R}^d} \int_0^{f(x)} g(x) \operatorname{d}\! x \operatorname{d}\! y = \int_{\mathbf{R}^d} g(x) f(x) \operatorname{d}\! x, \end{align*}\]

c’est à dire \(X\) de loi de densité \(f\).

A l’aide de cette proposition, on peut relire l’algorithme du théorème précédent Theorem 1.2 de la façon suivante:

• on simule \((Y_n, c U_n g(Y_n))_{n \ge 1}\) réalisations i.i.d. uniforme sur \(\mathcal{H}(c g)\)

• on garde les réalisations qui tombent dans \(\mathcal{H}(f)\)