Principes de simulation

Introduction
Tirage de nombres aléatoires
Confection d'échantillons
Simulation suivant une loi de probabilité

Tests

Introduction

Imaginons de calculer le nombre p en choisissant des points au hasard dans un carré de côté 2a. Dans ce carré est inscrit le cercle de rayon a. On compte le nombre total N de points et le nombre n de points situés dans le disque. Si l'on admet que tous les points du carré ont la même probabilité d'être choisis, la probabilité de choisir un point dans le cercle sera approchée par le rapport n/N.

Par ailleurs, la répartition des points étant supposée homogène, le nombre de points situés dans le carré est proportionnel à l'aire du carré 4a². De la même façon, le nombre de points situés dans le disque est proportionnel à l'aire du disque pa². On en déduit

Il suffit donc de faire un grand nombre d'expériences en utilisant, pour la commodité, un ordinateur qui tirera au hasard les coordonnées (x,y) d'un point ( l'origine des coordonnées étant prise au centre du cercle) en suivant la loi uniforme et en vérifiant si x² + y² <= a².

Ce petit exemple illustre le principe d'une simulation :

1) On ne fait pas une expérience réelle (cela prendrait beaucoup de temps), mais on utilise l'ordinateur pour effectuer les tirages (constitution d'échantillons aléatoires ) et comptabiliser les cas favorables et les cas possibles. Il s'agit donc d'une expérience "virtuelle".

2) On se rend vite compte que si on veut un résultat significatif, il faut un nombre de tirages assez élevé. C'est d'ailleurs une bonne raison d'utiliser un ordinateur.

Normalement, on s'attend à ce que le résultat converge vers la "bonne valeur". On peut alors se fixer un seuil de précision si l'on admet que le résultat converge ;  par exemple on peut arrêter le tirage lorsque la différence entre deux valeurs successives est inférieure à 0,001. Voici une nouvelle version de l'application précédente :

NB : la simulation utilise des tirages de points sur une surface donnée ; ce type de simulation est appelé "Monte-Carlo".

Autre simulation possible : toujours le calcul de p, mais avec la méthode de Buffon (1717). Il s'agit de jeter une aiguille de longueur L sur un parquet composé de lattes parallèles de même largeur d. On compte le nombre de fois N où l'on lance l'aiguille et le nombre de fois n où elle coupe une séparation de deux lattes. Buffon a montré que la probabilité de couper une séparation de lattes est (dans le cas où d > L) :

2L/pd

Autrement dit, pour N grand , le rapport n/N doit tendre vers cette valeur ce qui permet d'en déduire une valeur approchée de p :

Tirage de nombres aléatoires

Comme on le voit dans les exemples précédents, la simulation est basée sur l'utilisation de nombres aléatoires pour la fabrication d'échantillons artificiels.

Lorsque les ordinateurs n'existaient pas encore (ou étaient rares), il était possible (il est d'ailleurs toujours possible) d'utiliser des tables de nombres aléatoires obtenus par l'observation de phénomènes physiques obéissant au hasard. Une table bien connue de cette époque reculée était celle de la RAND Corporation qui avait été établie par observation du débit des tubes électroniques. L'utilisation de telles tables consistait à extraire une suite de nombres aléatoires, par exemple

58360937587312385305968851265431298745609854325789043268721095010285730265742892

On peut ensuite effectuer des groupements de deux chiffres par exemple pour obtenir des nombres au hasard entre 00 et 99 :

58 36 09 37 58 73 12 38 53 05 96 88 51 26 54 31 29 87 45 60 98 54 32 57 89 04 32 68 72 10 95 01 02 85 73 02 65 74 28 92

Évidemment, de nos jours, c'est l'ordinateur qui fournit directement ces nombres au hasard. Mais comme ils sont calculés, on ne peut pas dire qu'ils sont aléatoires, mais pseudo-aléatoires. Examinons quelques procédés de génération de nombres pseudo-aléatoires.

procédé du mid-square

On choisit un nombre initial u₀ de 2m chiffres. On l'élève au carré ce qui donne un nombre à 4m chiffres et on prend comme résultat le nombre u₁ formé des 2m chiffres centraux. On recommence le procédé sur u₁ ce qui conduit à u₂, etc...

exemple :   

ce qui donne la suite 73672726431057611891.....

procédé de Lehmer

On part de deux nombres au hasard u₀  de m chiffres et k de n chiffres. On les multiplie ce qui donne un nombre à m + n chiffres. On sépare les n chiffres de gauche ce qui constitue un nombre d₀ que l'on retranche nombre f₀ constitué des m chiffres de droite ; on obtient ainsi u₁. On réitère le processus avec u₁, etc...

exemple :    u₀ = 483215    k = 678

on obtient ainsi une suite 483215619443981935751265357161...

procédé de Fibonacci

On part avec trois nombres u₀, u₁ et k, ce dernier étant plus grand que les deux premiers. On génère une suite de nombres aléatoires avec l'algorithme suivant :

u_n+2 = u_n+1 + u_n + ak             avec a = 0 si u_n+1 + u_n <= k     et     a = -1 si u_n+1 + u_n > k

