4. Méthodes de simulation

4.1 Un plan d'expérience nécessite des "nombres au hasard"

Un plan d'expérience nécessite des "nombres au hasard". En effet :

1. L'indépendance est indispensable.
2. Le tirage au sort permet de garantir cette indépendance ... et de pouvoir la prouver !
3. Il est donc nécessaire de savoir "tirer au sort par ordinateur" selon des lois fixées par avance.

4.2 Générateur uniforme

Remark 4.2.1   L'objectif poursuivi est de créer un comportement aléatoire au sein d'un ordinateur, alors que cet ordinateur a tout au contraire été conçu pour être strictement déterministe. A cet effet, on utilise un dispositif fournissant des "nombres aléatoires". Un tel dispositif peut être matériel, de la pièce de pile ou face au compteur de désintégration radioactif en passant par une horloge mesurant le nombre de micro-secondes entre deux frappes au clavier.

Definition 4.2.2   On appelle "générateur" un dispositif logiciel utilisant un processus déterministe pour fournir des nombres ayant un comportement réputé être indiscernable du comportement d'un processus effectivement aléatoire.

Definition 4.2.3   On appelle "générateur aléatoire" (noté $\mathrm{ra}\left(\right)$ dans ce qui suit) un calcul déterministe fournissant des nombres appartenant à $\left[0,1\right]$ et censés ressembler aux nombres que fournirait une suite de variables aléatoires iid (indépendantes et identiquement distribuées) selon la loi uniforme.

Definition 4.2.4   Les générateurs censés imiter d'autres lois seront appelés : générateur discret (pour un nombre fini de valeurs possibles), générateur exponentiel, générateur normal, etc.

Remark 4.2.5   A nouveau : les "générateurs aléatoires" sont aussi peu aléatoires que possible. Cette désignation est donc tout à fait détestable, mais a l'avantage de la brièveté et de la consécration par l'usage. Nous dirons donc qu'un générateur est bon ou mauvais selon que les séquences obtenues possèdent ou non certaines propriétés.

Remark 4.2.6   Propriétés statistiques. Un certain nombre des propriétés requises sont de nature statistique. Ainsi sera-t-il exigé que la moyenne d'une séquence tant soit peu longue soit "raisonnablement proche" de l'espérance de la variable. Une propriété analogue est également exigée pour les variances : il suffirait sinon d'instancier $x$ par $\mathrm{E}\left(x\right)$ pour obtenir la bonne moyenne !

Remark 4.2.7   Propriété de période. Le générateur, prenant ses valeurs dans un ensemble fini et étant déterministe, est nécessairement périodique. Il est donc exigé que la période de ce générateur soit largement supérieure à la longueur de toutes les séquences d'instanciations que l'on pourra envisager d'utiliser.

Remark 4.2.8   Propriété de non-prévisibilité. Pour certaines applications (cryptographie, jeux de hasard, etc. ) il est indispensable que le prochain tirage ne soit pas prévisible à partir des tirages précédents. Pour d'autres (simulations sur ordinateur) cette propriété n'est pas utile.

Remark 4.2.9   Propriétés complémentaires. En premier lieu vient la rapidité d'exécution du générateur. En second lieu vient la qualité requise quant à la preuve des propriétés du générateur. Il existe en effet des générateurs "sans défauts connus", mais dont on ne connaît pas non plus de preuve théorique de leurs qualités !

Scilab 4.2.10   Le générateur "de base" de Scilab est rand(). Il fournit un nombre dans $\left[0,1\right]$. Les commandes rand('seed') et rand('seed', n) permettent respectivement, de lire et d'écrire l'état interne du générateur (seed).

Algorithm 4.2.11 (générateur congruentiel)   Le LISTING 1 donne le générateur utilisé par LGM (Park and Miller, 1988), le LISTING 2 donne celui utilisé par Maple et le LISTING 3 donne celui utilisé par Scilab. Dans les trois cas, la formule utilisée est :

$$seed\;\mapsto\;\left(a\times seed+b\right)\,\mod\, p$$

