2.11. Echantillonnage d’importance (Importance sampling)

La méthode d’échantillonnage d’importance ou importance sampling (IS) est une méthode de réduction de variance très utilisée en pratique et notamment dans les problèmes difficiles d’événements rares. Le but est de changer la loi sous-jacente à simuler pour favoriser l’apparition de réalisations qui ont un impact sur l’estimateur et donc de réduire la variance.

Remark 2.4 (Préliminaires sur l’erreur relative)

On veut calculer \(p = \mathbf{P}(X \ge K)\) pour \(X\) une variable aléatoire réelle donnée et \(K\) fixé, \(K\) grand de sorte que \(p\) soit petit par exemple \(p \simeq 10^{-6}\). L’estimateur Monte Carlo standard s’écrit alors

\[\begin{equation*} I_n = \frac{1}{n} \sum_{k=1}^n \mathbf{1}_{X_k \le K}, \end{equation*}\]

et si \(p \simeq 10^{-6}\), il y a très peu de réalisations de \(X_k\) qui sont plus grandes que \(K\)… en moyenne 1 réalisation pour \(10^6\) générées. La variance asymptotique est donnée par

\[\begin{equation*} \operatorname{var}(\mathbf{1}_{X \ge K}) = p - p^2, \end{equation*}\]

et l’intervalle de confiance est de la forme

\[\begin{equation*} p \in \Big[I_n \pm \frac{1.96 \sqrt{p-p^2}}{\sqrt{n}} \Big]. \end{equation*}\]

L’erreur absolue définie par \(|p-I_n|\) est donc majorer par \(\frac{1.96 \sqrt{p-p^2}}{\sqrt{n}}\). Cette erreur en \(\sqrt{p} \simeq 10^{-3}\) est très grande devant la quantité à calculer \(p \simeq 10^{-6}\).

Pour avoir un autre point de vu on peut aussi raisonner avec l’erreur relative définie par \(\big|\frac{p-I_n}{p}\big|\) qui vérifie

\[\begin{equation*} \big|\frac{p-I_n}{p}\big| \le \frac{1.96 \sqrt{p-p^2}}{\sqrt{n} p} \simeq \frac{1.96}{\sqrt{n p}}. \end{equation*}\]

Pour avoir une erreur relative de 20% (assez importante) il faut alors choisir \(\sqrt{np} \sim 10\) i.e. \(n = 100 p^{-1}\). Pour \(p \simeq 10^{-6}\), cela revient à prendre \(n = 10^8\)…

Example 2.2 (Préliminaires sur le cas Gaussien)

Supposons \(X = \varphi(G)\) avec \(\varphi:\mathbf{R} \to \mathbf{R}\) croissante et \(G \sim \mathcal{N}(0,1)\). Alors on a

\[\begin{equation*} \mathbf{P}\big(\varphi(G) \ge K\big) = \int_{\mathbf{R}} \mathbf{1}_{\varphi(x) \ge K} f(x) \operatorname{d}\!x \quad \text{avec $f(x) = \frac{1}{\sqrt{2 \pi}} e^{-\frac{x^2}{2}}$}. \end{equation*}\]

Pour rendre l’événement \(\{X \ge K\}\) plus probable, l’idée est de modifier la loi sous-jacente, ici la loi de \(G\) une gaussienne centrée réduite. Par exemple si on translate de \(\theta\) cette variable aléatoire on va considérer des événements de la forme \(\{f(G+\theta) \ge K\}\). On note \(G^{(\theta)} = G+\theta\) qui est une variable aléatoire de loi \(\mathcal{N}(\theta, 1)\) et on note \(f_\theta(x) = \frac{1}{\sqrt{2 \pi}} e^{-\frac{(x- \theta)^2}{2}}\) sa densité de sorte que

\[\begin{equation*} \mathbf{P}\big(\varphi(G) \ge K \big) = \int \mathbf{1}_{\varphi(x) \ge K} \frac{f(x)}{f_\theta(x)} f_\theta(x) \operatorname{d}\!x, \end{equation*}\]

et en réécrivant cette intégrale sous forme probabiliste on a

(2.19)\[ \mathbf{P}\big(\varphi(G) \ge K\big) = \mathbf{E}\Big[ \mathbf{1}_{\varphi(G^{(\theta)}) \ge K} \frac{f(G^{(\theta)})}{f_\theta(G^{(\theta)})} \Big].\]

Le ratio des densités est ici explicite \(\frac{f(x)}{f_\theta(x)} = e^{- \theta x + \frac{\theta^2}{2}}\) et comme \(G^{(\theta)} = G + \theta\) on a \(\frac{f(G + \theta)}{f_\theta(G + \theta)} = e^{-\theta G - \frac{\theta^2}{2}}\). D’où, dans ce cas Gaussien,