exemple :    u₀ = 05673      u₁ = 09251       k = 12654.     On obtient successivement

Ensuite, on élimine le premier et le dernier chiffre de chaque nombre obtenu, d'où la suite :

567925227152113000114114228......

Validité des procédés

Il existe bien d'autres procédés, mais il faut qu'ils soient acceptables, c'est à dire fournir des nombres équiprobables, ce que l'on peut vérifier en faisant, par exemple le test du Chi2.

Confection d'échantillons

Prenons un exemple pour être très concret. Imaginons un guichet devant lequel arrive des clients. On mesure tous les 15 minutes le nombre de clients qui arrivent et on obtient après 100 mesures les résultats statistiques suivants :

On peut considérer que cette statistique approche la loi de probabilité des arrivées supposée inconnue. La moyenne est l = 4,77 (en prenant comme unité le quart d'heure).

On mesure aussi, de la même manière le nombre de clients servis par quart d'heure. 100 mesures fournissent le tableau statistique suivant :

Comme précédemment, cette statistique approche la loi de probabilité des départs. On constate que la moyenne est m = 4,91

On souhaite étudier ce système d'attente et obtenir par simulation la longueur moyenne de la file d'attente des clients attendant d'être servis.

On commencera par compléter les deux tableaux ci-dessus en faisant correspondre à chaque fréquence une plage de nombres de 2 chiffres :

La signification de ces nombres est la suivante. Imaginons que l'on possède deux suites de nombres aléatoires (ou pseudo aléatoires fournis par un ordinateur) :

S1 : 90 27 14 39 52 29 24 79 .....

S2 : 72 71 67 53 43 97 30 98 60 ....

Pour le 1er quart d'heure, on tire 90 pour les arrivées et 72 pour les départs ; avec les tableaux cela signifie que 8 clients arrivent et que 6 clients partent pendant ce quart d'heure.
Pour le second quart d'heure, on tire 27 pour les arrivées et 71 pour les départs ; ces nombres correspondant à 4 arrivées et 6 départs.
On continue ainsi pour une journée entière supposée de 6 heures (24 périodes d'un quart d'heure) et on obtient le tableau suivant dans lequel figurent le nombre d'arrivées et le nombre théorique de services qui pourraient être rendus s'il y avait suffisamment de clients :

On constate que la file est en moyenne de 1,33 individus ce qui peut sembler assez satisfaisant. Un autre résultat qui peut être intéressant est le taux "d'oisiveté" du guichetier ; en effet, comme on le constate, à certaine période, il pourrait rendre plus de services qu'il n'y a de clients et on peut donc considérer qu'il est inactif une partie du temps. On constate que le guichetier a effectué 99 services et il aurait pu en effectuer 128, d'où un taux d'oisiveté de 33 % environ.

Une simulation de ce genre est à interpréter avec précaution

1) le nombre d'événements n'est pas très important et on a intérêt à effectuer plusieurs simulations (qui donneront des résultats différents) et à faire des moyenne sur l'ensemble des simulations.

Une moyenne sur 100 simulations a donné la valeur 7,74 comme moyenne de la longueur de la file d'attente (donc bien différente de la valeur trouvée à partir du tableau ci-dessus).

2) le système est fermé puisque le service ne s'effectue que sur une durée de 6 heures (supposées ici continues) ; il est bien évident qu'il existe des phénomènes de "bord" : il faut attendre le premier client pour le servir et la fermeture du service peut limiter le nombre des dernières entrées. Par suite, il est quelquefois convenable d'effectuer les calculs seulement en période de "croisière".

Simulation suivant une loi de probabilité

Il arrive que l'on connaisse les lois de probabilité régissant les phénomènes. Dans ce cas, il est possible, en utilisant l'ordinateur, de générer des nombres dont la distribution suit une loi de probabilité. Par exemple, supposons dans l'exemple ci-dessus que la loi des arrivées et celle des services soient des lois de Poisson :

la probabilité d'avoir n arrivées de clients est : avec l= 4,77

la probabilité d'avoir n services effectués est : avec m  = 4,91

Il faudra donc utiliser un programme générateur de nombres aléatoires suivant la distribution de Poisson. Rappelons que pour trouver un nombre aléatoire suivant la loi de Poisson, on détermine un nombre entier a tel que

où z est un nombre aléatoire compris entre 0 et 1 suivant la loi uniforme. On génère les nombres suivant la loi de Poisson de la manière suivante (en langage algorithmique). La fonction gpoiss() utilise les fonctions gunif0() et fac().

fonction int gpoiss(float m)
{
    int a = -1 ;
    float s = -1 ;
    float v = 0 ;
    float q = exp(-m) ;  -- m est la moyenne de la loi de Poisson
    float u = gunif0() ;   -- tirage d'un nombre compris entre 0 et 1 suivant la loi uniforme
    Tant que s <= u Faire
        a = a + 1 ;
        v = v + m^a/fac(a) ;  -- m^a signifie m à la puissance a et fac(a) est la factorielle de a 
        s = v*q ;
    FinTantque
    return a ;
}