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

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é et cela fonctionne parfaitement.

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

Premièrement, comme plusieurs réponses l'ont souligné, si les bits de poids faible de rand()ne sont pas convenablement uniformes, le résultat de l'opérateur de reste n'est pas non plus uniforme.

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

À titre d'exemple extrême, supposez que cela rand()produit des valeurs uniformément distribuées dans la plage [0..6]. Si vous regardez les restes de ces valeurs, lorsque rand()renvoie une valeur de la plage [0..5], le reste produit des résultats uniformément répartis dans la plage [0..5]. Lorsque rand()renvoie 6, rand() % 6renvoie 0, comme si rand()avait renvoyé 0. Vous obtenez donc une distribution avec deux fois plus de 0 que toute autre valeur.

Le second est le vrai problème avec rand() % 6.

Le moyen d'éviter ce problème consiste à ignorer les valeurs 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 vous rand()renvoie une valeur supérieure ou égale à ce multiple, vous la rejetez et appelez à nouveau `rand (), autant de fois que nécessaire.

Alors:

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

C'est une implémentation différente du code en question, destinée à montrer plus clairement ce qui se passe.

Il y a des profondeurs cachées ici:

1. L'utilisation du petit uin RAND_MAX + 1u. RAND_MAXest défini comme unint type et est souvent le plus grand possible int. Le comportement de RAND_MAX + 1serait indéfini dans les cas où vous déborderiez d'un signedtype. L'écriture 1uforce la conversion de type de RAND_MAXto unsigned, évitant ainsi le débordement.

2. L'utilisation de % 6 can (mais sur chaque implémentation de ce std::randque j'ai vu ne le fait pas ) introduit un biais statistique supplémentaire au-delà de l'alternative présentée. De tels cas où il % 6est dangereux sont les cas où le générateur de nombres a des plaines de corrélation dans les bits de poids faible, comme une implémentation IBM plutôt célèbre (en C) des randannées 1970, je pense, qui a inversé les bits haut et bas comme "un final fleurir". Une autre considération est que 6 est très petit cf. RAND_MAX, il y aura donc un effet minimal si ce RAND_MAXn'est pas un multiple de 6, ce qui n'est probablement pas le cas.

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

Cet exemple de code illustre qu'il std::rands'agit d'un cas de balderdash culte du fret hérité qui devrait faire lever vos sourcils à chaque fois que vous le voyez.

Il y a plusieurs problèmes ici:

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

Le premier problème est que le contrat supposé, des échantillons aléatoires uniformes indépendants dans chaque appel, n'est pas réellement ce que dit la documentation - et dans la pratique, les implémentations ont historiquement échoué à fournir même le plus simple simulacre d'indépendance. Par exemple, C99 §7.20.2.1 'La randfonction' dit, sans élaboration:

La randfonction calcule une séquence d'entiers pseudo-aléatoires compris entre 0 et RAND_MAX.

C'est une phrase dénuée de sens, car la pseudo-aléatoire est une propriété d'une fonction (ou d'une famille de fonctions ), pas d'un entier, mais cela n'empêche pas même les bureaucrates de l'ISO d'abuser du langage. Après tout, les seuls lecteurs qui en seraient contrariés savent mieux que de lire la documentation randpar crainte de voir leurs cellules cérébrales se décomposer.

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

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;
}

Cela 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 entiers pairs et impairs dans les appels consécutifs - après

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 les échantillons aléatoires indépendants sur toute distribution prise en charge sur des entiers pairs et impairs. Ainsi, un culte de la cargaison a émergé selon lequel il fallait se débarrasser des bits de poids faible pour chasser la bête insaisissable du «meilleur hasard». (Alerte spoiler: ce n'est pas un terme technique. Ceci est un signe que la prose que vous lisez ne sait pas de quoi elle parle, ou pense que vous n'avez aucune idée et doit être condescendante.)

Le deuxième problème est que même si chaque appel échantillonnait indépendamment d'une distribution aléatoire uniforme sur 0, 1, 2,… RAND_MAX, le résultat de rand() % 6ne serait pas uniformément distribué en 0, 1, 2, 3, 4, 5 comme un dé rouler, sauf si elle RAND_MAXest congruente à -1 modulo 6. Contre-exemple simple: SiRAND_MAX = 6, alors à partir de rand(), tous les résultats ont une probabilité égale de 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 consiste à utiliser un échantillonnage de rejet: tirez à plusieurs reprises un échantillon aléatoire uniforme indépendant sde 0, 1, 2,… RAND_MAX, et rejetez (par exemple) les résultats 0, 1, 2,…, ((RAND_MAX + 1) % 6) - 1- si vous obtenez l'un des ceux-là, recommencer; sinon, cédez 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 rand()que nous acceptons est divisible par 6, et chaque résultat possible s % 6est obtenu par le même nombre de résultats acceptésrand() , donc si rand()est uniformément distribué, il en est de même s. Il n'y a pas de limite sur le nombre d'essais, mais le nombre attendu est inférieur à 2 et la probabilité de succès augmente de façon exponentielle avec le nombre d'essais.