\[\begin{equation*} \mathbf{P}\big(\varphi(G) \ge K\big) = \mathbf{E}\Big[\mathbf{1}_{\varphi(G + \theta) \ge K} e^{- \theta G -\frac{\theta^2}{2}}\Big]. \end{equation*}\]

L’écriture (2.19) correspond à un changement de probabilité qu’on peut réécrire de la façon suivante: la variable aléatoire \(G\) est de loi \(\mathcal{N}(0,1)\) sous la probabilité \(\mathbf{P}\) et de loi \(\mathcal{N}(\theta, 1)\) sous la probabilité \(\mathbf{P}_\theta\) absolument continue (même équivalente) par rapport à \(\mathbf{P}\) de densité de Radon-Nikodym \(\frac{\operatorname{d}\!\mathbf{P}_\theta}{\operatorname{d}\!\mathbf{P}} = e^{\theta G - \frac{\theta^2}{2}}\). En effet, à l’aide de cette probabilité \(\mathbf{P}_\theta\) et de la notation \(\mathbf{E}_{\theta}[X] = \int X(\omega) \mathbf{P}_\theta(\operatorname{d}\!\omega)\), l’égalité (2.19) s’écrit

(2.20)\[\begin{equation} \mathbf{P} \big(\varphi(G) \ge K \big) = \mathbf{E}_{\theta} \Big[ \mathbf{1}_{\varphi(G) \ge K} e^{-\theta G + \frac{\theta^2}{2}} \Big]. \end{equation}\]

Soit \((\Omega, \mathcal{A}, \mathbf{P})\) un espace de probabilité et \(X\) une variable aléatoire réelle de loi \(\mu\) (sous \(\mathbf{P}\)) et de carré intégrable (sous \(\mathbf{P}\)) i.e. \(\mathbf{E}[X^2] < +\infty\). On veut calculer \(I = \mathbf{E}[X] = \int X(\omega) \mathbf{P}(\operatorname{d}\!\omega)\). On considère une probabilité auxiliaire \(\tilde{\mathbf{P}}\) telle que \(\mathbf{P}\) est absolument continue par rapport à \(\tilde{\mathbf{P}}\) de densité de Radon-Nikodym \(L\) c’est à dire \(L\) est une variable aléatoire positive telle que

(2.21)\[\begin{equation} \frac{\operatorname{d}\!\mathbf{P}}{\operatorname{d}\!\tilde{\mathbf{P}}} = L \quad \text{ou encore} \quad \mathbf{P}(\operatorname{d}\!\omega) = L(\omega) \tilde{\mathbf{P}}(\operatorname{d}\!\omega). \end{equation}\]

Pour que \(\tilde{\mathbf{P}}\) soit bien une probabilité il faut que la variable aléatoire \(L\) vérifie \(\mathbf{E}[L^{-1}] = 1\). En effet

\[\begin{equation*} \tilde{\mathbf{P}}(\Omega) = \int \tilde{\mathbf{P}}(\operatorname{d}\!\omega) = \int L^{-1}(\omega) \mathbf{P}(\operatorname{d}\!\omega) = \mathbf{E}[L^{-1}]. \end{equation*}\]

Alors en notant \(\tilde{\mathbf{E}}\) l’espérance sous la probabilité \(\tilde{\mathbf{P}}\) on a

\[\begin{equation*} I = \mathbf{E}[X] = \tilde{\mathbf{E}}[X L] = \int X(\omega) L(\omega) \tilde{\mathbf{P}}(\operatorname{d}\!\omega). \end{equation*}\]

D’autre part on suppose que \(L\) est variable aléatoire telle que \(\mathbf{E}[X^2 L] < +\infty\) de sorte que \(XL\) soit de carré intégrable sous \(\tilde{\mathbf{P}}\), en effet

\[\begin{equation*} \tilde{\mathbf{E}}[X^2 L^2] = \mathbf{E}[X^2 L] < +\infty. \end{equation*}\]

On suppose de plus qu’on sait simuler la variable aléatoire \(XL\) qui suit une loi \(\tilde \mu\) sous \(\tilde{\mathbf{P}}\) (loi à déterminer par le calcul, on verra des exemples plus loin). L’estimateur Monte Carlo basé sur cette nouvelle représentation s’écrit alors

(2.22)\[\begin{equation} \tilde I_n = \frac{1}{n} \sum_{k=1}^n X_k L_k \quad \text{avec $(X_k L_k)_{k \ge 1}$ suite i.i.d. de loi $\tilde \mu$.} \end{equation}\]

Example 2.3

Supposons \(X = g(Y_1, \dots, Y_d)\) avec \(Y_1, \dots, Y_d\) indépendantes de même loi de densité \(f\) (par rapport à la mesure de Lebesgue \(\lambda\)).

