Échantillonnage préférentiel

On souhaite estimer une quantité G, qui s'exprime sous la forme d'une intégrale :

${\displaystyle G=\int _{a}^{b}g(x)\,{\mbox{d}}x}$

On considère ici une intégration en dimension 1, mais on peut généraliser à une dimension quelconque.

Le principe de base des méthodes de Monte-Carlo est de voir l'intégrale précédente comme

${\displaystyle G=(b-a)\int _{a}^{b}g(x)f_{X}(x)\,{\mbox{d}}x=(b-a)\,\mathbb {E} (g(X))}$

où X est une variable aléatoire uniformément distribuée sur [a;b] et ${\displaystyle f_{X}(\cdot )={\frac {1}{b-a}}}$ sa densité.

Si on dispose d'un échantillon ${\displaystyle (x_{1},x_{2},\cdots ,x_{N})}$, indépendant et identiquement distribué (i.i.d.) selon ${\displaystyle {\mathcal {U}}([a;b])}$, on peut estimer G par :

${\displaystyle {\hat {g}}_{N}={\frac {(b-a)}{N}}\sum _{i=1}^{N}g(x_{i})}$

Il s'agit d'un estimateur de G non-biaisé (c'est-à-dire que ${\displaystyle \mathbb {E} ({\hat {g}}_{N})=G}$) et consistant (d'après la loi des grands nombres). Sa variance est :

${\displaystyle \sigma _{{\hat {g}}_{N}}^{2}={\frac {(b-a)^{2}\sigma _{g}^{2}}{N}}}$

avec ${\displaystyle \sigma _{g}^{2}}$ la variance de la variable aléatoire ${\displaystyle g(X)}$

${\displaystyle \sigma _{g}^{2}={\frac {1}{(b-a)}}\int _{a}^{b}g^{2}(x)\,{\mbox{d}}x-\left({\frac {1}{b-a}}\int _{a}^{b}g(x)\,{\mbox{d}}x\right)^{2}}$

L'idée principale de l'échantillonnage préférentiel est de remplacer dans la simulation la densité uniforme sur ${\displaystyle [a;b]}$ par une densité alternative (ou densité biaisée), notée ${\displaystyle f^{\ast }\,}$, qui tente d'imiter la fonction g. Ce faisant, on remplace les tirages uniformes, qui n'avantagent aucune région, par des tirages plus « fidèles ». Ainsi, l'échantillonnage est fait suivant l'importance de la fonction g : il est inutile de tirer dans les régions où g prend des valeurs non-significatives, pour, au contraire, se concentrer sur les régions de haute importance. On espère ainsi diminuer la variance ${\displaystyle \sigma _{g}^{2}}$. Autrement dit, si on se fixe un niveau d'erreur donné, l’échantillonnage préférentiel permet de diminuer théoriquement le nombre de simulations N par rapport à une méthode de Monte-Carlo classique.

L'intégrale à estimer est réécrite comme :

${\displaystyle G=\int _{a}^{b}{\frac {g(x)}{f^{\ast }(x)}}f^{\ast }(x)\,{\mbox{d}}x}$

ce qui revient à :

${\displaystyle G=\mathbb {E} ^{\ast }[w(X)]}$

où on a posé ${\displaystyle w(x)=g(x)/f^{\ast }(x)}$ (appelé ratio de vraisemblance) et où X est simulé selon la densité ${\displaystyle f^{\ast }}$. Il est facile de généraliser les résultats précédents : l'estimateur de G est

${\displaystyle {\tilde {g}}_{N}={\frac {1}{N}}\sum _{i=1}^{N}w(x_{i})}$

où ${\displaystyle (x_{1},x_{2},\cdots ,x_{N})}$ est un échantillon i.i.d. selon la densité ${\displaystyle f^{\ast }}$. La variance de l'estimateur est donnée par

${\displaystyle {\mbox{Var}}^{\ast }({\tilde {g}}_{N})={\frac {{\mbox{Var}}^{\ast }[w(X)]}{N}}}$

avec enfin

