Générateur de nombres aléatoires

La génération de nombres aléatoires

Le thème du site est axé autour de la physique numérique. Cette science fait une consommation immodérée de simulations en tous genres. Et parmi ces simulations, beaucoup traitent de phénomènes aléatoires: désintégration nucléaire, arrivée de clients dans une file d'attente, percolation, marche au hasard, etc. Il est donc nécessaire de disposer d'un outil qui permette à un programme de générer des nombres aléatoires. Pas simple lorsqu'on sait qu'un ordinateur et un programme sont des objets en principe déterministes!
Autant le dire tout de suite, il est impossible de créer avec un algorithme une suite de nombres parfaitement aléatoire. Je ne vais pas en développer ici les raisons et je vous renvoie à la littérature spécialisée. En fait, seuls des générateurs physiques peuvent générer du pur aléatoire, ceux basés sur les désintégrations nucléaires ou le bruit thermique par exemple. On se contentera donc de générer des suites pseudo-aléatoires avec un ordinateur.
Une suite de nombres pseudo-aléatoire est une suite périodique, de période aussi grande que possible, dont les nombres, à l'intérieur de chaque période, sont aussi peu corrélés entre eux que possible (i.e. le plus indépendants les uns de autres). Dans la suite du texte, je parlerai par paresse de "suite aléatoire": il faut comprendre bien sur "suite pseudo-aléatoire".
Deux physiciens, Park et Miller ("Random number generators: good ones are hard to find" in Communications of the ACM 31,1192-201) ont définis en 1988 les trois qualités d'un générateur de nombres aléatoires:

la suite de nombres produite doit être de période très grande. Classiquement, sur une CPU de 32 bits, la période choisie est de 2^31 - 1. Lorsqu'on dispose d'une CPU de 64 bits, la période choisie est de 2^63 - 1. Ce n'est pas infini mais c'est pas mal quand même....
la corrélation entre les nombres doit être la plus faible possible. Il existe un moyen pas très compliqué de calculer cette corrélation (voir par exemple le §1.2 de la troisième partie du livre "Physique numérique" de P.Depondt) mais à notre niveau, je préfère un moyen plus ludique: on plot les nombres aléatoires sur un plan xy et on cherche visuellement un motif quelconque qui traduirait un certain ordre dans la suite aléatoire. C'est ce que je fais dans le programme exemple ci-dessous.
le générateur doit être aussi rapide que possible. Dans certaines applications, comme l'intégration par la méthode de Monte-Carlo, les appels au générateur de nombres aléatoires sont massifs. Il ne s'agit donc pas que celui-ci mette une seconde à délivrer un nombre! Il faut dire que ce point est moins sensible aujourd'hui qu'en 1988! Avec une station quadruple core à 3 GHz, on a pas vraiment de problème...

Reste maintenant à trouver un algorithme susceptible de répondre à ces trois critères! Il en existe un certain nombre... Ceux que ça intéresse peuvent consulter par exemple "Numerical recipes" en C ou Fortran, au chapitre 7 dans les deux ouvrages.
Dans le cadre de cette expérimentation, j'utilise un algorithme de génération d'une distribution uniforme simple de mise en oeuvre et très efficace, que je décris au paragraphe suivant. Il s'agit de l'algorithme de Lehmer (ou schéma de congruence linéaire) qui date de 1951 et a été publié pour la première fois dans "Mathematical Models in large scale computing unit in Proceedings of the Second symposium on Large scale digital computing machinery . En passant, allez voir l'histoire de Dick Lehmer sur Wiki, elle est intéressante....

Distribution uniforme de nombres aléatoires et algorithme de Lehmer