Les choix intéressants sont d'une part $p$ premier, $b=0$ et $a$ générateur de $\left(\mathbb{Z}/p\,\mathbb{Z}\right)^{*}$ conduisant à une période $p-1$ et d'autre part $p=2^{m}$, $b$ impair et $a$ d'ordre $p/4$ conduisant à une période $p/2$. Comme il y a "beaucoup" d'éléments d'ordre maximal vis à vis du module $p$, les $a$ ne sont pas difficiles à trouver. Pour une discussion plus précise, se reporter à (Douillet, 2004).

Algorithm 4.2.12 (générateur standard)   Pour générer des nombres décimaux dans l'intervalle $\left[0,\,1\right]$ à partir d'un générateur congruentiel, une division est nécessaire. Un choix habile du module permet d'éviter le coût d'une division en virgule flottante. Un choix possible est $p$ puissance de deux : pour des nombres écrits en binaire, il suffit de déplacer la virgule. Lorsque les nombres sont stockés en codage décimal, prendre $p$ proche d'une puissance de dix conduit au même résultat. Ainsi le générateur uniforme principal de Maple est LISTING 4.

4.3 Choix uniforme d'un individu parmi $n$

Algorithm 4.3.1   Pour obtenir des pseudo-instanciations $\left(x_{j}\right)$ d'une variable entière $\mathbf{x}$ uniformément répartie dans l'intervalle $\left[1\,..\, n\right]$, on utilise la partie entière par excès de $n\times\mathrm{ra}\left(\right)$, soit $x=\mathrm{ceil}\left(n\times\mathrm{ra}\left(\right)\right)$ (LISTING 5).

Remark 4.3.2   L'algorithme LISTING 5 est soumis aux limitations : soit $m>n^{2}$ soit $m>\max\left(n,\, j^{3}\right)$. Pour Maple, cela fait $n<10^{6}$ ou $j<10^{4}$.

4.4 Méthode de tir

Algorithm 4.4.1 (cas discret)   Soit à instancier une variable entière $X\in\left[1,\, n\right]$ proportionnellement à des chances $c\left(x\right)$ fixées d'avance. Posons $M\doteq\max\left\{ c\left(x\right)\mid1\leq x\leq n\right\}$.
Nous procédons en trois étapes : (a) tirage d'un candidat $x$, uniformément parmi les valeurs possibles pour la variable ; (b) tirage d'un verdict, uniformément dans l'intervalle $\left[0,\, M\right]$ ; (c) acceptation ou rejet : si $c\left(x\right)\leq M$ on accepte la valeur $x$ et, sinon, on recommence jusqu'à obtenir une valeur acceptable.

Remark 4.4.2   Lorsque le nombre $n$ est petit, et les chances normalisées pour que $M=1$, on peut accélérer le processus en n'appelant qu'une seule fois le générateur aléatoire (LISTING 6).

Algorithm 4.4.3 (support fini)   Cette méthode se généralise au cas d'une variable à support borné $X\in\left[a,\, b\right]$ distribuée selon une densité bornée $\forall x,\, c\left(x\right)\leq M$. L'algorithme équivaut à tirer uniformément un point dans le rectangle $\left[a,\, b\right]\times\left[0,\, M\right]$ et à rejeter les points situés au-dessus du graphe de la fonction $c$.

Definition 4.4.4   On appelle facteur de tir la proportion de réussite par essai. On a clairement :

$$\phi=\frac{\int_{a}^{b}c\left(x\right)\, \mathrm{d}x}{\int_{a}^{b}M\, \mathrm{d}x}$$

qui est le rapport entre la surface utile (située sous le graphe) et la surface totale du rectangle.

Remark 4.4.5   Il est clair qu'un choix $M>sup\left(c\right)$ serait également possible, mais aurait pour seul résultat de ralentir le processus. On supposera toujours que le diagramme de tir a été optimisé pour que $M=sup\left(c\right)$.

Proposition 4.4.6   Le nombre moyen d'essais par valeur obtenue vaut $1/\phi$.

Exercise 4.4.7   Retrouver ce résultat en remarquant que $Pr\left(x=k\right)=\phi\,\left(1-\phi\right)^{k-1}$.

FIG. 4.1: Un diagramme de tir.

Example 4.4.8   Considérons le diagramme de tir donné FIG. 4.1. Les chances correspondantes sont :

$$\left[0.8,\,0.7,\,0.5,\,0.4,\,0.1\right]$$

Utilisé "tel que", le facteur de tir est $0.5$. Il est clair qu'une normalisation des chances optimise le processus. En utilisant