${\displaystyle {\mbox{Var}}^{\ast }[w(X)]={\mbox{Var}}^{\ast }\left[{\frac {g(X)}{f^{\ast }(X)}}\right]=\int _{a}^{b}\left[{\frac {g(x)}{f^{\ast }(x)}}\right]^{2}f^{\ast }(x)\,{\mbox{d}}x-G^{2}}$

Dès lors, le problème est de se concentrer sur l'obtention d'une densité biaisée ${\displaystyle f^{\ast }\,}$ telle que la variance de l'estimateur EP soit moindre que celle de la méthode de Monte-Carlo classique. La densité minimisant la variance (qui la rend nulle sous certaines conditions) est appelée densité biaisée optimale. Cette dernière est égale à

${\displaystyle f^{\ast }(x)={\frac {|g(x)|}{\displaystyle \int _{a}^{b}|g(x)|\,{\mbox{d}}x}}}$

mais ce choix est inopérant, car on recherche justement le dénominateur. Toutefois, on peut s'attendre à réduire la variance en choisissant une densité ${\displaystyle f^{\ast }}$ reproduisant (imitant) la fonction g.

Pour estimer l'intégrale ${\displaystyle G=\int _{a}^{b}g(x)\,{\mbox{d}}x}$ on peut également se passer de tout le formalisme probabiliste précédent. Au lieu d'utiliser des variables aléatoires, on se sert de suites à faible discrépance (suites de Sobol par exemple). En dimension 1 l'approche la plus simple est

${\displaystyle G=\int _{a}^{b}g(x)\,{\mbox{d}}x\quad \simeq \quad {\frac {b-a}{N}}\sum _{i=1}^{N}g\left(a+(b-a){\frac {i}{N}}\right)}$

De même qu'en Monte Carlo usuel, cette approximation de l'intégrale converge d'autant plus vite que la fonction g est proche d'une constante. Si g est rigoureusement constante il suffit de prendre N = 1 pour avoir l'intégrale exacte. Réduire la variance de g est par conséquent crucial ici aussi ; dans ce but, l'échantillonnage préférentiel s'utilise comme suit :

${\displaystyle G=\int _{a}^{b}g(x)\,{\mbox{d}}x=\int _{a}^{b}{\frac {g(x)}{f^{*}(x)}}f^{*}(x)\,{\mbox{d}}x\quad =\quad \int _{F(a)}^{F(b)}{\frac {g(F^{-1}(y))}{f^{*}(F^{-1}(y))}}\,{\mbox{d}}y}$

où l'on a fait le changement de variable y = F(x) avec ${\displaystyle F\,'(x)=f^{*}(x)>0}$. Il apparaît clairement que si ${\displaystyle f^{*}\simeq k.g}$ alors la fonction à intégrer à droite est proche d'une constante (et donc de faible variance).

Pour faire le lien avec l'interprétation probabiliste de la section précédente, on remarque que ${\displaystyle f^{*}}$ est définie à un facteur K près qui disparaît dans le quotient. On peut donc imposer que ${\displaystyle \int _{a}^{b}f^{*}(x)\,{\mbox{d}}x=1}$, ce qui en fait une densité de probabilité sur [a, b]. Le changement de variable s'interprète alors naturellement comme un changement de probabilité et on a la simplification :

${\displaystyle \int _{F(a)}^{F(b)}{\frac {g(F^{-1}(y))}{f^{*}(F^{-1}(y))}}\,{\mbox{d}}y\quad =\quad \int _{0}^{1}{\frac {g(F^{-1}(y))}{f^{*}(F^{-1}(y))}}\,{\mbox{d}}y}$

Cette technique se généralise immédiatement en dimension quelconque.

Détail de la méthode

On souhaite estimer la quantité suivante :

${\displaystyle G=\int _{0}^{1}x^{4}(1-x)^{2}\,{\mbox{d}}x}$

qui se trouve être la fonction bêta de paramètres (5;3), qui vaut G = 1/105 = 0,0095238095238095. Cela correspond au cas général avec a=0, b=1 et ${\displaystyle g(x)=x^{4}(1-x)^{2}}$.

