Algorithmes Stochastiques

J.Ph. Chancelier¹

Version pdf de ce document

1 Recherche de zéros

On cherche ici à trouver numériquement $u^{\star} \in {\mathbb{R}}^d$ tel que $\Phi(u^\star)= \alpha$ où $\Phi$ est une fonction strictement croissante donnée et $\alpha$ un réel donné. On suppose par ailleurs que la fonction $\Phi$ est donnée par :

$$\Phi(u) = {\mathbb{E}}\left[ f(u,\xi)\right]$$

où $\xi$ est une variable aléatoire que l'on sait simuler et $f$ une fonction que l'on sait évaluer. Un exemple typique est $f(u,\xi) = \Phi(u) + \xi$. Pour $u$ donné on ne connaît $\Phi(u)$ qu'a travers une mesure $f(u,\xi)$ qui est entachée d'une erreur de mesure $\xi$ de loi ${\cal N}(0,1)$.

Un algorithme stochastique pour trouver un zéro de $\Phi(u)- \alpha$ sera le suivant. On se donne $U_0$ et une suite de variables aléatoire $\xi_1,\ldots, \xi_n$ indépendantes et de même loi que $\xi$. On met à jour $U_n$ l'estimation de $u^{\star}$ par :

$$U_{n+1} = U_n - \gamma_n ( f(U_n,\xi_n) - \alpha)$$

où $\gamma_n$ est une suite déterministe décroissant vers 0 et telle que :

$$\sum \gamma_n = \infty$$   et$$\quad \sum \gamma_n^2 < \infty$$

La suite $U_n$ converge alors p.s (voir le cours) vers $u^{\star}$ si l'on suppose par exemple :

• $f(u,\xi)$ est bornée.
• $\Phi(u).(u- u^{\star}) \ge c {\Vert u - u^{\star} \Vert}^2$, $c > 0$

Écrire une fonction Scilab qui mets en oeuvre cet algorithme. La fonction dont on cherche à calculer le zéro sera passée en argument. Faire tourner l'algorithme sur une fonction de votre choix ($d=1$, ou $2$, ...).

2 Un calcul de quantiles

On regarde dans cette section le cas particulier où on cherche à évaluer un quantile. Soit $\Phi(u)$ la fonction de répartition d'une loi ${\cal L}$. On suppose ici $\Phi(u)$ continue. Calculer un quantile d'ordre $\alpha$ de ${\cal L}$ consiste à calculer $u_{\alpha}$ vérifiant :

$${\mathbb{P}}\left[ X \le u_{\alpha}\right] = \alpha$$

où $X$ est une variable aléatoire de loi ${\cal L}$. Soit $f(u,x) = {\mathbb{I}}_{\{ x \le u \}}$ cela s'écrit encore :

$${\mathbb{E}}\left[ f( u_{\alpha}, X)\right] = \alpha$$

On peut donc utiliser l'algorithme décrit dans la section précédente. une autre méthode est aussi envisageable. Soient $X_1,\ldots, X_n$ $n$ tirages indépendants de loi ${\cal L}$. On peut les réordonner en ordre croissant $Y_1, \ldots , Y_n$ et estimer le quantile d'ordre $\alpha$ par $Y_j$ avec

$$j/n \le \alpha \le (j+1)/n$$

Cela revient à calculer le quantile d'ordre $\alpha$ de la fonction de répartition empirique :

$$\Phi_n(u) = \frac{1}{n} \sum_{k=1}^n {\mathbb{I}}_{ Y_k \le u}$$

Mettre en oeuvre les deux algorithmes par exemple pour la loi normale ou pour d'autres loi de votre choix (grand). Pour l'algorithme sur la fonction de répartition empirique on pourra utiliser dsearch pour rajouter un élément dans un tableau déjà trié.

Faire une animation graphique montrant l'évolution des deux algorithme aux cours des itérations. Pour certaines lois donc la loi normale on peut calculer les quantiles dans Scilab (cdfnor). On pourra utiliser cela pour superpose sur les courbes la solution $u_{\alpha}$.

 
driver('X11') 
xset('pixmap',1) // on utilise un mémoire graphique pour dessiner 
for i=1:n 
   .......
   xset('wwpc')  // efface la mémoire graphique 
   // commande graphiques à introduire ici 
   xset('wshow') // affiche le contenu de la mémoire graphique à l'écran.
end

3 Recherche d'un maximum (Kiefer-Wolfowitz)

On suppose cette fois que l'on veut maximiser une fonction concave $\Phi(u) = {\mathbb{E}}\left[ f(u,\xi)\right]$. On pourrait utiliser des méthodes d'optimisation déjà vue, par exemple une méthode de gradient sur la fonction $\Phi{u}$ mais cela demanderait à chaque itération l'évaluation d'une espérance ${\mathbb{E}}\left[ \frac{\partial}{\partial_u} f(u,\xi)\right]$. On va donc plutôt utiliser un algorithme stochastique pour trouver un zéro de ${\mathbb{E}}\left[ \frac{\partial}{\partial_u} f(u,\xi)\right]$ et plutôt que de calculer explicitement le gradient on va l'évaluer par :

$$\frac{1}{2}(\Phi(u+c) - \Phi(u-c))$$

et $c$ va tendre vers 0 au cours de l'algorithme itératif. Soient $\xi_{1}^1,\xi_{1}^2,\ldots,\xi_n^1 , \xi_n^2$ une suite de variables i.i.d de même loi que $\xi$ la remise à jour de $U_n$ se fait par :

$$U_{n+1} = U_n + (\gamma_n/c_n) ( f(U_n+c_n,\xi^1_n) - f(U_n-c_n,\xi^2_n) )$$

où $\gamma_n$ est une suite déterministe décroissant vers 0 et telle que :

$$\sum \gamma_n = \infty \quad, \sum \gamma_n c_n < \infty$$   et$$\quad \sum (\gamma_n/c_n)^2 < \infty$$

$U_n$ converge alors p.s vers $u^{\star}$ pour $\Phi$ de classe ${\cal C}^2$ strictement concave et telle que $\vert\Phi''(u)\vert \le K(1+\vert u\vert)$ et ${\mathbb{E}}\left[f^2(u,\xi)\right] \le K(1+u^2)$.

Mettre en oeuvre cet algorithme dans Scilab.

Notes

... Chancelier¹

Cermics, École Nationale des Ponts et Chaussées, 6 et 8 avenue Blaise Pascal, 77455, Marne la Vallée, Cedex 2