$$\left[1,\,0.875,\,0.625,\,0.500,\,0.125\right]$$

le facteur de tir passe à $0.625$ : il faut en moyenne $1.6$ essais par nombre engendré.

4.5 Algorithme de Walker

Algorithm 4.5.1   On modifie la méthode de tir, de façon à obtenir un facteur de tir égal à 1 : lorsque le verdict est inférieur à $c\left(x\right)$, le résultat est $x$ (comme précédemment). Dans le cas contraire, le résultat est donné par une table auxiliaire (cf LISTING 7).

1. La variable $\mathbf{u}$ uniformément répartie dans $\left[0,1\right]$ nous sert en premier lieu à choisir uniformément une méthode de tir parmi les $n$ possibles, puis la "partie inutilisée" de $u$, c'est à dire la partie fractionnaire de $n\, u$, fournit le verdict utilisé par la méthode de tir sélectionnée. Toute l'affaire consiste à construire les deux tables $gap$ et $aux$, contenant l'une les seuils de rejet (nombres de $\left[0,1\right]$) et l'autre des entiers dans $\left[1..n\right]$ indiquant qui est choisi à la place du candidat par défaut lorsque celui-ci est rejeté.
2. Illustrons ce problème en considérant un exemple, avec $n=9$ et
   $${\displaystyle chances=\left[\frac{370}{480},\,\frac{983}{480},\,... ...\,\frac{153}{480},\,\frac{121}{480},\,\frac{335}{480},\,\frac{684}{480}\right]}$$
   ce vecteur ayant été normé pour que la somme de ses éléments soit égale à $n$. Nous allons trier les éléments de $chances$ par valeur croissante, tout en mémorisant l'indice d'origine, ce qui conduit à la liste :
   $${\displaystyle L=\left[\begin{array}{ccccccccc} \frac{121}{480} &... ...{480} & \frac{983}{480}\\ 7 & 5 & 4 & 8 & 3 & 1 & 9 & 6 & 2\end{array}\right]}$$
   A chaque étape de l'algorithme, nous allons supprimer le premier élément de $L$ et modifier le dernier (et donc le replacer à l'endroit nécessaire pour que $L$ reste triée par ordre croissant). Par construction, le premier élément est inférieur ou égal à 1, et le dernier est supérieur ou égal à 1. Nous posons donc $gap\left(7\right)=\frac{121}{480}$ et $aux\left(7\right)=2$, nous supprimons le couple $\left[\frac{121}{480},\,7\right]$ et nous remplaçons la valeur $\frac{983}{480}$ par $\frac{983}{480}-1+\frac{121}{480}>0$. Autrement dit, le couple $\left[\frac{983}{480},\,2\right]$ est remplacé par le couple $\left[\frac{624}{480},\,2\right]$, qui vient s'intercaler entre $\left[\frac{370}{480},\,1\right]$ et $\left[\frac{684}{480},\,9\right]$. Nous obtenons donc une nouvelle liste de $m=6$ nombres positifs et de somme $m$ : il est clair que le procédé peut être répété tant que $L$ n'est pas vide. On obtient le LISTING 8 et les listes :
   $$\left[\begin{array}{c} gap\\ aux\end{array}\right]=\left[\begin{... ...7}{96} & \frac{143}{160}\\ 9 & 2 & 2 & 6 & 6 & 2 & 2 & 9 & 6\end{array}\right]$$

   
   
   

3. Remarque. Le codage donné ci-dessus pour l'initialisation de l'algorithme de Walker n'a pas été optimisé outre mesure. En effet, c'est le nombre $N$ d'instanciations fournies par le générateur qui est susceptible de devenir important, bien plus que le nombre $n$ d'objets à choisir. Le fait d'utiliser une liste pour représenter $L$, et de la re-trier ( $\mathrm{sort}$) à chaque fois que l'on a supprimé le premier élément ( $\mathrm{subsop}\left(1=NULL,\,L\right)$) donne une complexité (approximativement) quadratique en $n$ : une meilleure implémentation peut se révéler utile lorsque $n$ augmente. En pareil cas, il peut également devenir nécessaire d'utiliser une nouvelle instanciation du générateur uniforme comme verdict dans l'algorithme principal.