Contenu et design par Dominique Lefebvre -
www.tangenteX.com octobre 2009 
Contact
: 
TangenteX.com
La physique numérique 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. 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:
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, elle est intéressante....
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 = (axi + b) mod c
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!
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 (xi, xi+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 !
J'ai implémenté cet algorithme en C, en développant dans
l'environnement Dev 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
//
// http://www.tangenteX.com">www.tangenteX.com
//**************************************************************************
#include <cstdkib>
#include <iostream>
#include <math.h>
#include <time.h>
using namespace std;
// déclaration des constantes
#define NBTIRAGE 10000 v // déclaration des variables globales
long seed;v
//**************************************************************
// 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"); v // fermeture du pipe
pclose(pipe);v return; v }
//***************************************************************
// 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); v 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 v //
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:
Pour télécharger le code source : RandomGenerator.cpp et le projet Dev C++ RandomGenerator . Le fichier de commande de GnuPlot est Random.p
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
//
// http://www.tangenteX.com">www.tangenteX.com
//**************************************************************************
#include <cstdkib>
#include <iostream>
#include <math.h>
#include <time.h>
using namespace std;
// déclaration des constantes
#define NBTIRAGE 10000 v // déclaration des variables globales
long seed;v
//***************************************************************************
// 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(); v 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
Contenu et design par Dominique Lefebvre -
www.tangenteX.com octobre 2009 Contact
: