Dlaczego rand ()% 6 jest obciążony?

Czytając, jak używać std :: rand, znalazłem ten kod na cppreference.com

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

Co jest nie tak z wyrażeniem po prawej stronie? Wypróbowałem i działa idealnie.

Istnieją dwa problemy z rand() % 6( 1+nie dotyczy żadnego problemu).

Po pierwsze, jak wskazało kilka odpowiedzi, jeśli niskie bity rand()nie są odpowiednio jednorodne, wynik operatora reszty również nie jest jednolity.

Po drugie, jeśli liczba odrębnych wartości utworzonych przez rand()nie jest wielokrotnością 6, to reszta da więcej wartości niskich niż wysokich. To prawda, nawet jeśli rand()zwraca idealnie rozłożone wartości.

Jako skrajny przykład udawaj, że rand()generuje równomiernie rozłożone wartości w zakresie [0..6]. Jeśli spojrzysz na reszty dla tych wartości, gdy rand()zwrócisz wartość z zakresu [0..5], reszta daje równomiernie rozłożone wyniki w zakresie [0..5]. Kiedy rand()zwraca 6, rand() % 6zwraca 0, tak jakby rand()zwróciło 0. W ten sposób otrzymujesz rozkład z dwukrotnie większą liczbą zer niż każda inna wartość.

Drugi to prawdziwy problem rand() % 6.

