Échantillonnage préférentiel

L'échantillonnage préférentiel, en anglais importance sampling, est une méthode de réduction de la variance qui peut être utilisée dans la méthode de Monte-Carlo. L'idée sous-jacente à l'échantillonnage préférentiel, EP dans la suite, est que certaines valeurs prises par une variable aléatoire dans une simulation ont plus d'effet que d'autres sur l'estimateur recherché. Si ces valeurs importantes se réalisent plus souvent, la variance de notre estimateur peut être réduite.

Par conséquent la méthode de l'EP est de choisir une loi qui « encourage » les valeurs importantes. L'utilisation d'une loi biaisée conduira à un estimateur biaisé si nous l'appliquons directement aux simulations. Cependant, les différentes simulations sont pondérées afin de corriger ce biais ; l'estimateur EP est alors sans biais. Le poids qui est donné à chaque simulation est le ratio de vraisemblance, qui est la densité de Radon-Nikodym de la vraie loi par rapport à la loi biaisée.

Le point fondamental dans l'implémentation d'une simulation utilisant l'EP est le choix de la loi biaisée. Choisir ou créer une bonne loi biaisée est l'art des EP. L'avantage peut alors être une énorme économie de temps de calculs alors que l'inconvénient pour une mauvaise loi peut être des calculs plus longs qu'une simple simulation de Monte-Carlo.

Théorie

En Monte-Carlo

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

$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

$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 $f_{X}(\cdot )={\frac {1}{b-a}}$ sa densité.

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

${\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 $\mathbb {E} ({\hat {g}}_{N})=G$) et consistant (d'après la loi des grands nombres). Sa variance est :

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

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

$\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}$

Principe de l'échantillonnage préférentiel

L'idée principale de l'échantillonnage préférentiel est de remplacer dans la simulation la densité uniforme sur $[a;b]$ par une densité alternative (ou densité biaisée), notée $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 $\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 :

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

ce qui revient à :

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

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

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

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

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

avec enfin

${\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 $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 à

$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é $f^{\ast }$ reproduisant (imitant) la fonction g.

En quasi Monte Carlo

Pour estimer l'intégrale $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

$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 :

$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 $F\,'(x)=f^{*}(x)>0$. Il apparaît clairement que si $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 $f^{*}$ est définie à un facteur K près qui disparaît dans le quotient. On peut donc imposer que $\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 :

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

Application : estimation d'une probabilité

Considérons que nous voulons estimer par simulation la probabilité p_t d'un événement ${X\geq t\ }$ où X est une variable aléatoire de densité de probabilité $f$. Ce problème se ramène à la présentation générale dans le sens où il met en œuvre une intégrale à estimer. Un échantillon $\left(X_{i}\right)_{i\in \{1,\dots ,K\}}$ identiquement et indépendamment distribué (i.i.d.) est tiré dans cette loi. On note k_t le nombre de réalisations supérieures ou égales à t lors de K essais. La variable k_t est une variable aléatoire suivant une loi binomiale de paramètres K et p_t :

$P(k_{t}=k)={K \choose k}p_{t}^{k}(1-p_{t})^{K-k},\,\quad \quad k=0,1,\dots ,K$

ce qui signifie notamment que $\mathbb {E} (k_{t})=Kp_{t}$: la fréquence empirique $k_{t}/K$ converge donc vers sa probabilité associée p_t (loi des grands nombres).

L'échantillonnage préférentiel entre en jeu ici pour diminuer la variance de l'estimation Monte-Carlo de la probabilité p_t. En effet, p_t est donnée par

${\begin{aligned}p_{t}&={\mathbb {E} }[I_{\{x\geqslant t\}}]\\&=\int I_{\{x\geqslant t\}}{\frac {f(x)}{f^{\ast }(x)}}f^{\ast }(x)\,dx\\&={\mathbb {E} _{*}}[I_{\{x\geqslant t\}}w(X)]\end{aligned}}$

où, on a encore posé

$w(\cdot )\equiv {\frac {f(\cdot )}{f^{\ast }(\cdot )}}$

La dernière égalité de l'équation précédente suggère l'estimateur de $p_{t}$ suivant :

${\hat {p}}_{t}={\frac {1}{K}}\,\sum _{i=1}^{K}I_{\{x\geqslant t\}}w(X_{i}),\,\quad \quad X_{i}\sim f^{\ast }$

C'est un estimateur EP de p_t qui est sans biais. Ceci étant défini, la procédure d'estimation est de générer un échantillon i.i.d. à partir de la densité $f^{\ast }\,$ et pour chaque réalisation dépassant t de calculer son poids W. Le résultat sera la moyenne obtenue avec K tirages. La variance de cet estimateur est :

${\begin{aligned}{\mbox{Var}}^{\ast }{\hat {p}}_{t}&={\frac {1}{K}}{\mbox{Var}}^{\ast }[I_{\{x\geqslant t\}}w(X)]\\&={\frac {1}{K}}\left[\mathbb {E} ^{*}[I_{\{x\geqslant t\}}w^{2}(X)]-p_{t}^{2}\right]\\&={\frac {1}{K}}\left[\mathbb {E} [I_{\{x\geqslant t\}}w(X)]-p_{t}^{2}\right]\end{aligned}}$

Là encore il faudra profiler au mieux la densité $f^{\ast }$ afin de diminuer la variance.

Exemples numériques

Intégration de la fonction bêta par l'usage d'une loi triangulaire

Détail de la méthode

On souhaite estimer la quantité suivante :

$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 $g(x)=x^{4}(1-x)^{2}$.

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

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

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

et l'estimateur de sa variance :

${\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 du graphe de la fonction bêta, on peut remplacer la loi uniforme standard par la loi triangulaire ${\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

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

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

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

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

et l'estimateur de sa variance est

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

Comparaison de la méthode Monte-Carlo et de l'échantillonnage préférentiel

	Monte-Carlo classique			Échantillonnage préférentiel
n	estimateur	biais	variance	estimateur	biais	variance
500	0,009843	-3,19E-004	1,32E-007	0,009712	-1,88E-004	2,50E-008
1000	0,009735	-2,12E-004	6,53E-008	0,009680	-1,57E-004	1,26E-008
2500	0,009628	-1,04E-004	2,60E-008	0,009576	-5,18E-005	5,36E-009
5000	0,009717	-1,93E-004	1,31E-008	0,009542	-1,83E-005	2,71E-009
10000	0,009634	-1,10E-004	6,52E-009	0,009544	-1,99E-005	1,35E-009

On espère améliorer encore les performances en considérant une densité $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.

Intégration d'une gaussienne

La difficulté d'intégration d'une telle fonction est que la variable d'intégration prend des valeurs sur $\mathbb {R}$. Dans ce cas, utiliser une densité de probabilité uniforme est impossible, car la probabilité d'occurrence pour $x\in \mathbb {R}$ serait $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 $-\infty$ à $+\infty$, lorsque la loi employée permet elle-même de fournir des valeurs dans $\mathbb {R}$.

Voir aussi

Liens internes

• 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)