test des générateurs de nombres pseudo-aléatoires

[leon1789]

Bonjour à tous, surtout aux gens spécialistes des probabilités et de tests des générateurs de nombres pseudo-aléatoires

A la fin, je déclare "douteux" un générateur de nombres pseudo aléatoires. Ai-je raison, ai-je une justification correcte ? c'est ma question...

Venons en aux faits. Voici un code C que j'explique ci-dessous :

#include <stdio.h>
#include <math.h>
 
int main()
{
  int face=0 ;
  int fois=0 ;
 
// srand (time (NULL)) ; // si on veut indexer la graine sur le temps machine
 
  for (face=0; face<16; face++)
  {
    float var=0.;

    for (fois=0; fois<1000; fois++)
    {
// on réalise 1000 fois l'expérience suivante :
// on lance 16000 fois un dé à 16 faces et
// on compte le nombre de fois où on obtient "face"
      int cpt=0 ;
      int i=0 ;
      for (i=0; i<16000; i++) if ( rand()%16 == face) cpt++ ;
 
      float moy=1000.;
      var +=(cpt-moy)*(cpt-moy);
    }

    float emq=sqrt(var/999.) ;
 
    printf("face %i, ecart-type = %0.3f\n",face,emq) ;
  }
  return 0;
}

Rapide explication du programme :
- On demande 16000 tirages aléatoires d'un nombre entre 0 et 15, suivant la loi uniforme par hypothèse.
- On compte le nombre de fois où le nombre 0 est apparu (variable cpt) : l'espérance est évidemment 16000 / 16 = 1000.
- On calcule la somme des écarts quadratiques entre la variable cpt et l'espérance sur 1000 tirages (variable var).
- On recommence 15 fois les trois lignes précédentes avec les nombres 1,2,...,15 (variable face).
- On estime alors l'écart-type de la loi sous-jacente (variable emq) ;

De manière théorique, il est clair que la variable cpt suit la loi binomiale B( 16000 , 1/16 ) donc la valeur de la variable emq doit tourner autour de l'écart-type de cette loi binomiale, à savoir avec p=1/16, autrement dit 30.6 environ.

On teste en machine en lançant le programme, on obtient ce résultat :

face 0, ecart-type = 29.575
face 1, ecart-type = 32.129
face 2, ecart-type = 29.689
face 3, ecart-type = 31.405
face 4, ecart-type = 29.712
face 5, ecart-type = 32.305
face 6, ecart-type = 29.750
face 7, ecart-type = 31.380
face 8, ecart-type = 29.748
face 9, ecart-type = 32.076
face 10, ecart-type = 30.025
face 11, ecart-type = 31.255
face 12, ecart-type = 29.619
face 13, ecart-type = 32.197
face 14, ecart-type = 29.621
face 15, ecart-type = 31.225

Les écart-types tournent autour de 30.6 comme prévu. Mais la dispersion des écart-types m'interroge... Je poursuis donc avec une étude mathématique pour déterminer si la dispersion des écart-types est douteuse ou pas.
Comme je vous le dis depuis le début, la valeur théorique de l'écart-type que l'on calcul est \( \sqrt{16000.p.(1-p)} \)~ 30.6 . De manière équivalente, la valeur théorique de la variance est