Sposobem na uniknięcie tego problemu jest odrzucenie wartości, które powodowałyby niejednorodne duplikaty. Obliczasz największą wielokrotność liczby 6, która jest mniejsza lub równa RAND_MAX, i za każdym razem, gdy rand()zwraca wartość, która jest większa lub równa tej wielokrotności, odrzucasz ją i ponownie wywołujesz `rand () tyle razy, ile potrzeba.

Więc:

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

To inna implementacja omawianego kodu, mająca na celu wyraźniejsze pokazanie, co się dzieje.

Są tu ukryte głębiny:

1. Zastosowanie małych plików uw RAND_MAX + 1u. RAND_MAXjest zdefiniowany jako inttyp i często jest największy z możliwych int. Zachowanie RAND_MAX + 1byłoby niezdefiniowane w takich przypadkach, gdy przepełnienie signedtypu. Pisanie 1uwymusza konwersję typu RAND_MAXna unsigned, zapobiegając w ten sposób przepełnieniu.

2. Użycie % 6 can (ale w każdej implementacji std::rand, którą widziałem , nie wprowadza żadnych dodatkowych błędów statystycznych poza przedstawioną alternatywą). Takie sytuacje, w których % 6jest niebezpieczny, to przypadki, w których generator liczb ma równiny korelacji w bitach niskiego rzędu, takie jak dość znana implementacja IBM (w języku C) z rand, jak sądzę, lat siedemdziesiątych XX wieku, która odwróciła wysokie i niskie bity jako „ostateczne zakrętas". Kolejną kwestią jest to, że 6 jest bardzo małe, por. RAND_MAX, więc efekt będzie minimalny, jeśli RAND_MAXnie będzie wielokrotnością liczby 6, co prawdopodobnie nie jest.

Podsumowując, ostatnio, ze względu na jego podatność, użyłbym % 6. Nie jest prawdopodobne, aby wprowadził jakiekolwiek anomalie statystyczne poza tymi, które wprowadza sam generator. Jeśli nadal masz wątpliwości, przetestuj swój generator, aby sprawdzić, czy ma on odpowiednie właściwości statystyczne dla Twojego przypadku użycia.

Ten przykładowy kod pokazuje, że std::randjest to przypadek legendarnego baldaszka kultowego cargo, który powinien podnosić brwi za każdym razem, gdy go widzisz.

Jest tu kilka problemów:

Ludzie kontraktowi zwykle zakładają - nawet biedne, nieszczęsne dusze, które nie wiedzą nic lepszego i nie będą myśleć o tym dokładnie w ten sposób - są takie, że randpróbki z jednolitego rozkładu liczb całkowitych w 0, 1, 2,… RAND_MAX,, a każde wywołanie daje niezależną próbkę.

Pierwszy problem polega na tym, że założona umowa, niezależne, jednolite losowe próbki w każdym zaproszeniu, nie jest w rzeczywistości tym, co mówi dokumentacja - aw praktyce historycznie wdrożenia nie zapewniły nawet najdrobniejszego symulakru niezależności. Na przykład C99 §7.20.2.1 „ randFunkcja” mówi bez rozwinięcia:

randFunkcja wylicza sekwencję liczb pseudo-losowych w przedziale od 0 do RAND_MAX.

To bezsensowne zdanie, ponieważ pseudolosowość jest właściwością funkcji (lub rodziny funkcji ), a nie liczby całkowitej, ale to nie powstrzymuje nawet biurokratów ISO przed nadużywaniem języka. W końcu jedyni czytelnicy, którzy byliby tym zdenerwowani, wiedzą, że lepiej nie czytać dokumentacji randze strachu przed rozpadem ich komórek mózgowych.

Typowa historyczna implementacja w C działa tak:

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

Ma to niefortunną właściwość, że nawet jeśli pojedyncza próbka może być równomiernie rozłożona w ramach jednolitego losowego ziarna (co zależy od określonej wartości RAND_MAX), naprzemiennie zmienia się między parzystymi i nieparzystymi liczbami całkowitymi w kolejnych wywołaniach - po

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

wyrażenie (a & 1) ^ (b & 1)daje 1 ze 100% prawdopodobieństwem, co nie ma miejsca w przypadku niezależnych prób losowych w dowolnym rozkładzie obsługiwanym przez parzyste i nieparzyste liczby całkowite. W ten sposób pojawił się kult cargo, w którym należy odrzucić mniej znaczące bity, aby ścigać nieuchwytną bestię o „lepszej losowości”. (Uwaga spoilera: to nie jest termin techniczny. To znak, że ktokolwiek czytasz prozę, albo nie wie, o czym mówi, albo myśli , że nie masz pojęcia i musi być protekcjonalny.)

Drugi problem polega na tym, że nawet gdyby każde wywołanie próbowało niezależnie od jednolitego losowego rozkładu na 0, 1, 2,… RAND_MAX,, wynik rand() % 6nie byłby równomiernie rozłożony na 0, 1, 2, 3, 4, 5 jak kostka rzuć, chyba że RAND_MAXjest przystająca do -1 modulo 6. Prosty kontrprzykład: Jeśli RAND_MAX= 6, to z rand(), wszystkie wyniki mają równe prawdopodobieństwo 1/7, ale z rand() % 6, wynik 0 ma prawdopodobieństwo 2/7, podczas gdy wszystkie inne wyniki mają prawdopodobieństwo 1/7 .

Właściwym sposobem na to jest próbkowanie odrzucenia: wielokrotnie losuj niezależną, jednolitą próbkę losową sz 0, 1, 2,… RAND_MAXi odrzucaj (na przykład) wyniki 0, 1, 2,… - ((RAND_MAX + 1) % 6) - 1jeśli otrzymasz jeden z te, zacznij od nowa; w przeciwnym razie wydajność s % 6.

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

W ten sposób zbiór wyników z rand(), który akceptujemy, jest podzielny po równo przez 6, a każdy możliwy wynik z s % 6jest uzyskiwany przez tę samą liczbę zaakceptowanych wyników z rand(), więc jeśli rand()jest równomiernie rozłożony, to tak jest s. Nie ma ograniczeń co do liczby prób, ale oczekiwana liczba jest mniejsza niż 2, a prawdopodobieństwo sukcesu rośnie wykładniczo wraz z liczbą prób.

Wybór , który Efekty o rand()odrzuceniu ma znaczenia, pod warunkiem, że mapa równą liczbę nich do każdej liczby całkowitej poniżej 6. Kod na cppreference.com sprawia, że inny wybór, bo od pierwszego problemu wyżej, że nic nie jest gwarantowane o dystrybucja lub niezależność wyników rand(), aw praktyce bity niskiego rzędu wykazywały wzorce, które nie „wyglądają wystarczająco losowo” (nie wspominając o tym, że następny wynik jest deterministyczną funkcją poprzedniego).

Ćwiczeń dla czytelnika: udowodnić, że kod w cppreference.com uzyskuje się równomierne rozprowadzenie na matrycy rolek, jeżeli rand()wydajność rozkład jednolity o 0, 1, 2, ..., RAND_MAX.

Ćwiczenie dla czytelnika: Dlaczego wolałbyś odrzucić jeden lub drugi podzbiór? Jakie obliczenia są potrzebne dla każdego procesu w dwóch przypadkach?

Trzeci problem polega na tym, że przestrzeń nasienna jest tak mała, że ​​nawet jeśli ziarno jest równomiernie rozłożone, przeciwnik uzbrojony w wiedzę o twoim programie i jednym wyniku, ale nie w ziarnie, może łatwo przewidzieć ziarno i późniejsze wyniki, co sprawia, że ​​nie wydają się takie w końcu losowe. Więc nawet nie myśl o używaniu tego do kryptografii.

Możesz przejść fantazyjną, nadmiernie inżynierską trasę i std::uniform_int_distributionklasę C ++ 11 z odpowiednim losowym urządzeniem i ulubionym losowym silnikiem, takim jak zawsze popularny twister Mersenne, std::mt19937aby grać w kości ze swoim czteroletnim kuzynem, ale nawet to nie będzie być zdolnym do generowania materiału klucza kryptograficznego - a twister Mersenne jest również strasznym świrem kosmicznym ze stanem wielokilobajtowym siejącym spustoszenie w pamięci podręcznej procesora z nieprzyzwoitym czasem konfiguracji, więc jest zły nawet dla np. równoległych symulacji Monte Carlo z odtwarzalne drzewa obliczeń podrzędnych; jego popularność prawdopodobnie wynika głównie z chwytliwej nazwy. Ale możesz go użyć do rzucania zabawkowymi kośćmi, jak na tym przykładzie!

Innym podejściem jest użycie prostego kryptograficznego generatora liczb pseudolosowych z małym stanem, na przykład prostego szybkiego usuwania klucza PRNG , lub po prostu szyfrowania strumieniowego, takiego jak AES-CTR lub ChaCha20, jeśli masz pewność ( np. W symulacji Monte Carlo dla badania w naukach przyrodniczych), że nie ma żadnych negatywnych konsekwencji w przewidywaniu przeszłych skutków, jeśli państwo kiedykolwiek zostanie zagrożone.

W żadnym wypadku nie jestem doświadczonym użytkownikiem C ++, ale chciałem sprawdzić, czy inne odpowiedzi mówią, std::rand()/((RAND_MAX + 1u)/6)że jestem mniej stronniczy, niż 1+std::rand()%6jest w rzeczywistości. Dlatego napisałem program testowy do zestawienia wyników dla obu metod (od dawna nie pisałem C ++, sprawdź to). Link do uruchomienia kodu znajduje się tutaj . Jest również odtwarzany w następujący sposób:

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

Następnie wziąłem wynik tego i użyłem chisq.testfunkcji w R, aby uruchomić test Chi-kwadrat, aby sprawdzić, czy wyniki różnią się znacznie od oczekiwanych. To pytanie dotyczące wymiany stosów dotyczy bardziej szczegółowo korzystania z testu chi-kwadrat do testowania uczciwości matrycy: Jak sprawdzić, czy kość jest uczciwa? . Oto wyniki kilku przebiegów:

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

W trzech przebiegach, które przeprowadziłem, wartość p dla obu metod była zawsze większa niż typowe wartości alfa używane do testowania istotności (0,05). Oznacza to, że nie uważalibyśmy żadnego z nich za stronnicze. Co ciekawe, rzekomo bezstronna metoda ma konsekwentnie niższe wartości p, co wskazuje, że w rzeczywistości może być bardziej stronnicza. Z zastrzeżeniem, że zrobiłem tylko 3 przebiegi.

AKTUALIZACJA: Kiedy pisałem swoją odpowiedź, Konrad Rudolph opublikował odpowiedź, która ma to samo podejście, ale daje zupełnie inny wynik. Nie mam reputacji, by komentować jego odpowiedź, więc odniosę się do tego tutaj. Po pierwsze, najważniejsze jest to, że kod, którego używa, używa tego samego ziarna dla generatora liczb losowych za każdym razem, gdy jest uruchamiany. Jeśli zmienisz ziarno, w rzeczywistości otrzymasz różnorodne wyniki. Po drugie, jeśli nie zmienisz ziarna, ale zmienisz liczbę prób, otrzymasz również różnorodne wyniki. Spróbuj zwiększyć lub zmniejszyć o rząd wielkości, aby zobaczyć, o co mi chodzi. Po trzecie, dochodzi do obcięcia lub zaokrąglenia liczb całkowitych, gdy oczekiwane wartości nie są całkiem dokładne. Prawdopodobnie to nie wystarczy, aby coś zmienić, ale jest.

Podsumowując, po prostu przypadkiem otrzymał właściwe ziarno i liczbę prób, które mogły uzyskać fałszywy wynik.

Można myśleć o generatorze liczb losowych jako o pracy na strumieniu cyfr binarnych. Generator zamienia strumień na liczby, dzieląc go na kawałki. Jeśli std:randfunkcja działa z RAND_MAXwartością 32767, to używa 15 bitów w każdym wycinku.

Kiedy weźmiemy moduły liczby od 0 do 32767 włącznie, znajdziemy 5462 '0 i' 1, ale tylko 5461 '2,' 3, '4 i' 5 '. Stąd wynik jest tendencyjny. Im większa wartość RAND_MAX, tym mniejsze będzie odchylenie, ale jest nieuniknione.

To, co nie jest odchylane, to liczba z zakresu [0 .. (2 ^ n) -1]. Możesz wygenerować (teoretycznie) lepszą liczbę z zakresu 0..5, wyodrębniając 3 bity, zamieniając je na liczbę całkowitą z zakresu 0..7 i odrzucając 6 i 7.

Można mieć nadzieję, że każdy bit w strumieniu bitów ma równe szanse na bycie „0” lub „1”, niezależnie od tego, gdzie znajduje się w strumieniu lub wartości innych bitów. W praktyce jest to wyjątkowo trudne. Wiele różnych implementacji programowych PRNG oferuje różne kompromisy między szybkością a jakością. Liniowy generator kongruencjalny, taki jak std::randoferuje największą prędkość przy najniższej jakości. Generator kryptograficzny zapewnia najwyższą jakość przy najniższej prędkości.