«Densité» de la loi beta ${\displaystyle g(x)=x^{4}(1-x)^{2}}$

On simule un échantillon ${\displaystyle (y_{1},\cdots ,y_{n})}$ selon la loi uniforme standard ${\displaystyle {\mathcal {U}}[0;1]}$ pour obtenir l'estimateur Monte-Carlo classique :

${\displaystyle {\hat {g}}_{1}={\frac {1}{n}}\sum _{i=1}^{n}g(y_{i})}$

et l'estimateur de sa variance :

${\displaystyle {\hat {\sigma }}_{g_{1}}^{2}={\frac {1}{n}}\left[{\frac {1}{n}}\sum _{i=1}^{n}g^{2}(y_{i})-{\hat {g}}_{1}^{2}\right]}$

S'inspirant de la forme générale de la fonction bêta, on peut remplacer la loi uniforme standard par la loi triangulaire ${\displaystyle {\mathcal {T}}(a=0,c=2/3,b=1)}$.

Elle ressemble à un triangle basé sur le segment [0;1] et « culminant » en (2/3;2). Sa densité est

${\displaystyle f^{\ast }(x)={\begin{cases}3x&x\in [0;2/3]\\6(1-x)&x\in [2/3;1]\end{cases}}}$

Loi de densité ${\displaystyle f^{\ast }(x)}$ pour l'échantillonnage préférentiel de ${\displaystyle g(x)}$

On simule un échantillon ${\displaystyle (z_{1},z_{2},\cdots ,z_{n})}$ dans cette loi, par la méthode de la transformée inverse, et, en posant ${\displaystyle w(x)=g(x)/f^{\ast }(x)}$, l'estimateur EP est donné par

${\displaystyle {\hat {g}}_{2}={\frac {1}{n}}\sum _{i=1}^{n}w(z_{i})}$

et l'estimateur de sa variance est

${\displaystyle {\hat {\sigma }}_{g_{2}}^{2}={\frac {1}{n}}\left[{\frac {1}{n}}\sum _{i=1}^{n}w^{2}(z_{i})-{\hat {g}}_{2}^{2}\right]}$

Dans le tableau, on constate que l'utilisation de l'EP permet systématiquement de réduire la variance de l'estimation par rapport à l'estimation Monte-Carlo classique de même taille (c'est-à-dire à n donné). On constate aussi que la variance d'estimation est proportionnelle à 1/n : en passant de n = 1000 à n = 10 000 (multiplication par 10 de la taille), on réduit d'un facteur 10 la variance.

On espère améliorer encore les performances en considérant une densité ${\displaystyle f^{\ast }}$ « plus proche » de la densité f (la fonction g renormalisée). Le principal problème sera d'obtenir des simulations. Dans les cas les plus simples, comme la loi triangulaire, la méthode de la transformée inverse pourra suffire ; dans les cas plus complexes, il faudra avoir recours à la méthode de rejet.

La difficulté d'intégration d'une telle fonction est que la variable d'intégration prend des valeurs sur ${\displaystyle \mathbb {R} }$. Dans ce cas, utiliser une densité de probabilité uniforme est impossible, car la probabilité d'occurrence pour ${\displaystyle x\in \mathbb {R} }$ serait ${\displaystyle p(x)=\lim _{l\to \infty }{\frac {x}{l}}=0}$. Or, dans une intégration Monte-Carlo brute, il n'y a pas de connaissance a priori de la forme de la fonction, et toutes les valeurs de x doivent avoir une probabilité égale. Ainsi, l'échantillonnage préférentiel est une méthode permettant d'intégrer des fonctions pour une variable d'intégration comprise entre ${\displaystyle -\infty }$ à ${\displaystyle +\infty }$, lorsque la loi employée permet elle-même de fournir des valeurs dans ${\displaystyle \mathbb {R} }$.

• Méthode de Monte-Carlo
• (en) J. Morio et M. Balesdent, Estimation of Rare Event Probabilities in Complex Aerospace and Other Systems : a practical approach, Cambridge, Elsevier Science, 2015, 216 p. (ISBN 978-0-08-100091-5)