Pourquoi rand()%6 est-il biaisé ?

Question

En lisant comment utiliser std::rand, j'ai trouvé ce code sur cppreference.com

int x = 7;
while(x > 6) 
    x = 1 + std::rand()/((RAND_MAX + 1u)/6);  // Note: 1+rand()%6 is biased

Qu'est-ce qui ne va pas avec l'expression de droite ? Je l'ai essayée et elle fonctionne parfaitement.

Answer 1

Il y a deux problèmes avec rand() % 6 (le 1+ n'affecte aucun des deux problèmes).

Tout d'abord, comme plusieurs réponses l'ont souligné, si les bits de poids faibles de rand() ne sont pas uniformes de manière appropriée, le résultat de l'opérateur de reste n'est pas non plus uniforme.

Deuxièmement, si le nombre de valeurs distinctes produit par rand() n'est pas un multiple de 6, alors le reste produira plus de valeurs basses que de valeurs hautes. C'est vrai même si rand() renvoie des valeurs parfaitement distribuées.

A titre d'exemple extrême, supposons que rand() produit des valeurs uniformément distribuées dans l'intervalle [0..6] . Si vous regardez les restes pour ces valeurs, quand rand() renvoie une valeur dans l'intervalle [0..5] le reste produit des résultats uniformément répartis dans l'intervalle [0..5] . Quand rand() rendements 6, rand() % 6 renvoie 0, tout comme si rand() avait donné 0. On obtient donc une distribution avec deux fois plus de 0 que toute autre valeur.

Le second est le réel problème avec rand() % 6 .

Pour éviter ce problème, il faut jeter qui produiraient des doublons non uniformes. Vous calculez le plus grand multiple de 6 qui est inférieur ou égal à RAND_MAX et chaque fois que rand() retourne une valeur supérieure ou égale à ce multiple, vous le rejetez et appelez `rand() à nouveau, autant de fois que nécessaire.

Donc :

int max = 6 * ((RAND_MAX + 1u) / 6)
int value = rand();
while (value >= max)
    value = rand();

Il s'agit d'une implémentation différente du code en question, destinée à montrer plus clairement ce qui se passe.

Answer 2

Il y a des profondeurs cachées ici :

1. L'utilisation de la petite u en RAND_MAX + 1u . RAND_MAX est défini comme étant un int et est souvent le plus grand type de int . Le comportement de RAND_MAX + 1 serait indéfini dans de telles circonstances que vous débordiez d'un signed type. Écriture 1u force la conversion de type de RAND_MAX a unsigned ce qui évite le débordement.

2. L'utilisation de % 6 peut (mais sur chaque implémentation de std::rand J'ai vu n'a pas ) introduisent un biais statistique supplémentaire par rapport à l'alternative présentée. Les cas où % 6 est dangereux sont les cas où le générateur de nombres a des plaines de corrélation dans les bits de bas ordre, comme une implémentation IBM assez célèbre (en C) de rand dans les années 1970, je crois, qui inversait les bits de poids fort et de poids faible pour donner une touche finale. Une autre considération est que 6 est très petit cf. RAND_MAX Il y aura donc un effet minimal si RAND_MAX n'est pas un multiple de 6, ce qui n'est probablement pas le cas.

En conclusion, ces jours-ci, en raison de sa tractabilité, j'utiliserais % 6 . Il est peu probable qu'il introduise des anomalies statistiques autres que celles introduites par le générateur lui-même. Si vous avez encore des doutes, test votre générateur pour voir s'il possède les propriétés statistiques appropriées pour votre cas d'utilisation.

Answer 3

Cet exemple de code illustre que std::rand est un cas de balivernes de culte du cargo hérité qui devrait vous faire lever les sourcils à chaque fois que vous le verrez.

Il y a plusieurs problèmes ici :

Le contrat que les gens supposent généralement - même les pauvres âmes infortunées qui ne savent pas mieux et ne pensent pas précisément en ces termes - est que rand des échantillons de la distribution uniforme sur les entiers en 0, 1, 2, , RAND_MAX et chaque appel donne lieu à un indépendant échantillon.

Le premier problème est que le contrat supposé, à savoir des échantillons aléatoires uniformes indépendants à chaque appel, n'est pas réellement ce que dit la documentation - et dans la pratique, les implémentations n'ont historiquement pas réussi à fournir même le plus petit simulacre d'indépendance. Par exemple, dans le document C99 §7.20.2.1, "L'autorité de surveillance de l'environnement". rand fonction" dit, sans élaborer :

El rand calcule une séquence d'entiers pseudo-aléatoires dans l'intervalle 0 à RAND_MAX .

Cette phrase n'a aucun sens, car le caractère pseudo-aléatoire est une propriété d'une fonction (ou famille de fonctions ), et non d'un entier, mais cela n'empêche pas les bureaucrates de l'ISO d'abuser de ce langage. Après tout, les seuls lecteurs que cela pourrait déranger savent qu'il vaut mieux ne pas lire la documentation de la norme rand de peur que les cellules de leur cerveau ne se décomposent.

Une implémentation historique typique en C fonctionne comme suit :

static unsigned int seed = 1;

static void
srand(unsigned int s)
{
    seed = s;
}

static unsigned int
rand(void)
{
    seed = (seed*1103515245 + 12345) % ((unsigned long)RAND_MAX + 1);
    return (int)seed;
}

Ceci a la propriété malheureuse que même si un seul échantillon peut être uniformément distribué sous une graine aléatoire uniforme (qui dépend de la valeur spécifique de RAND_MAX ), il alterne entre les nombres entiers pairs et impairs dans des appels consécutifs.

int a = rand();
int b = rand();

l'expression (a & 1) ^ (b & 1) donne 1 avec une probabilité de 100%, ce qui n'est pas le cas pour indépendant échantillons aléatoires sur n'importe quelle distribution supportée sur les entiers pairs et impairs. C'est ainsi qu'est apparu un culte du fret selon lequel il fallait se débarrasser des bits d'ordre inférieur pour chasser la bête insaisissable du "meilleur hasard". (Alerte spoiler : ce n'est pas un terme technique. C'est un signe que la personne dont vous lisez la prose ne sait pas de quoi elle parle, ou pense que le terme "meilleur hasard" n'a pas de sens. vous sont ignorants et doivent être traités avec condescendance).

Le deuxième problème est que même si chaque appel a été échantillonné indépendamment à partir d'une distribution aléatoire uniforme. sur 0, 1, 2, , RAND_MAX le résultat de rand() % 6 ne serait pas uniformément distribué en 0, 1, 2, 3, 4, 5 comme un jet de dé, sauf si RAND_MAX est congruente à -1 modulo 6. Contre-exemple simple : Si RAND_MAX = 6, alors de rand() toutes les issues ont la même probabilité 1/7, mais à partir de rand() % 6 le résultat 0 a une probabilité de 2/7 tandis que tous les autres résultats ont une probabilité de 1/7.

La bonne façon de procéder est l'échantillonnage par rejet : à plusieurs reprises tirer un échantillon aléatoire uniforme indépendant s de 0, 1, 2, , RAND_MAX et rejeter (par exemple) les résultats 0, 1, 2, , ((RAND_MAX + 1) % 6) - 1 Si vous obtenez l'un de ces résultats, recommencez ; sinon, rendez-le. s % 6 .

unsigned int s;
while ((s = rand()) < ((unsigned long)RAND_MAX + 1) % 6)
    continue;
return s % 6;

De cette façon, l'ensemble des résultats de rand() que nous acceptons est divisible de manière égale par 6, et chaque résultat possible de s % 6 est obtenu par le même nombre de accepté les résultats de rand() donc si rand() est uniformément distribué, alors il en est de même pour s . Il n'y a pas de lié sur le nombre d'essais, mais le nombre attendu est inférieure à 2, et la probabilité de réussite croît de façon exponentielle avec le nombre d'essais.

Le choix de dont les résultats de rand() que vous rejetez n'a pas d'importance, à condition que vous en fassiez correspondre un nombre égal à chaque nombre entier inférieur à 6. Le code sur cppreference.com fait un différents en raison du premier problème ci-dessus, à savoir que rien n'est garanti quant à la distribution ou à l'indépendance des résultats de l'opération. rand() Dans la pratique, les bits de poids faible présentent des motifs qui ne semblent pas suffisamment aléatoires (sans compter que la sortie suivante est une fonction déterministe de la précédente).

Exercice pour le lecteur : Prouvez que le code de cppreference.com produit une distribution uniforme sur les jets de dé si rand() donne une distribution uniforme sur 0, 1, 2, , RAND_MAX .

Exercice pour le lecteur : Pourquoi préférer l'un ou l'autre des sous-ensembles à rejeter ? Quel calcul est nécessaire pour chaque essai dans les deux cas ?

Un troisième problème est que l'espace des graines est si petit que même si la graine est distribuée uniformément, un adversaire connaissant votre programme et un résultat mais pas la graine peut facilement prédire la graine et les résultats suivants, ce qui les fait paraître moins aléatoires après tout. Donc ne pensez même pas à l'utiliser pour la cryptographie.

Vous pouvez opter pour la voie de la sur-ingénierie fantaisiste et l'approche C++11. std::uniform_int_distribution avec un dispositif aléatoire approprié et votre moteur aléatoire préféré, comme le toujours populaire twister de Mersenne. std::mt19937 pour jouer aux dés avec votre cousin de quatre ans, mais même cela n'est pas adapté à la génération de clés cryptographiques - et le tordeur de Mersenne est également très gourmand en espace, avec un état de plusieurs kilo-octets qui fait des ravages dans le cache de votre processeur et un temps de configuration obscène, donc il n'est pas bon même pour.., par exemple Ce logiciel permet de réaliser des simulations de Monte Carlo parallèles avec des arbres de sous-calculs reproductibles ; sa popularité est probablement due à son nom accrocheur. Mais vous pouvez l'utiliser pour lancer des dés comme cet exemple !

Une autre approche consiste à utiliser un générateur de nombres pseudo-aléatoires cryptographique simple avec un petit état, tel qu'un simple effacement rapide des clés PRNG ou simplement un chiffrement par flux tel que AES-CTR ou ChaCha20 si vous êtes sûr de vous ( par exemple dans une simulation de Monte Carlo pour la recherche en sciences naturelles) qu'il n'y a pas de conséquences négatives à prévoir les résultats passés si l'état est jamais compromis.

Answer 4

Je ne suis en aucun cas un utilisateur expérimenté de C++, mais j'étais intéressé de voir si les autres réponses relatives à std::rand()/((RAND_MAX + 1u)/6) étant moins biaisé que 1+std::rand()%6 est en fait vrai. J'ai donc écrit un programme de test pour tabuler les résultats des deux méthodes (je n'ai pas écrit de C++ depuis des lustres, merci de le vérifier). Un lien pour exécuter le code se trouve aquí . Il est également reproduit comme suit :

// Example program
#include <cstdlib>
#include <iostream>
#include <ctime>
#include <string>

int main()
{
    std::srand(std::time(nullptr)); // use current time as seed for random generator

    // Roll the die 6000000 times using the supposedly unbiased method and keep track of the results

    int results[6] = {0,0,0,0,0,0};

    // roll a 6-sided die 20 times
    for (int n=0; n != 6000000; ++n) {
        int x = 7;
        while(x > 6) 
            x = 1 + std::rand()/((RAND_MAX + 1u)/6);  // Note: 1+rand()%6 is biased

        results[x-1]++;
    }

    for (int n=0; n !=6; n++) {
        std::cout << results[n] << ' ';
    }

    std::cout << "\n";

    // Roll the die 6000000 times using the supposedly biased method and keep track of the results

    int results_bias[6] = {0,0,0,0,0,0};

    // roll a 6-sided die 20 times
    for (int n=0; n != 6000000; ++n) {
        int x = 7;
        while(x > 6) 
            x = 1 + std::rand()%6;

        results_bias[x-1]++;
    }

    for (int n=0; n !=6; n++) {
        std::cout << results_bias[n] << ' ';
    }
}

J'ai ensuite pris le résultat de cette opération et utilisé la fonction chisq.test dans R pour exécuter un test du chi-deux afin de voir si les résultats sont significativement différents de ceux attendus. Cette question de stackexchange explique plus en détail comment utiliser le test du chi carré pour tester l'équité du dé : Comment puis-je tester si un dé est juste ? . Voici les résultats de quelques essais :

> ?chisq.test
> unbias <- c(100150, 99658, 100319, 99342, 100418, 100113)
> bias <- c(100049, 100040, 100091, 99966, 100188, 99666 )

> chisq.test(unbias)

Chi-squared test for given probabilities

data:  unbias
X-squared = 8.6168, df = 5, p-value = 0.1254

> chisq.test(bias)

Chi-squared test for given probabilities

data:  bias
X-squared = 1.6034, df = 5, p-value = 0.9008

> unbias <- c(998630, 1001188, 998932, 1001048, 1000968, 999234 )
> bias <- c(1000071, 1000910, 999078, 1000080, 998786, 1001075   )
> chisq.test(unbias)

Chi-squared test for given probabilities

data:  unbias
X-squared = 7.051, df = 5, p-value = 0.2169

> chisq.test(bias)

Chi-squared test for given probabilities

data:  bias
X-squared = 4.319, df = 5, p-value = 0.5045

> unbias <- c(998630, 999010, 1000736, 999142, 1000631, 1001851)
> bias <- c(999803, 998651, 1000639, 1000735, 1000064,1000108)
> chisq.test(unbias)

Chi-squared test for given probabilities

data:  unbias
X-squared = 7.9592, df = 5, p-value = 0.1585

> chisq.test(bias)

Chi-squared test for given probabilities

data:  bias
X-squared = 2.8229, df = 5, p-value = 0.7273

Dans les trois essais que j'ai effectués, la valeur p des deux méthodes était toujours supérieure aux valeurs alpha généralement utilisées pour tester la signification (0,05). Cela signifie que nous ne pourrions pas considérer qu'aucune des deux méthodes n'est biaisée. Il est intéressant de noter que la méthode supposée non biaisée présente des valeurs p systématiquement inférieures, ce qui indique qu'elle pourrait en fait être plus biaisée. La mise en garde étant que je n'ai fait que 3 essais.

MISE À JOUR : Pendant que j'écrivais ma réponse, Konrad Rudolph a publié une réponse qui adopte la même approche, mais obtient un résultat très différent. Je n'ai pas la réputation de commenter sa réponse, je vais donc l'aborder ici. Tout d'abord, l'élément principal est que le code qu'il utilise utilise la même graine pour le générateur de nombres aléatoires à chaque fois qu'il est exécuté. Si vous changez la graine, vous obtenez en fait une variété de résultats. Deuxièmement, si vous ne changez pas la graine, mais que vous modifiez le nombre d'essais, vous obtenez également des résultats variés. Essayez d'augmenter ou de diminuer d'un ordre de grandeur pour voir ce que je veux dire. Troisièmement, il y a une certaine troncature ou un arrondi des nombres entiers qui fait que les valeurs attendues ne sont pas tout à fait exactes. Ce n'est probablement pas suffisant pour faire une différence, mais c'est là.

En gros, pour résumer, il s'est avéré que la graine et le nombre d'essais étaient les bons pour obtenir un faux résultat.

Answer 5

On peut considérer qu'un générateur de nombres aléatoires travaille sur un flux de chiffres binaires. Le générateur transforme le flux en chiffres en le découpant en morceaux. Si le std:rand fonctionne avec une RAND_MAX de 32767, alors il utilise 15 bits dans chaque tranche.

Lorsque l'on prend les modules d'un nombre compris entre 0 et 32767 inclus, on trouve 5462 '0' et '1' mais seulement 5461 '2', '3', '4' et '5'. Le résultat est donc faussé. Plus la valeur de RAND_MAX est grande, moins il y aura de biais, mais il est inévitable.

Ce qui n'est pas biaisé est un nombre dans la plage [0..(2^n)-1]. Vous pouvez générer un nombre (théoriquement) meilleur dans la plage 0..5 en extrayant 3 bits, en les convertissant en un nombre entier dans la plage 0..7 et en rejetant 6 et 7.

On espère que chaque bit du flux binaire a une chance égale d'être un "0" ou un "1", indépendamment de sa position dans le flux ou des valeurs des autres bits. Cela est exceptionnellement difficile dans la pratique. Les nombreuses implémentations des PRNG logiciels offrent différents compromis entre vitesse et qualité. Un générateur congruentiel linéaire tel que std::rand offre la vitesse la plus rapide pour la qualité la plus basse. Un générateur cryptographique offre la meilleure qualité pour la vitesse la plus faible.