On veut calculer \(I = \mathbf{E}[X] = \mathbf{E}[\phi(Y_1, \dots, Y_d)]\). Soit \(\tilde f\) une densité auxiliaire de même support que \(f\). On considère la variable aléatoire \(L = \prod_{i=1}^d \frac{f(Y_i)}{\tilde f(Y_i)}\) et la probabilité \(\tilde{\mathbf{P}}\) telle que \(\frac{\operatorname{d}\!\mathbf{P}}{\operatorname{d}\!\tilde{\mathbf{P}}} = L\) alors

\[\begin{equation*} \mathbf{E}\big[ g(Y_1, \dots, Y_d) \big] = \tilde{\mathbf{E}} \Big[g(Y_1, \dots, Y_d) \prod_{i=1}^d \frac{f(Y_i)}{\tilde f(Y_i)} \Big], \end{equation*}\]

avec \((Y_1,\dots, Y_d)\) de loi de densité \(\prod \tilde f(y_k)\) sous \(\tilde{\mathbf{P}}\).

Proposition 2.8

Soit \(X\) une variable aléatoire réelle de carré intégrable \(\mathbf{E}[X^2] < +\infty\) et telle que \(\mathbf{P}(X = 0) = 0\). Soit \(\mathbf{P}^*\) la probabilité définie par \(\frac{\operatorname{d}\!\mathbf{P}^*}{\operatorname{d}\!\mathbf{P}} = \frac{|X|}{\mathbf{E}[|X|]}\) i.e. \(\mathbf{P}\) est absolument continue par rapport à \(\mathbf{P}^*\) de densité de Radon-Nikodym \(L^* = \frac{\mathbf{E}[|X|]}{|X|}\).

1. Alors l’estimateur Monte Carlo basé sur la représentation \(I = \mathbf{E}_{\mathbf{P}^*}[X L^*]\) a la plus petite variance asymptotique possible (parmi les changements de probabilité possibles).

2. Si \(X \ge 0\) \(\mathbf{P}\)–\(p.s.\) la variance est nulle.

Proof. La variance asymptotique de l’estimateur \(I^*_n = \frac{1}{n} \sum_{k=1}^n X_k L^*_k\) est donnée par \(\operatorname{var}_{\mathbf{P}^*}(X L^*) = \mathbf{E}_{\mathbf{P}^*}[(X L^*)^2] - \mathbf{E}[X]^2\).

• Par définition de \(L^*\) on a \(\mathbf{E}_{\mathbf{P}^*}[(X L^*)^2] = \mathbf{E}[X^2 L^*] = \mathbf{E}[|X|]^2\) et donc

(2.23)\[ \operatorname{var}_{\mathbf{P}^*}(XL^*) = \mathbf{E}[|X|]^2 - \mathbf{E}[X]^2\]

Pour toute probabilité \(\tilde{\mathbf{P}}\) telle que \(\operatorname{d}\!\mathbf{P} = \tilde{L} \operatorname{d}\!\tilde{\mathbf{P}}\) avec \(\tilde{L}\) densité de Radon-Nikodym on a par l’inégalité de Jensen

\[\begin{equation*} \mathbf{E}[|X|]^2 = \tilde{\mathbf{E}}[|X| L]^2 \le \tilde{\mathbf{E}}[(X L)^2], \end{equation*}\]

d’où \(\operatorname{var}_{\mathbf{P}^*}(XL^*) \le \operatorname{var}_{\tilde{\mathbf{P}}}(X \tilde{L})\) et l’optimalité de \(\mathbf{P}^*\).

• D’après (2.23) si \(X \ge 0\) \(\mathbf{P}\)–\(p.s.\) alors \(\operatorname{var}_{\mathbf{P}^*}(X L^*) = 0\).

Remark 2.5

Il est impossible en pratique d’utiliser ce changement de probabilité optimal car \(\mathbf{P}^*\) dépend de \(\mathbf{E}[|X|]\) quantité que l’on ne connait pas. Par contre on peut chercher une probabilité auxiliaire \(\tilde{\mathbf{P}}\) proche de \(\mathbf{P}^*\).

Remark 2.6

Si on relit cette proposition dans le cas à densité \(X = g(Y)\) avec \(g:\mathbf{R} \to \mathbf{R}_+\) et \(Y\) de densité \(f\) (exemple précédent avec \(d=1\) pour simplifier), cela consiste à considérer une densité auxilisaire optimale \(f^*\) qui vérifie

\[\begin{equation*} L^* = \frac{f(Y)}{f^*(Y)} = \frac{\mathbf{E}[g(Y)]}{g(Y)}. \end{equation*}\]

Ainsi la densité optimale \(f^*\) doit être proportionnelle au produit \(g f\) i.e.