Mais quelle est la dispersion théorique des mesures autour de cette valeur de la variance ? On peut très bien répondre à cette question. Je ne vais pas détailler les calculs, disons simplement qu'on utilise la loi du khi² pour prouver ce que je vais écrire ( confer https://fr.wikipedia.org/wiki/Loi_du_%CF%87%C2%B2 ). Il en résulte que l'écart-type de la variance est
( confer https://fr.wikipedia.org/wiki/%C3%89...pes_empiriques )

Ainsi, on peut dire que notre variance suit une loi normale de moyenne 937.5 et d'écart-type 42

Par ailleurs, il est d'usage de dire que les mesures supérieures à espérance + 2 écart-types sont peu probables (probabilité de 0.025 , confer https://fr.wikipedia.org/wiki/%C3%89...on_diagram.svg )
Ainsi, les mesures de notre variance ont toutes les chances de varier entre 937.5+2x42 et 937.5+2x42, autrement dit entre 853 et 1022.

Enfin, avec un coup de racine carrée, les mesures de notre écart-type ont beaucoup de chance de varier en dessous de

Or quand on regarde les résultats de la machine, cette valeur 32 est dépassée 4 fois, soit 25% du temps ... Comparer cela aux 2.5% : c'est 10 fois plus !
J'en conclus que le générateur de nombres pseudo aléatoires est douteux.
Etes-vous d'accord ?

Merci de m'avoir lu jusqu'ici...
-----

[Matmat]

Mais il faut tenir compte que c'est sur un échantillon très faible (16) que tu obtiens ton 25% .

[Médiat]

Bonjour,

Connaissez-vous : https://www.random.org/analysis/Analysis2001.pdf

[fregoli]

Rien de plus compliqué, pour une machine, donc déterministe par définition, de générer des séries de nombres aléatoires, donc une phénomène indéterministe par définition, sauf à y ajouter des dispositifs qui vont tenir compte de phénomènes extérieurs imprévisibles par la machine (fluctuation de courant, bruits divers (au sens perturbation externe d'un système), ...).
Qu'un programme informatique ne fournisse pas des séries parfaitement aléatoires ne me surprend pas outre mesure.

[A voir en vidéo sur Futura]

[leon1789]

Merci Médiat pour le document : j'y jette un oeil dès que je l'aurais téléchargé (prob en cet instant).

Mais il faut tenir compte que c'est sur un échantillon très faible (16) que tu obtiens ton 25% .

Je me suis aussi fait la remarque, l'échantillon est relativement petit certes, mais la fréquence observée est déjà énorme. Petit calcul facile :
Proba de dépassement sur un essai : 2.5 % = 1/40 (sauf erreur de ma part)
Proba d'au moins 4 dépassements sur 16 essais : avec la loi binomiale B(16, 1/40) , on a P(X >= 4 ) ~ 1/1790
J'ai pas de chance alors ?

[leon1789]

(...)Qu'un programme informatique ne fournisse pas des séries parfaitement aléatoires ne me surprend pas outre mesure.

ok, mais là, c'est quand même rapidement imparfait, non ?

[Matmat]

Merci Médiat pour le document : j'y jette un oeil dès que je l'aurais téléchargé (prob en cet instant).


Je me suis aussi fait la remarque, l'échantillon est relativement petit certes, mais la fréquence observée est déjà énorme. Petit calcul facile :
Proba de dépassement sur un essai : 2.5 % = 1/40 (sauf erreur de ma part)
Proba d'au moins 4 dépassements sur 16 essais : avec la loi binomiale B(16, 1/40) , on a P(X >= 4 ) ~ 1/1790
J'ai pas de chance alors ?

pour moi 25% n'est pas étonnant car c'est comme si tu évaluais la répartition pour une variable de Bernoulli d'éspérance 1/40 dont l'écart type vaut racine de rac(39)/40 = 0.156 en fait

[gg0]

Bonjour Leon1789.

Dans ton programme, tu ne calcules pas l'écart type, mais la racine carrée de la somme des carrés des écarts à un nombre, 1000. Donc une valeur systématiquement supérieure à l'écart type. Pour estimer l'écart type, il vaudrait mieux utiliser l'écart type d'échantillon (s) qui a l'avantage d'être sans biais et convergent. mais ça complique le programme !

Cordialement.

[minushabens]

il vaudrait mieux utiliser l'écart type d'échantillon (s) qui a l'avantage d'être sans biais et convergent.

asymptotiquement sans biais je pense.

[leon1789]

pour moi 25% n'est pas étonnant car c'est comme si tu évaluais la répartition pour une variable de Bernoulli d'éspérance 1/40 dont l'écart type vaut racine de rac(39)/40 = 0.156 en fait

oui, j'évalue une variable de Bernoulli (p=1/40) mais 16 fois, pas qu'une seule fois. C'est pas négligeable.
Il faudrait poursuivre le test pour constater que le pourcentage de dépassement reste très élevé. Cela éviterait ta contestation visant davantage de prudence.

[leon1789]

Bonjour Leon1789.

Dans ton programme, tu ne calcules pas l'écart type, mais la racine carrée de la somme des carrés des écarts à un nombre, 1000. Donc une valeur systématiquement supérieure à l'écart type. Pour estimer l'écart type, il vaudrait mieux utiliser l'écart type d'échantillon (s) qui a l'avantage d'être sans biais et convergent. mais ça complique le programme !

Cordialement.

Salut gg0,
je ne comprends pas ta remarque : je calcule emq := sqrt(var/999.) ; , c'est bien l'estimateur "sans biais" (ou presque) de l'écart-type (999 = 1000-1),
var étant la somme des carrés des écarts à l'espérance.
D'ailleurs, je viens de me rendre compte qu'il faut diviser par 1000 ... mais bon, ça ne change pas grand chose je pense.

[gg0]

Après réflexion, 1000 étant la moyenne "vraie", tu utilises bien un estimateur de l'écart type à condition de diviser par 1000. C'est l'estimateur lorsque la moyenne est connue.

je suis allé un peu vite !

Cordialement.

[leon1789]

ok trop vite gg0,
moi aussi quand j'ai programmé, j'ai divisé trop vite par 1000-1 par habitude, alors que c'est effectivement par 1000 qu'il faut le faire ici. Mais ça ne change quasi rien en fait.

Plus important, au niveau du raisonnement pour établir les fluctuations de l'écart-type , est-ce que ça colle ? Je pense que oui, mais je ne suis pas probabiliste et l'erreur est humaine.
J'ai vu un certain nombre de tests qui confirment, mais deux ou trois qui infirment, les résultats.

[leon1789]

Mais il faut tenir compte que c'est sur un échantillon très faible (16) que tu obtiens ton 25% .

Je me suis aussi fait la remarque, l'échantillon est relativement petit certes, mais la fréquence observée est déjà énorme. Petit calcul facile :
Proba de dépassement sur un essai : 2.5 % = 1/40 (sauf erreur de ma part)
Proba d'au moins 4 dépassements sur 16 essais : avec la loi binomiale B(16, 1/40) , on a P(X >= 4 ) ~ 1/1790

Avec quasiment la même proba, c'est comme si on jetait 12 fois une pièce de monnaie et que l'on tombait du même coté à chaque fois... Tu ne penses pas que cette pièce serait truquée ? En tout cas, elle serait vite interdite dans les casinos

[mike.p]

Salut,

avec un petit port en javascript ( facile et dispo pour tous sans compilateur ) , j'obtiens des valeurs comprises entre 29,325 et 31,628

 Cliquez pour afficher

face 0, ecart-type = 31,628
face 1, ecart-type = 30,275
face 2, ecart-type = 30,480
face 3, ecart-type = 30,610
face 4, ecart-type = 31,354
face 5, ecart-type = 30,557
face 6, ecart-type = 31,389
face 7, ecart-type = 31,424
face 8, ecart-type = 31,499
face 9, ecart-type = 30,570
face 10, ecart-type = 31,249
face 11, ecart-type = 29,515
face 12, ecart-type = 31,435
face 13, ecart-type = 31,392
face 14, ecart-type = 30,350
face 15, ecart-type = 29,325

adaptation JS utilisée :

 Cliquez pour afficher

[code]
function mainLeon()
{
var face , fois , moy=1000.;

// srand (time (NULL)) ; // si on veut indexer la graine sur le temps machine
//// inutile en JS

for (face=0; face<16; face++)
{
var data = 0.0;

for(fois=0; fois<1000; fois++)
{
// on réalise 1000 fois l'expérience suivante :
// on lance 16000 fois un dé à 16 faces et
// on compte le nombre de fois où on obtient "face"
var cpt=0 ;
for (var i=0; i<16000; i++)
{
if ( (Math.floor(Math.random()*16)% 16 ) == face)
cpt++ ;
}
data +=(cpt-moy)*(cpt-moy);
}

var emq=Math.sqrt(data/moy) ;
console.log("face %i, ecart-type = %0.3f\n",face,emq) ;
// document.write("face "+face+", ecart-type = "+emq+"<br/>") ;
}
}
[code]

Pour firefox. Pour lire les résultats activez la console ( Developpement, Console Web).
Avec d'autres navigateurs, il faut remplacer le console.log() par autre chose, au pire, par la ligne commentée qui suit.

[gg0]

En fait,

tout dépend du simulateur de hasard utilisé par le logiciel. Je me souviens que sur Scilab, certains moyens étaient préconisés comme moins "réguliers". Et j'ai vu une simulation de marche au hasard sur Amstrad tourner au déplacement systématique vers le coins haut droit de l'écran (même en changeant de germe).

Léon, as-tu refait l'expérience (refaire fonctionner un programme, ce n'est généralement pas trop coûteux) ?

Cordialement.

[Matmat]

Il doit y avoir une erreur de raisonnement car si j'ai bien compris l'écart type pour cpt est 30.6 donc si ton cpt est presque tout le temps ( sauf 2.5% du temps ) entre 1000 +-(30.6*2) alors ça va .... or ta face 5 ( c'est la pire ) qui a un emq=32.305 correspond à un cpt qui vaut en moyenne 1042 et c'est bien inférieur à 1061.2

[leon1789]

tout dépend du simulateur de hasard utilisé par le logiciel.

Je suppose utiliser un générateur qui suit la loi uniforme. C'est le plus commun, mais il y a des logiciels comme Scilab qui proposent des alternatives, en effet.


(Dans mon raisonnement, loi uniforme, puis loi binomiale approchée par loi normale, puis khi² approché par loi normale .)

Léon, as-tu refait l'expérience (refaire fonctionner un programme, ce n'est généralement pas trop coûteux) ?

Pour tout vous dire, ce n'est pas mon ordi qui a donné ces résultats : c'est une autre personne qui me les a donné. Du coup, je ne peux pas refaire les tests a priori.

Quand mon ordi tourne (en C ou en Maple), ça donne des tests conformes, un peu comme ceux que mike.p vient de poster en javascript.

[leon1789]

Il doit y avoir une erreur de raisonnement car si j'ai bien compris l'écart type pour cpt est 30.6 donc si ton cpt est presque tout le temps ( sauf 2.5% du temps ) entre 1000 +-(30.6*2) alors ça va .... or ta face 5 ( c'est la pire ) qui a un emq=32.305 correspond à un cpt qui vaut en moyenne 1042 et c'est bien inférieur à 1061.2

je ne comprends pas comment tu obtiens la moyenne de cpt_1, ... cpt_1000 connaissant l'écart-type de cette série.

Mais je pense comprendre que tu me soumets l'idée d'un cpt constant à 1032 (ce qui est dans les clous 1000+-2*30.6), on obtient un emq de 32 que je déclare hors des clous.
Je vais y réfléchir... c'est une question de proba

[leon1789]

Il doit y avoir une erreur de raisonnement car si j'ai bien compris l'écart type pour cpt est 30.6 donc si ton cpt est presque tout le temps ( sauf 2.5% du temps ) entre 1000 +-(30.6*2) alors ça va .... or ta face 5 ( c'est la pire ) qui a un emq=32.305 correspond à un cpt qui vaut en moyenne 1042 et c'est bien inférieur à 1061.2

Une moyenne cpt de 1042 ou de 1032, c'est très très très peu probable (donc je doute de ton contre-raisonnement) : en effet, la variable cpt suit la binomiale B(16000, 1/16) d'espérance 1000 et d'écart-type 30.6, donc la moyenne sur 1000 variables iid suit une loi très proche de la loi normale (TCL) d'espérance 1000 et d'écart-type . Cela implique que la moyenne de 1000 cpt se situe quasi-surement entre entre 998 et 1004 (j'ai pris 4 écart-types, c'est vraiment très sûr !)

[Matmat]

C'est peut être moi qui fait un raisonnement faux , je me suis basé sur le programme qui lorsque face==5 a renvoyé emq=32.305 , donc le cumul des var a été 32.305² que j'ai multiplié par 999 pour en déduire 1042569 de faces 5 sur 16000000 possibles

[leon1789]

C'est peut être moi qui fait un raisonnement faux , je me suis basé sur le programme qui lorsque face==5 a renvoyé emq=32.305 , donc le cumul des var a été 32.305² que j'ai multiplié par 999 pour en déduire 1042569 de faces 5 sur 16000000 possibles

Je ne comprends pas ce qui est souligné.

Cela étant, juste un exemple : 999 cpt valant 999 et 1 cpt valant 2020 donnent une moyenne de 1000 et un emq de 32.3 .
Donc un emq de 32.3 n'implique pas nécessairement une moyenne de 1042 (ou 1032) .

[Matmat]

tu as raison ce que je dis est faux ( j'avais mal compris le programme , il ne peut pas être passé 1042569 de fois sur face 5 puisqu'il y passe exactement 1000000 de fois et qu'il cherche les coincidence avec i ) ... excuse moi .

[gg0]

Leon1789 :

Les générateurs aléatoires sont divers, ils changent suivant les logiciels. Donc ce qu'on obtient avec l'un ne se reproduira probablement pas avec un autre. Sans parler du germe qui, s'il n'est pas changé, fera se reproduire les mêmes séquences quand on recommence.
C'est d'ailleurs pourquoi on utilise parfois des générateurs aléatoires (bruit de fond d'un système électronique, tirages de pile ou face, etc.) externes.

Si tes calculs sont corrects (je n'ai pas eu le temps de vérifier, je ne connais pas la loi de la variance d'un résultat d'une multinomiale), soit le générateur est biaisé, soit le tirage est un peu exceptionnel.

Cordialement.

[leon1789]

gg0,
les générateurs différent suivant les logiciels, je suis d'accord évidemment. Mais la loi qu'ils sont sensés suivre est toujours (au moins dans leur paramétrage par défaut) la loi uniforme. Donc les résultats produits par des générateurs doivent "normalement" suivre les propriétés qui en découlent mathématiquement.

Je ne comprends pas pourquoi tu parles de loi multinomiale dans le cas que je présente : je précise qu'il y a 16 000 tirages pour chaque "face" (ce qui demande 16 * 16 000 tirages aléatoires, je n'ai pas fait dans l'économie )
Si j'avais fait 16 000 tirages (seulement) et que j'en avais déduit l'histogramme des effectifs sur les 16 "faces", oui, il y aurait loi multinomiale sur l'ensemble des 16 "faces", mais dans le programme, ce n'est pas le cas (chaque est vraiment indépendante des autres)

[gg0]

Léon,

je suis d'accord qu'un générateur pseudo-aléatoire devrait simuler la loi uniforme sur un certain intervalle d'entiers, mais cette simulation n'est pas aléatoire. Ce qui veut dire que selon la façon dont le générateur a été produit, la qualité est plus ou moins bonne. Et les défauts sont cumulatifs. Comme on ne sait rien du générateur utilisé, il est difficile de conclure. D'ailleurs, le fait que tu testes les "faces" dans l'ordre met peut-être plus en évidence les défauts de ce générateur.
Effectivement, le mot "multinomiale" n'est pas adapté; C'est la loi binomiale qui est simulée. Par contre, l'une des difficultés qu'on rencontre avec les simulateurs est que les tirages successifs ne sont pas réellement indépendants (connaissant le premier et le procédé, on trouve le deuxième). Ce qui fait que l'on ne teste pas vraiment le hasard parfait, mais l'équirépartition des valeurs. Ce que tu as fait avec les moyennes. On teste souvent aussi l'équirépartition de couples de valeurs.

Sans aucun renseignement sur le générateur utilisé, je ne me prononcerai pas, mais il y a suspicion.

Cordialement.

[leon1789]

Merci pour tes réponses.

je suis d'accord qu'un générateur pseudo-aléatoire devrait simuler la loi uniforme sur un certain intervalle d'entiers, mais cette simulation n'est pas aléatoire. Ce qui veut dire que selon la façon dont le générateur a été produit, la qualité est plus ou moins bonne.

oui, et c'est justement cette qualité qu'on essaie d'évaluer à travers ce simple test.

Et les défauts sont cumulatifs. Comme on ne sait rien du générateur utilisé, il est difficile de conclure.

On sait que le générateur est sensé simuler la loi uniforme, ça suffit pour faire des maths (obtenir les constantes théoriques) et tester le générateur par rapport à ces constantes.
Tu n'es pas d'accord ?

Par contre, l'une des difficultés qu'on rencontre avec les simulateurs est que les tirages successifs ne sont pas réellement indépendants (connaissant le premier et le procédé, on trouve le deuxième).

certes, il n'y a pas d'aléatoire dans la machine, mais un bon générateur doit simuler correctement cette indépendance théorique entre deux tirages. L'algo à l'intérieur du générateur est développé pour en donner l'illusion, sinon je change de crèmerie ! (l'indépendance des tirages est très importante, et mon programme ne le teste pas...)

Ce qui fait que l'on ne teste pas vraiment le hasard parfait, mais l'équirépartition des valeurs. Ce que tu as fait avec les moyennes. On teste souvent aussi l'équirépartition de couples de valeurs.

Dans le programme présenté, on teste à la fois la moyenne (1000 théoriquement), mais aussi l'écart-type (30.6 théoriquement) et finalement la dispersion de l'écart-type : c'est justement ce dernier point qui me pousse à douter du générateur dans sa fonction de simuler "correctement" la loi uniforme.

Sans aucun renseignement sur le générateur utilisé

Que voudrais-tu comme renseignement(s) sur le générateur, mis à part qu'il est vendu pour simuler la loi uniforme ? Je te pose cette question, non pas pour te renseigner (je n'ai aucun autre renseignement sur le générateur) mais pour comprendre quel autre renseignement te semble utile.

[gg0]

Léon,

On sait que le générateur est sensé simuler la loi uniforme, ça suffit pour faire des maths

On sait aussi que certains générateurs qui ont été utilisés sont inefficaces. Ou peu utilisables, comme le rand de Scilab. Ou celui de l'Amstrad CPC64.

mais un bon générateur doit simuler correctement cette indépendance théorique entre deux tirages

Dans tout ce que j'ai lu à ce propos, cette question n'est pas évoquée. Je ne sais d'ailleurs pas vraiment comment on pourrait le tester.
Justement, si le rand du programme est celui de Scilab, il est connu (et dit dans l'aide de scilab que je lisais il y a 10 ans) qu'il vaut mieux utiliser grand.

Un moyen de voir si un générateur est nettement imparfait est de générer une marche au hasard pixel par pixel en partant du centre de l'écran. Avec un générateur peu efficace, on sort rapidement de l'écran même en changeant le germe.

Cordialement.

[leon1789]

On sait aussi que certains générateurs qui ont été utilisés sont inefficaces. Ou peu utilisables, comme le rand de Scilab. Ou celui de l'Amstrad CPC64.

c'est justement détecter cela (ou ceux-là ) de manière simple qui m'intéresse.

Dans tout ce que j'ai lu à ce propos, cette question n'est pas évoquée. Je ne sais d'ailleurs pas vraiment comment on pourrait le tester.

un premier truc, pas suffisant certes, ce serait calculer la co-variance.

[gg0]

ce serait calculer la co-variance.

Tu veux sans doute parler d'une forme d'auto-corrélation ? C'est une piste, effectivement.

Cordialement.