Je vous renvoie à votre cours de statistique pour établir la définition d'une distribution uniforme de variables aléatoires. Je vais rappeler ici quelques unes de ses propriétés qui nous intéresseront pour tester notre algorithme.
Pourquoi avoir choisi une distribution uniforme? Parce que c'est la plus simple à analyser et à mettre en oeuvre, parce qu'à partir de celle-ci on peut construire d'autres distributions, mais aussi et surtout parce que c'est la plus courante dans les phénomènes physiques.
Nous cherchons donc un algorithme qui génère des nombres aléatoires selon une distribution uniforme dans le domaine [0,1[. Pourquoi ce domaine? Par simplicité et parce que c'est suffisant. Les propriétés de la distribution uniforme sur ce domaine m'intéressent particulièrement. Sa moyenne est égale à:

\begin{align} <x> = \int_{0}^{1}xdx \end{align}

et sa variance est:

\begin{align} \sigma^2 = \int_{0}^{1}x^2dx - <x>^2 \end{align}

En calculant ces valeurs, on obtient une moyenne de 0.5 et un écart type de 0.289. Ces valeurs sont caractéristiques de la distribution uniforme sur le domaine choisi.
Voyons maintenant le moyen d'obtenir la génération d'une suite pseudo-aléatoire de distribution uniforme. Comme dit plus haut, j'utilise le schéma de congruence linéaire qui est on ne peut plus simple, regardez :

x _i+1 = (ax_i + b) mod c

Une suite au niveau d'une Première S! Ceux qui voudront comprendre le fonctionnement de l'algo peuvent se plonger dans l'article de Lehmer ou dans "Computational physics" de Landau §10. Les autres me feront confiance.
Evidemment, vous aurez compris que le point sensible de l'algorithme tient dans le choix des paramètres a,b et c. Le paramètre c est le paramètre de normalisation, qui permet de restreindre les valeurs au domaine [0,1[. Il vaut donc 2^31 - 1 pour une machine 32 bits (ou 2^63 - 1 pour une 64 bits).
Le choix des paramètres a et b est plus large. J'ai choisi les valeurs les plus communes dans la littérature (celles de Park and Miller): a = 7^5 et b = 0. Vous pouvez modifier ces valeurs pour voir ce qui se passe, en matière de "randomness" sur la distribution uniforme!

Comment vérifier le caractère aléatoire de la distribution uniforme

J'ai choisi les moyens les plus simples! Je calcule la moyenne et l'écart type de la distribution pour un nombre de tirages donné et je les compare aux valeurs théoriques. Les plus courageux pourront améliorer le programme en calculant la corrélation entre les nombres générés et vérifier qu'elle est aussi petite que possible (voir la référence ci-dessus).
La méthode que je préfère, même si ce n'est pas la plus précise, c'est la méthode visuelle. A partir de la série générée, je trace chaque point sur un plan en choisissant pour coordonnées du point (x_i, x_i+10). Le i+10 est arbitraire, ça pourrait être i+20 ou i+50 sans problème! Je fais plusieurs tracés en faisant varier le nombre de tirages. Et je cherche dans les figures obtenues des motifs, des zones non homogènes ou des tâches particulières. Si la corrélation entre nombres est très faible, avec un nombre de tirages grand (10^4 ou 10^5), j'arrive à une couverture quasi-uniforme du domaine. L'algorithme de Lehmer est assez bon!

Le programme C et les résultats

J'ai implémenté cet algorithme en C, en développant dans l'environnement Dev C++, selon mes habitudes (voir la page d'initiation au C). J'utilise Gnuplot pour le tracé xy de la distribution.
Le code source du programme RandomGenerator:

//*****************************************************************************

// Générateur de nombres aléatoires

// Dominique Lefebvre - Octobre 2009 V1.0

// TangenteX

// tangenteX.com

//******************************************************************************

#include <cstdkib>

#include <iostream>

#include <math.h>

#include <time.h>

using namespace std;

// déclaration des constantes

#define NBTIRAGE 10000

// déclaration des variables globales

long seed;

//******************************************************************************

// Génération de nombres aléatoires distribués uniformément entre 0 et 1

// Selon la formule de Lehmer et l'algorithme et les paramètres de

// Parker and Miller (1988):

// x(i+1) = a*x(i) mod c

// Sa période est fixée par c, soit 2^31 - 1

//******************************************************************************

double Random(void)

{

const double a = 16807.0; // 7^5

>const double c = 2147483647.0; // 2^31 - 1

double x,y;

x = a*seed;

seed = (long)(fmod(x,c));

y = seed/c; // normalisation par c

return (y);

}

//******************************************************************************

// Génération de la constante d'initialisation du générateur Random

// à partir de la date courante selon Anderson (1990)

//******************************************************************************

long GenerationSeed(void)

{

struct tm *temps;

time_t maintenant;

int t0,t1,t2,t3,t4,t5;

long lseed;

// récupération de l'heure actuelle

time(&maintenant);

// enregistrement de l'heure actuelle

temps = localtime(&maintenant);

// conversion en structure tm

t0 = temps->tm_sec;

t1 = temps->tm_min;

t2 = temps->tm_hour;

t3 = temps->tm_mday;

t4 = temps->tm_mon;

t5 = temps->tm_year;

// calcul du seed d'après l'algo d'Anderson (seed varie entre 0 et 2^31-1

lseed = t5 + 70*(t4 + 12*(t3 + 31*(t2 + 23*(t1 + 59*t0))));

// s'assurer que seed est impair

if ((lseed%2) == 0) lseed = lseed-1;

return lseed;

}

//******************************************************************************

// Affichage des points aléatoires dans un carré [0..1][0..1] en utilisant

// GNUPlot

// Le fichier Random.p contient les commandes Gnuplot de tracé

// Les données sont contenues dans le fichier Donnees.dat

//******************************************************************************

void PlotSerie(void)

{

FILE *pipe;

// ouverture d'un pipe vers Gnuplot et lancement du tracé

pipe = popen("gnuplot -persist Random.p","w");

// fermeture du pipe

pclose(pipe);

return;

}

//******************************************************************************

// Programme principal

//******************************************************************************

int main(int argc, char *argv[])

{

FILE *Donnees;

char Buffer[20];

int i;

double Serie[NBTIRAGE],Somme, Moyenne, EcartType;

// initialisation du générateur

printf("Generateur de nombres aleatoires selon une distribution uniforme"

" - Dominique Lefebvre Octobre 2009\n\n");

seed = GenerationSeed();

// génération d'une série aléatoire de distribution uniforme

for(i=0;i<NBTIRAGE;i++) Serie[i] = Random();

// Vérification du caractère aléatoire de la série

// pour une distribution aléatoire uniforme entre 0 et 1, la moyenne de

// la distribution doit être de 0.5 et son écart 0,289 (1/12)

Somme = 0.;

for(i=0;i<NBTIRAGE;i++) Somme = Somme + Serie[i];

Moyenne = Somme/NBTIRAGE;

Somme = 0.;

for(i=0;i<NBTIRAGE;i++) Somme = Somme + (Serie[i] - Moyenne)*(Serie[i] - Moyenne);

EcartType = sqrt(Somme/NBTIRAGE);

printf("Moyenne : %.5f Ecart type: %.5f\n\n",Moyenne, EcartType);

printf("Pour une distribution aleatoire uniforme entre 0 et 1, "

"pour un nombre infini de tirages:\n");

printf("- la moyenne vaut 0.5\n");

printf("- l'ecart-type vaut 1/12 soit 0.289\n");

// stockage des données dans le fichier Donnees.dat pour affichage

// ce fichier est dans le même répertoire que le programme RandomGenerator

Donnees = fopen("Donnees.dat","w");

for(i=0;i<NBTIRAGE-10;i++) fprintf(Donnees,"%.5f %.5f\n",Serie[i],Serie[i+10]);

fclose(Donnees);

// Visualisation de la distribution aléatoire uniforme

PlotSerie();

// Fin du programme

return EXIT_SUCCESS;

}

Les points particuliers de ce programme:

le nombre de tirages, c'est à dire le nombre d'éléments de la série, est fixé par la constante NBTIRAGE (original, non!).
la routine GenerationSeed : elle génére le germe (seed) nécessaire à l'algo de Lehmer. La valeur de seed doit être différente avant chaque construction de distribution. J'obtiens cette garantie en construisant mon seed à partir de la date/heure de l'instant de lancement. Cette idée est due à S.L Anderson (1990), un camarade de SIAM (Society for Industrial and Applied Mathematics) dans son article "Random number generators on vector supercomputers and other advanced architectures" in "SIAM Review 32, 221-51.
La communication avec GnuPlot se fait par un pipe, technique que je n'avais pas encore présenté sur ce forum. Il faut bien sur que le programme wgnuplot.exe soit installé sur votre machine.
Rien d'autre de particulier, programme très simple qui doit tourner du premier coup.

Pour télécharger le code source : RandomGenerator.cpp et le projet Dev C++ RandomGenerator . Le fichier de commande de GnuPlot est Random.p

Un programme d'application curieux : l'estimation de Pi

Pour s'amuser, voici un petit programme, tiré du premier, qui permet de calculer une valeur estimée de pi à partir d'une série de nombres aléatoires. Je vous laisse réfléchir à la méthode utilisée, qui est très commune dans certaines techniques d'intégration (la célèbre méthode de Monte-Carlo par exemple!).

//*****************************************************************************

// Calcul du nombre Pi par tirages aléatoires

// Dominique Lefebvre - Octobre 2009 V1.0

// TangenteX

// tangenteX.com

//******************************************************************************

#include <cstdkib>

#include <iostream>

#include <math.h>

#include <time.h>

using namespace std;

// déclaration des constantes

#define NBTIRAGE 100000

// déclaration des variables globales

long seed;

//******************************************************************************

// Génération de nombres aléatoires distribués uniformément entre 0 et 1

// Selon la formule de Lehmer et l'algorithme et les paramètres de

// Parker and Miller (1988):

// x(i+1) = a*x(i) mod c

// Sa période est fixée par c, soit 2^31 - 1

//******************************************************************************

double Random(void)

{

const double a = 16807.0; // 7^5

const double c = 2147483647.0; // 2^31 - 1

double x,y;

x = a*seed;

seed = (long)(fmod(x,c));

y = seed/c; // normalisation par c

return (y);

}

//******************************************************************************

// Génération de la constante d'initialisation du générateur Random

// à partir de la date courante selon Anderson (1990)

//******************************************************************************

long GenerationSeed(void)

{

struct tm *temps;

time_t maintenant;

int t0,t1,t2,t3,t4,t5;

long lseed;

// récupération de l'heure actuelle

time(&maintenant);

// enregistrement de l'heure actuelle

temps = localtime(&maintenant);

// conversion en structure tm

t0 = temps->tm_sec;

t1 = temps->tm_min;

t2 = temps->tm_hour;

t3 = temps->tm_mday;

t4 = temps->tm_mon;

t5 = temps->tm_year;

// calcul du seed d'après l'algo d'Anderson (seed varie entre 0 et 2^31-1

lseed = t5 + 70*(t4 + 12*(t3 + 31*(t2 + 23*(t1 + 59*t0))));

// s'assurer que seed est impair

if ((lseed%2) == 0) lseed = lseed-1;

return lseed;

}

//******************************************************************************

// Programme principal

//******************************************************************************

int main(int argc, char *argv[])

{

int i, pi4;

double x,y,pi;

// initialisation du générateur

printf("Calcul du nombre PI par tirages aleatoires"

" - Dominique Lefebvre Octobre 2009\n\n");

seed = GenerationSeed();

// calcul de PI

pi4 = 0;

for(i=0;i<NBTIRAGE;i++)

{

x = Random();

y = Random();

if ((x*x + y*y) < 1) pi4++;

}

pi = 4.0*pi4/NBTIRAGE;

printf("Valeur estimee de PI: %.6f pour %i tirages\n",pi,NBTIRAGE);

// Fin du programme

return EXIT_SUCCESS;

}

Faites tourner ce programme! Vous verrez que si le nombre de tirages est grand, la valeur de pi obtenue n'est pas complètement idiote!
Le source: CalculPi.cpp