Le choix dont les résultats de rand()vous rejetez est sans importance, à condition que vous associez un nombre égal d'entre eux à chaque entier inférieur à 6. Le code à cppreference.com fait un autre choix, en raison du premier problème ci-dessus que rien est garanti sur le la distribution ou l'indépendance des sorties de rand(), et en pratique, les bits de poids faible présentaient des modèles qui ne «semblaient pas assez aléatoires» (sans oublier que la sortie suivante est une fonction déterministe de la précédente).

Exercice pour le lecteur: Démontrer que le code à cppreference.com produit une distribution uniforme sur les rouleaux de matrice se rand()produit une distribution uniforme sur 0, 1, 2, ..., RAND_MAX.

Exercice pour le lecteur: Pourquoi préféreriez-vous que l'un ou l'autre sous-ensemble soit rejeté? Quel calcul est nécessaire pour chaque essai dans les deux cas?

Un troisième problème est que l'espace de départ est si petit que même si la graine est uniformément distribuée, un adversaire armé de la connaissance de votre programme et d'un résultat, mais pas de la graine, peut facilement prédire la graine et les résultats ultérieurs, ce qui les fait paraître non. aléatoire après tout. Alors ne pensez même pas à l'utiliser pour la cryptographie.

Vous pouvez emprunter la voie sophistiquée et la std::uniform_int_distributionclasse C ++ 11 avec un appareil aléatoire approprié et votre moteur aléatoire préféré comme le toujours populaire Mersenne Twister std::mt19937pour jouer aux dés avec votre cousin de quatre ans, mais même cela ne va pas être apte à générer du matériel de clé cryptographique - et le twister de Mersenne est également un espace terrible avec un état de plusieurs kilo-octets qui ravage le cache de votre processeur avec un temps de configuration obscène, il est donc mauvais, même pour, par exemple , des simulations de Monte Carlo parallèles avec arbres reproductibles de sous-calculs; sa popularité découle probablement principalement de son nom accrocheur. Mais vous pouvez l'utiliser pour lancer des dés jouets comme cet exemple!

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

Je ne suis en aucun cas un utilisateur expérimenté de C ++, mais j'étais intéressé de voir si les autres réponses concernant le fait d' std::rand()/((RAND_MAX + 1u)/6)être moins biaisé que ce sont 1+std::rand()%6réellement vraies. J'ai donc écrit un programme de test pour tabuler les résultats pour les deux méthodes (je n'ai pas écrit C ++ depuis longtemps, veuillez le vérifier). Un lien pour exécuter le code se trouve ici . 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 la sortie de ceci et utilisé la chisq.testfonction dans R pour exécuter un test du chi carré pour voir si les résultats sont significativement différents de ceux attendus. Cette question de stackexchange va plus en détail sur l'utilisation du 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 typiques utilisées pour tester la signification (0,05). Cela signifie que nous ne considérons ni l'un ni l'autre comme étant biaisé. Il est intéressant de noter que la méthode supposée non biaisée a systématiquement des valeurs p plus faibles, ce qui indique qu'elle pourrait en fait être plus biaisée. La mise en garde étant que je n'ai fait que 3 courses.

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, alors je vais en parler ici. Tout d'abord, l'essentiel 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 changez le nombre d'essais, vous obtenez également une variété de résultats. Essayez d'augmenter ou de diminuer d'un ordre de grandeur pour voir ce que je veux dire. Troisièmement, il y a une troncature ou un arrondi d'entiers lorsque les valeurs attendues ne sont pas tout à fait exactes. Ce n'est probablement pas suffisant pour faire une différence, mais c'est là.

Fondamentalement, en résumé, il a juste eu la bonne graine et le bon nombre d'essais qu'il pourrait obtenir un faux résultat.

On peut penser à un générateur de nombres aléatoires comme travaillant sur un flux de chiffres binaires. Le générateur transforme le flux en nombres en le découpant en morceaux. Si la std:randfonction fonctionne avec un RAND_MAXde 32767, alors elle utilise 15 bits dans chaque tranche.

Quand on prend les modules d'un nombre compris entre 0 et 32767 inclus, on trouve que 5462 «0» et «1» mais seulement 5461 «2», «3», «4» et «5». Par conséquent, le résultat est biaisé. Plus la valeur RAND_MAX est élevée, moins il y aura de biais, mais c'est inéluctable.

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

On espère que chaque bit du train de bits a une chance égale d'être un «0» ou un «1» indépendamment de l'endroit où il se trouve dans le flux ou des valeurs des autres bits. Ceci est exceptionnellement difficile en pratique. Les nombreuses implémentations différentes des PRNG logiciels offrent différents compromis entre vitesse et qualité. Un générateur congruentiel linéaire tel que std::randoffre la vitesse la plus rapide pour une qualité la plus basse. Un générateur cryptographique offre la plus haute qualité pour la vitesse la plus basse.