Por que rand ()% 6 é tendencioso?

Ao ler como usar std :: rand, encontrei este código em cppreference.com

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

O que há de errado com a expressão à direita? Tentei e funciona perfeitamente.

Existem dois problemas com rand() % 6(o 1+não afeta nenhum dos problemas).

Em primeiro lugar, como várias respostas apontaram, se os bits baixos de rand()não forem apropriadamente uniformes, o resultado do operador restante também não será uniforme.

Em segundo lugar, se o número de valores distintos produzidos por rand()não for um múltiplo de 6, o restante produzirá mais valores baixos do que valores altos. Isso é verdade mesmo se rand()retornar valores perfeitamente distribuídos.

Como um exemplo extremo, imagine que rand()produz valores uniformemente distribuídos no intervalo [0..6]. Se você olhar para os restantes para esses valores, quando rand()retorna um valor no intervalo [0..5], o restante produz resultados uniformemente distribuídos no intervalo [0..5]. Quando rand()retorna 6, rand() % 6retorna 0, como se rand()tivesse retornado 0. Assim, você obtém uma distribuição com o dobro de 0's que qualquer outro valor.

O segundo é o verdadeiro problema com rand() % 6.

A maneira de evitar esse problema é descartar valores que produziriam duplicatas não uniformes. Você calcula o maior múltiplo de 6 que é menor ou igual a RAND_MAXe sempre que rand()retorna um valor maior ou igual a esse múltiplo, você o rejeita e chama `rand () novamente, quantas vezes forem necessárias.

Assim:

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

Essa é uma implementação diferente do código em questão, destinada a mostrar mais claramente o que está acontecendo.

Existem profundidades escondidas aqui:

1. O uso do pequeno uem RAND_MAX + 1u. RAND_MAXé definido como um inttipo e geralmente é o maior possível int. O comportamento de RAND_MAX + 1seria indefinido em tais casos, pois você estaria transbordando um signedtipo. A gravação 1uforça a conversão de tipo de RAND_MAXpara unsigned, evitando assim o estouro.

2. O uso de % 6 lata (mas em todas as implementações do std::randque eu vi não ) introduzir qualquer viés estatístico adicional acima e além da alternativa apresentada. Os casos em que % 6é perigoso são os casos em que o gerador de número tem planos de correlação nos bits de ordem inferior, como uma implementação bastante famosa da IBM (em C) rand, eu acho, da década de 1970, que inverteu os bits superior e inferior como "um final florescer". Uma consideração adicional é que 6 é muito pequeno cf. RAND_MAX, então haverá um efeito mínimo se RAND_MAXnão for um múltiplo de 6, o que provavelmente não é.

Em conclusão, hoje em dia, devido à sua tratabilidade, eu usaria % 6. Não é provável que introduza quaisquer anomalias estatísticas além daquelas introduzidas pelo próprio gerador. Se você ainda estiver em dúvida, teste seu gerador para ver se ele possui as propriedades estatísticas apropriadas para seu caso de uso.

Este código de exemplo ilustra que std::randé um caso de balderdash de culto de carga legado que deve fazer você levantar as sobrancelhas toda vez que você vê-lo.

Existem várias questões aqui:

O pessoal do contrato geralmente assume - mesmo as pobres almas infelizes que não conhecem nada melhor e não pensam nisso precisamente nestes termos - é que as randamostras da distribuição uniforme nos inteiros em 0, 1, 2, ... RAND_MAX,, e cada chamada produz uma amostra independente .

O primeiro problema é que o contrato assumido, amostras aleatórias uniformes independentes em cada chamada, não é realmente o que diz a documentação - e, na prática, as implementações historicamente falharam em fornecer nem mesmo o mais básico simulacro de independência. Por exemplo, C99 §7.20.2.1 'A randfunção' diz, sem elaboração:

A randfunção calcula uma sequência de inteiros pseudoaleatórios no intervalo de 0 a RAND_MAX.

Esta é uma frase sem sentido, porque a pseudo-aleatoriedade é uma propriedade de uma função (ou família de funções ), não de um número inteiro, mas isso não impede nem mesmo os burocratas da ISO de abusar da linguagem. Afinal, os únicos leitores que ficariam chateados sabem que não é melhor ler a documentação randpor medo de que suas células cerebrais se deteriorem.

Uma implementação histórica típica em C funciona assim:

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

Isso tem a propriedade infeliz de que , embora uma única amostra possa ser uniformemente distribuída sob uma semente aleatória uniforme (que depende do valor específico de RAND_MAX), ela alterna entre inteiros pares e ímpares em chamadas consecutivas - após

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

a expressão (a & 1) ^ (b & 1)produz 1 com 100% de probabilidade, o que não é o caso para amostras aleatórias independentes em qualquer distribuição suportada em inteiros pares e ímpares. Assim, surgiu um culto à carga em que se deve descartar os bits de ordem inferior para perseguir a besta indescritível de 'melhor aleatoriedade'. (Alerta de spoiler: este não é um termo técnico. Este é um sinal de que qualquer prosa que você está lendo não sabe do que está falando ou pensa que você é um ignorante e deve ser condescendente.)

O segundo problema é que mesmo se cada chamada fosse amostrada independentemente de uma distribuição aleatória uniforme em 0, 1, 2, ... RAND_MAX, o resultado de rand() % 6não seria uniformemente distribuído em 0, 1, 2, 3, 4, 5 como um dado role, a menos que RAND_MAXseja congruente com -1 módulo 6. Contra-exemplo simples: Se RAND_MAX= 6, então de rand(), todos os resultados têm probabilidade igual 1/7, mas de rand() % 6, o resultado 0 tem probabilidade 2/7, enquanto todos os outros resultados têm probabilidade 1/7 .

A maneira certa de fazer isso é com a amostragem de rejeição: retire repetidamente uma amostra aleatória uniforme independente sde 0, 1, 2, ... RAND_MAX, e rejeite (por exemplo) os resultados 0, 1, 2, ..., ((RAND_MAX + 1) % 6) - 1- se você obtiver um dos aqueles, recomece; caso contrário, rendimento s % 6.

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

Dessa forma, o conjunto de resultados rand()que aceitamos é igualmente divisível por 6, e cada resultado possível de s % 6é obtido pelo mesmo número de resultados aceitos de rand(), portanto, se rand()for uniformemente distribuído, então o é s. Não há limite para o número de tentativas, mas o número esperado é menor que 2, e a probabilidade de sucesso aumenta exponencialmente com o número de tentativas.

A escolha de quais resultados rand()você rejeita é imaterial, desde que você mapeie um número igual deles para cada número inteiro abaixo de 6. O código em cppreference.com faz uma escolha diferente , por causa do primeiro problema acima - que nada é garantido sobre o distribuição ou independência de saídas de rand(), e na prática, os bits de ordem inferior exibiram padrões que não 'parecem aleatórios o suficiente' (não importa que a próxima saída seja uma função determinística da anterior).

Exercício para o leitor: Prove que o código em cppreference.com produz uma distribuição uniforme nas jogadas de dados se rand()produz uma distribuição uniforme em 0, 1, 2,… RAND_MAX,.

Exercício para o leitor: Por que você prefere rejeitar um ou outro subconjunto? Qual cálculo é necessário para cada tentativa nos dois casos?

Um terceiro problema é que o espaço da semente é tão pequeno que, mesmo que a semente seja distribuída uniformemente, um adversário armado com conhecimento do seu programa e um resultado, mas não a semente, pode prever prontamente a semente e os resultados subsequentes, o que os faz parecer que não afinal aleatório. Portanto, nem pense em usar isso para criptografia.

Você pode seguir o caminho sofisticado da superengenharia e as std::uniform_int_distributionaulas de C ++ 11 com um dispositivo aleatório apropriado e seu mecanismo aleatório favorito como o sempre popular tornado de Mersenne std::mt19937para jogar dados com seu primo de quatro anos, mas mesmo isso não vai estar apto para gerar material de chave criptográfica - e o Mersenne twister é um terrível devorador de espaço também com um estado de vários kilobytes causando estragos no cache de sua CPU com um tempo de configuração obsceno, por isso é ruim mesmo para, por exemplo , simulações de Monte Carlo paralelas com árvores reproduzíveis de subcomputações; sua popularidade provavelmente decorre principalmente de seu nome atraente. Mas você pode usá-lo para rolar dados de brinquedo como este exemplo!

Outra abordagem é usar um gerador de números pseudo-aleatórios criptográficos simples com um estado pequeno, como um PRNG de apagamento rápido de chave simples ou apenas uma cifra de fluxo, como AES-CTR ou ChaCha20, se você estiver confiante ( por exemplo , em uma simulação de Monte Carlo para pesquisa em ciências naturais) de que não há consequências adversas em prever resultados passados ​​se o estado for comprometido.

Não sou um usuário experiente de C ++ de forma alguma, mas estava interessado em ver se as outras respostas sobre std::rand()/((RAND_MAX + 1u)/6)ser menos tendencioso do que 1+std::rand()%6realmente são verdadeiras. Então, escrevi um programa de teste para tabular os resultados de ambos os métodos (não escrevo C ++ há anos, verifique). Um link para executar o código pode ser encontrado aqui . Também é reproduzido da seguinte forma:

// 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] << ' ';
    }
}

Em seguida, peguei o resultado disso e usei a chisq.testfunção em R para executar um teste de qui-quadrado para ver se os resultados são significativamente diferentes do esperado. Esta questão de troca de pilha apresenta mais detalhes sobre o uso do teste qui-quadrado para testar a justiça do dado: Como posso testar se um dado é justo? . Aqui estão os resultados de algumas execuções:

> ?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

Nas três execuções que fiz, o valor-p para ambos os métodos foi sempre maior do que os valores alfa típicos usados ​​para testar a significância (0,05). Isso significa que não consideraríamos nenhum deles tendencioso. Curiosamente, o método supostamente não tendencioso tem valores de p consistentemente mais baixos, o que indica que ele pode realmente ser mais tendencioso. A ressalva é que eu fiz apenas 3 execuções.

ATUALIZAÇÃO: Enquanto eu escrevia minha resposta, Konrad Rudolph postou uma resposta que segue a mesma abordagem, mas obtém um resultado muito diferente. Não tenho reputação de comentar sua resposta, então vou abordá-la aqui. Primeiro, o principal é que o código que ele usa usa a mesma semente para o gerador de números aleatórios toda vez que é executado. Se você mudar a semente, você realmente obterá uma variedade de resultados. Em segundo lugar, se você não mudar a semente, mas mudar o número de tentativas, também obterá uma variedade de resultados. Tente aumentar ou diminuir em uma ordem de magnitude para ver o que quero dizer. Terceiro, há algum truncamento ou arredondamento de número inteiro em que os valores esperados não são muito precisos. Provavelmente não é o suficiente para fazer a diferença, mas está lá.

Basicamente, em resumo, ele apenas obteve a semente certa e o número de tentativas que pode estar obtendo um resultado falso.

Pode-se pensar em um gerador de números aleatórios trabalhando em um fluxo de dígitos binários. O gerador transforma o fluxo em números, dividindo-o em pedaços. Se a std:randfunção estiver funcionando com um RAND_MAXde 32767, ela estará usando 15 bits em cada fatia.

Quando se pega os módulos de um número entre 0 e 32.767, inclusive, descobre-se 5462 '0's e' 1's, mas apenas 5461 '2's,' 3's, '4's e' 5's. Portanto, o resultado é tendencioso. Quanto maior for o valor de RAND_MAX, menor será o viés, mas é inevitável.

O que não é tendencioso é um número no intervalo [0 .. (2 ^ n) -1]. Você pode gerar um número (teoricamente) melhor no intervalo 0..5 extraindo 3 bits, convertendo-os em um número inteiro no intervalo 0..7 e rejeitando 6 e 7.

Espera-se que cada bit no fluxo de bits tenha uma chance igual de ser '0' ou '1', independentemente de onde ele esteja no fluxo ou dos valores de outros bits. Isso é excepcionalmente difícil na prática. As muitas implementações diferentes de PRNGs de software oferecem diferentes compromissos entre velocidade e qualidade. Um gerador congruencial linear, como std::randoferece a velocidade mais rápida para a qualidade mais baixa. Um gerador criptográfico oferece a mais alta qualidade com a menor velocidade.