[ edito | sommaire | contribution | contacts ]
Jean Debord
JDebord@compuserve.com http://ourworld.compuserve.com/homepages/JDebord/ |
| > I. Théorie |
| > I.A. Introduction |
Générer des nombres aléatoires sur ordinateur revient à créer une suite d'entiers :
où f est une fonction qui doit être choisie judicieusement pour que la répartition des nombres In ne puisse pas être distinguée de ce que donnerait le hasard. On parle alors de nombres pseudo-aléatoires.
La suite peut fournir M nombres aléatoires dans l'intervalle [0, M-1]. M dépend du type des entiers :
La valeur de départ I0, appelée graine (seed), doit être fournie par l'utilisateur. La même graine donnera toujours la même suite de nombres. D'autre part, la suite se reproduit au bout d'un nombre de valeurs appelé période. Cette période doit être la plus grande possible.
Si la fonction f a été correctement choisie, chaque nombre a la même probabilité d'apparaître. On dit que les nombres aléatoires suivent une loi uniforme.
On montre que les autres lois de probabilité (binomiale, normale...) peuvent être simulées à partir de cette loi uniforme.
Nous décrirons ici des méthodes permettant de générer des nombres aléatoires répartis selon :
une loi uniforme
une loi normale (courbe de Gauss)
une loi multinormale (loi normale à plusieurs variables)
| > I.B. Loi uniforme |
| > I.B.1. Générateurs congruentiels |
C'est la base des générateurs implantés dans les langages de programmation. Ils sont du type :
avec b = 216 par exemple.
Ces générateurs ont plusieurs inconvénients qui les font déconseiller dans la pratique.
| > I.B.2. Générateurs de type "Multiplication avec retenue" (Multiply With Carry, MWC) |
Ces générateurs, introduits par George Marsaglia, sont du type :
a est le multiplicateur, b la base, cn la retenue (on peut prendre c0 = 0). Pour b = 216, si a est choisi parmi les valeurs du tableau suivant, la période est (a * 215 - 1).
----------------------------------------------------------- 18000 18030 18273 18513 18879 19074 19098 19164 19215 19584 19599 19950 20088 20508 20544 20664 20814 20970 21153 21243 21423 21723 21954 22125 22188 22293 22860 22938 22965 22974 23109 23124 23163 23208 23508 23520 23553 23658 23865 24114 24219 24660 24699 24864 24948 25023 25308 25443 26004 26088 26154 26550 26679 26838 27183 27258 27753 27795 27810 27834 27960 28320 28380 28689 28710 28794 28854 28959 28980 29013 29379 29889 30135 30345 30459 30714 30903 30963 31059 31083 ----------------------------------------------------------- | l |
Un générateur 32 bits peut être obtenu en concaténant 2 nombres de 16 bits. La période est alors égale au produit des périodes des deux générateurs.
Ces générateurs ont de meilleures propriétés que les précédents et sont donc recommandés.
| >I.C. Loi normale |
| >I.C.1. Introduction |
La loi normale, notée N(m, s), est caractérisée par sa moyenne m et son écart-type s (ou sa variance s2). Elle est telle que :
La loi N(0, 1) est appelée loi normale réduite. Elle est représentée sur la courbe suivante (courbe de Gauss) :
Figure 1. Loi normale réduite.
| >I.C.2. Simulation d'une loi normale |
On utilise la méthode de Box-Muller : si x1 et x2 sont deux nombres aléatoires uniformes sur ]0,1[, les nombres y1 et y2 définis par :
suivent la loi normale réduite.
On en déduit que les nombres z1 = m + sy1 et z2 = m + sy2 suivent la loi normale N(m, s).
| >I.D. Loi multinormale |
| >I.D.1. Introduction |
La loi multinormale correspond à l'extension de la loi normale au cas de n variables x1, x2,... xn. Elle est caractérisée par un vecteur de moyennes m et une matrice de variance-covariance V. On la note N(m, V).
Le coefficient de corrélation rij s'en déduit par :
Ce coefficient est toujours compris entre -1 et 1. Il est positif si xi et xj ont tendance à varier dans le même sens, négatif dans le cas contraire.
| >I.D.2. Simulation d'une loi multinormale |
Pour simuler une loi multinormale N(m, V) de dimension n, on applique l'algorithme suivant :
Soit u un vecteur constitué de n nombres aléatoires indépendants distribués selon la loi normale réduite.| >I.E. Loi multi-lognormale |
| >I.E.1. Introduction |
On dit que le vecteur x suit la loi multi-lognormale LN(m, V), où m désigne le vecteur de moyennes et V la matrice de variance-covariance, si le vecteur x° tel que :
suit une loi multinormale.
| >I.E.2. Simulation d'une loi multi-lognormale |
Si x suit la loi multi-lognormale LN(m, V), de dimension n, on définit la variable auxiliaire x° distribuée selon la loi multinormale L(m°, V°) telle que :
Le vecteur x tel que xi = exp(x°i) suit alors la loi multi-lognormale LN(m, V).
| > II. Programmation en Turbo Pascal |
La programmation des calculs précédents peut être réalisée au moyen de la bibliothèque TP MATH.
| > II.A. L'unité RANDOM.PAS |
L'unité RANDOM.PAS, contenue dans TPMATH1.ZIP fournit les fonctions et procédures suivantes :
function IRanMar : LongInt;
Retourne un entier aléatoire de 32 bits obtenu en concaténant 2 entiers de 16 bits. Le premier bit est interprété par le Turbo Pascal comme étant le bit de signe ; les valeurs possibles sont donc :
soit :
Les multiplicateurs utilisés pour les 2 générateurs 16 bits sont a1 = 18000 et a2 = 30903. La période du générateur 32 bits est donc :
Nous avons vérifié que ce générateur passe les tests DIEHARD de G. Marsaglia.
procedure RMarIn(Seed1, Seed2 : Integer)
Initialise les 2 générateurs 16 bits utilisés par IRanMar. L'initialisation par défaut correspond à RMarIn(1802, 9373).
function RanMar : Float;
Retourne un réel aléatoire compris dans l'intervalle [0, 1[. Ce réel est calculé à partir d'un appel à IRanMar, selon :
function RanGauss(Mu, Sigma : Float) : Float;
Retourne un nombre aléatoire tiré de la loi normale de moyenne Mu et d'écart-type Sigma.
procedure RanMult(M : PVector; L : PMatrix;
N : Integer; X : PVector);Retourne un vecteur aléatoire X tiré de la loi multinormale à N dimensions ayant pour vecteur de moyennes M ; L est le facteur de Cholesky de la matrice de variance-covariance (la décomposition doit être faite par le programme appelant).
| > II.B. Programmes de démonstration |
Les programmes suivants sont disponibles dans TPMATH2.ZIP :
| > II.B.1. Test du générateur |
Le programme RANTEST.PAS a pour seul but de vérifier que le générateur de nombres aléatoires fonctionne avec votre compilateur. Pour cela, le programme tire 20000 nombres entiers, avec l'initialisation par défaut RMarIn(1802, 9373), affiche les 6 nombres suivants et les compare à leurs valeurs théoriques.
Les résultats obtenus doivent être les suivants :
Test of Marsaglia random number generator
---------------------------------------------
Correct Actual
---------------------------------------------
921625997 921625997 OK
1094293978 1094293978 OK
115775252 115775252 OK
499820504 499820504 OK
-1929018715 -1929018715 OK
2008943384 2008943384 OK
---------------------------------------------
| l |
| > II.B.2. Simulation d'une loi normale |
Le programme RANGAUS.PAS simule le tirage d'un échantillon aléatoire à partir d'une population gaussienne. Il estime la moyenne et la variance de cette population au moyen des formules classiques :
où n désigne la taille de l'échantillon.
Ces estimations sont faites par les fonctions Average et EstVar de l'unité STAT.PAS
Le programme calcule ensuite un intervalle de confiance à 95% pour la moyenne, selon la formule :
valable pour les "grands" échantillons (n > 30), pour lesquels la moyenne suit une distribution normale.
Cet intervalle de confiance a 95% de chances de recouvrir la moyenne de la population.
Le programme affiche :
Avec l'initialisation par défaut du générateur de nombres aléatoires, les résultats sont les suivants :
Population mean = 10.0000 Population SD = 2.0000 Sample size = 100 Sample mean = 10.0005 Sample SD = 1.9670 95% CI of mean = [ 9.6150 , 10.3860 ] | l |
| > II.B.3. Simulation d'une loi multinormale |
Le programme RANMUL.PAS simule une loi multinormale. Le vecteur de moyennes et la matrice de variance-covariance sont lus dans un fichier ASCII dont la structure est la suivante :
Le fichier par défaut (RANMUL.DAT) est un exemple avec N = 2.
La procédure ReadParam lit le fichier de données. La matrice de variance-covariance est recalculée à partir des écart-types et des coefficients de corrélation :
procedure ReadParam(FileName : String; var Name : String;
var N : Integer; var M : PVector;
var V : PMatrix);
var
F : Text; { Fichier de données }
I, J : Integer; { Variables de boucle }
S : PVector; { Ecart-types }
R : Float; { Coefficient de corrélation }
begin
Assign(F, FileName);
Reset(F);
Readln(F, Name);
Readln(F, N);
DimVector(M, N);
DimVector(S, N);
DimMatrix(V, N, N);
{ Lire les moyennes et écart-types.
Calculer les variances }
for I := 1 to N do
begin
Read(F, M^[I], S^[I]);
V^[I]^[I] := Sqr(S^[I]);
end;
{ Lire les coefficients de corrélation.
Calculer les covariances }
for I := 2 to N do
for J := 1 to Pred(I) do
begin
Read(F, R);
V^[I]^[J] := R * S^[I] * S^[J];
V^[J]^[I] := V^[I]^[J];
end;
Close(F);
DelVector(S, N);
end;
| l |
Le programme principal effectue la décomposition de Cholesky de la matrice de variance-covariance et appelle la procédure RanMult autant de fois que nécessaire. Les vecteurs aléatoires ainsi générés sont écrits dans le fichier de sortie, à raison d'un vecteur par ligne. Le fichier de sortie par défaut est RANMUL.OUT
{ Décomposer la matrice de variance-covariance }
Cholesky(V, 1, N, L);
{ NSIM est le nombre de simulations }
for I := 1 to NSIM do
begin
{ Tirer un vecteur aléatoire X }
RanMult(M, L, N, X);
{ Ecrire le vecteur dans le fichier de sortie }
end;
| l |
| > II.B.4. Simulation d'une loi multi-lognormale |
Le programme RANMULL.PAS simule une loi multi-lognormale. Les fichiers d'entrée et de sortie ont le même format que pour le programme précédent.