\[\begin{equation*} \forall y \in \mathbf{R}, \quad f^*(y) = \frac{1}{\mathbf{E}[g(Y)]} f(y) g(y). \end{equation*}\]

2.11.1. Changement de mesure exponentielle

Souvent on fait de l’échantillonnage d’importance en considérant une classe de probabilité indexée par un paramètre \(\theta\) puis on cherche la meilleure probabilité dans cette classe, c’est à dire le meilleur paramètre \(\theta\). Une classe de probabilité très utilisée est celle appelée changement de mesure exponentielle (ou exponential tilting en anglais, aka transformée d’Esscher en actuariat / finance).

Soit \((\Omega, \mathcal{A}, \mathbf{P})\) un espace de probabilité et \(X\) un vecteur aléatoires à valeurs dans \(\mathbf{R}^d\) de fonction de répartition \(F(x_1, \dots, x_d) = \mathbf{P}\big(X_1 \le x_1, \dots, X_d \le x_d\big)\). On note \(\psi:\mathbf{R}^d \to \mathbf{R}\) la log-Laplace (ou fonction génératrice des cumulants) la fonction définie par

\[\begin{equation*} \forall \theta = (\theta_1, \dots, \theta_d) \in \mathbf{R}^d, \quad \psi(\theta) = \log\Big(\mathbf{E} \big[e^{\langle \theta, X \rangle}\big] \Big), \end{equation*}\]

avec \(\langle \theta, x \rangle = \sum_{j=1}^d \theta_j x_j\). Soit \(\Theta = \{\theta \in \mathbf{R}^d, \psi(\theta) < +\infty\}\) ensemble non vide (car il contient \(0\)). Alors pour \(\theta \in \Theta\) on considère la probabilité auxiliaire \(\mathbf{P}_\theta\)

\[\begin{equation*} \frac{\operatorname{d}\!\mathbf{P}}{\operatorname{d}\!\mathbf{P}_\theta} = L_\theta \quad \text{avec} \quad L_\theta = e^{- \langle \theta, X \rangle + \psi(\theta)}. \end{equation*}\]

La quantité déterminisite \(e^{\psi(\theta)}\) est la constante de normalisation qui fait de \(\mathbf{P}_\theta\) une probabilité. En effet on vérifie que \(\mathbf{P}_\theta(\Omega) = \mathbf{E}[L_\theta^{-1}] = \mathbf{E}[e^{\langle \theta, X \rangle - \psi(\theta)}] = 1\).

La loi de \(X\) sous \(\mathbf{P}_\theta\) est donnée par sa fonction de répartition \(F_\theta\) qui vérifie donc

\[\begin{equation*} \begin{split} \forall x\in \mathbf{R}^d, \quad F_\theta(x) &= \mathbf{E} \big[ \mathbf{1}_{X_1 \le x_1, \dots, X_d \le x_d} L_\theta^{-1} \big], \\ &= \mathbf{E} \big[ \mathbf{1}_{X \le x} e^{\langle \theta, X \rangle - \psi(\theta)} \big]. \end{split} \end{equation*}\]

En particulier si \(X\) a une loi à densité \(f\) sous \(\mathbf{P}\), alors \(X\) a une loi à densité \(f_\theta\) sous \(\mathbf{P}_\theta\) avec

\[\begin{equation*} \forall x \in \mathbf{R}^d, \quad f_\theta(x) = e^{\langle \theta, x \rangle - \psi(\theta)} f(x). \end{equation*}\]

Example 2.4 (Formule de Cameron-Martin (dimension finie))

Comme application de ce cadre général de changement de mesure exponentielle, on revient au cas Gaussien dans \(\mathbf{R}^d\). Soit \(X = (X_1, \dots, X_d)\) un vecteur gaussien centré réduit sous \(\mathbf{P}\), \(X \sim \mathcal{N}(0, \operatorname{Id}_d)\) (composantes indépendantes de variance 1). Alors la log–Laplace est \(\psi(\theta) = \frac{|\theta|^2}{2}\) avec \(|\theta| = \sqrt{\sum_{j=1}^d \theta_j^2}\). Ainsi \(\Theta = \mathbf{R}^d\) et pour tout \(\theta \in \mathbf{R}^d\) on considère la probabilité auxiliaire \(\mathbf{P}_\theta\) définie par

\[\begin{equation*} \operatorname{d}\!\mathbf{P}_\theta = e^{\langle \theta, X \rangle - \psi(\theta)} \operatorname{d}\!\mathbf{P}. \end{equation*}\]

Alors la loi de \(X\) sous \(\mathbf{P}_\theta\) est un vecteur gaussien centré en \(\theta\) de composantes indépendantes de variance 1 i.e. \(X \sim \mathcal{N}(\theta, \operatorname{Id}_d)\) (vérifier ce résultat par le calcul et faire le lien avec l’exemple préliminaire).