Generowanie losowej liczby całkowitej z zakresu

Potrzebuję funkcji, która wygeneruje losową liczbę całkowitą w podanym zakresie (w tym wartości graniczne). Nie mam nieuzasadnionych wymagań dotyczących jakości / losowości, mam cztery wymagania:

• Potrzebuję tego, żeby był szybki. Mój projekt musi generować miliony (a czasem nawet dziesiątki milionów) liczb losowych, a moja obecna funkcja generatora okazała się wąskim gardłem.
• Potrzebuję, aby był w miarę jednolity (użycie rand () jest całkowicie w porządku).
• zakresy min-max mogą wynosić od <0, 1> do <-32727, 32727>.
• musi być zaszczepiany.

Obecnie mam następujący kod w C ++:

output = min + (rand() * (int)(max - min) / RAND_MAX)

Problem w tym, że nie jest tak naprawdę jednolity - max jest zwracany tylko wtedy, gdy rand () = RAND_MAX (dla Visual C ++ jest to 1/32727). Jest to poważny problem w przypadku małych zakresów, takich jak <-1, 1>, gdzie ostatnia wartość prawie nigdy nie jest zwracana.

Więc złapałem długopis i papier i wymyśliłem następującą formułę (która opiera się na sztuczce zaokrąglania liczb całkowitych (int) (n + 0,5)):

Ale to nadal nie daje mi jednolitej dystrybucji. Powtarzane serie z 10000 próbek dają mi stosunek 37:50:13 dla wartości -1, 0,1.

Czy mógłbyś zaproponować lepszą formułę? (lub nawet cała funkcja generatora liczb pseudolosowych)

Jest to szybkie, nieco lepsze niż twoje, ale wciąż niejednorodne rozwiązanie rozproszone

output = min + (rand() % static_cast<int>(max - min + 1))

Z wyjątkiem sytuacji, gdy rozmiar zakresu jest potęgą 2, ta metoda daje tendencyjne niejednorodne rozłożone liczby niezależnie od jakości rand(). Aby uzyskać kompleksowy test jakości tej metody, przeczytaj to .

Najprostszą (a przez to najlepszą) odpowiedzią w C ++ (przy użyciu standardu 2011) jest

#include <random>

std::random_device rd;     // only used once to initialise (seed) engine
std::mt19937 rng(rd());    // random-number engine used (Mersenne-Twister in this case)
std::uniform_int_distribution<int> uni(min,max); // guaranteed unbiased

auto random_integer = uni(rng);

Nie ma potrzeby ponownego wynajdywania koła. Nie musisz się martwić o uprzedzenia. Nie musisz martwić się o wykorzystanie czasu jako losowego ziarna.

Jeśli Twój kompilator obsługuje C ++ 0x i używanie go jest opcją dla Ciebie, to nowy standardowy <random>nagłówek prawdopodobnie spełni Twoje potrzeby. Ma wysoką jakość, uniform_int_distributionktóra akceptuje minimalne i maksymalne granice (włącznie z potrzebami) i możesz wybierać spośród różnych generatorów liczb losowych, aby podłączyć się do tej dystrybucji.

Oto kod, który generuje milion losowych liczb intrównomiernie rozmieszczonych w [-57, 365]. Użyłem nowych funkcji standardowego <chrono>czasu, ponieważ wspomniałeś, że wydajność jest dla ciebie głównym problemem.

#include <iostream>
#include <random>
#include <chrono>

int main()
{
    typedef std::chrono::high_resolution_clock Clock;
    typedef std::chrono::duration<double> sec;
    Clock::time_point t0 = Clock::now();
    const int N = 10000000;
    typedef std::minstd_rand G;
    G g;
    typedef std::uniform_int_distribution<> D;
    D d(-57, 365);
    int c = 0;
    for (int i = 0; i < N; ++i) 
        c += d(g);
    Clock::time_point t1 = Clock::now();
    std::cout << N/sec(t1-t0).count() << " random numbers per second.\n";
    return c;
}

Dla mnie (2,8 GHz Intel Core i5) to wypisuje:

2.10268e + 07 liczb losowych na sekundę.

Możesz zaszczepić generator, przekazując int do jego konstruktora:

    G g(seed);

Jeśli później okaże się, że intnie obejmuje zakresu potrzebnego do dystrybucji, można temu zaradzić, zmieniając coś uniform_int_distributionpodobnego (np. Na long long):

    typedef std::uniform_int_distribution<long long> D;

Jeśli później okaże się, że minstd_randgenerator nie jest wystarczająco wysokiej jakości, można go łatwo wymienić. Na przykład:

    typedef std::mt19937 G;  // Now using mersenne_twister_engine

Posiadanie oddzielnej kontroli nad generatorem liczb losowych i rozkładem losowym może być dość wyzwalające.

Obliczyłem również (nie pokazano) pierwsze 4 „momenty” tego rozkładu (używając minstd_rand) i porównałem je z wartościami teoretycznymi, próbując określić ilościowo jakość rozkładu:

min = -57
max = 365
mean = 154.131
x_mean = 154
var = 14931.9
x_var = 14910.7
skew = -0.00197375
x_skew = 0
kurtosis = -1.20129
x_kurtosis = -1.20001

( x_Przedrostek oznacza „oczekiwany”)

Podzielmy problem na dwie części:

• Wygeneruj liczbę losową nz zakresu od 0 do (maks-min).
• Dodaj min do tej liczby

Pierwsza część jest oczywiście najtrudniejsza. Załóżmy, że wartość zwracana przez rand () jest idealnie jednolita. Użycie modulo doda odchylenie do pierwszych (RAND_MAX + 1) % (max-min+1)liczb. Więc gdybyśmy mogli magicznie zmienić się RAND_MAXna RAND_MAX - (RAND_MAX + 1) % (max-min+1), nie byłoby już żadnych uprzedzeń.

Okazuje się, że możemy skorzystać z tej intuicji, jeśli jesteśmy skłonni dopuścić pseudo-niedeterminizm do czasu działania naszego algorytmu. Za każdym razem, gdy rand () zwraca zbyt dużą liczbę, po prostu prosimy o inną liczbę losową, aż otrzymamy wystarczająco małą.

Czas działania jest teraz rozłożony geometrycznie , z wartością oczekiwaną, 1/pgdzie pjest prawdopodobieństwo uzyskania wystarczająco małej liczby przy pierwszej próbie. Ponieważ RAND_MAX - (RAND_MAX + 1) % (max-min+1)jest zawsze mniejsze niż (RAND_MAX + 1) / 2, wiemy o tym p > 1/2, więc oczekiwana liczba iteracji zawsze będzie mniejsza niż dwa dla dowolnego zakresu. Powinno być możliwe wygenerowanie dziesiątek milionów liczb losowych w mniej niż sekundę na standardowym procesorze z tą techniką.

EDYTOWAĆ:

Chociaż powyższe jest technicznie poprawne, odpowiedź DSimona jest prawdopodobnie bardziej przydatna w praktyce. Nie powinieneś sam wdrażać tego. Widziałem wiele implementacji próbkowania odrzucania i często bardzo trudno jest sprawdzić, czy jest poprawne, czy nie.

Co powiesz na Mersenne Twister ? Implementacja przyspieszenia jest raczej łatwa w użyciu i dobrze przetestowana w wielu rzeczywistych aplikacjach. Sam używałem go w kilku akademickich projektach, takich jak sztuczna inteligencja i algorytmy ewolucyjne.

Oto ich przykład, w którym wykonują prostą funkcję rzucania sześciościenną kostką:

#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_int.hpp>
#include <boost/random/variate_generator.hpp>

boost::mt19937 gen;

int roll_die() {
    boost::uniform_int<> dist(1, 6);
    boost::variate_generator<boost::mt19937&, boost::uniform_int<> > die(gen, dist);
    return die();
}

Aha, a tutaj jeszcze trochę stręczycielstwa tego generatora, na wypadek gdybyś nie był przekonany, że powinieneś go używać na znacznie gorszym rand():

Mersenne Twister to generator „liczb losowych” wymyślony przez Makoto Matsumoto i Takuji Nishimurę; ich strona internetowa zawiera liczne implementacje algorytmu.

Zasadniczo Mersenne Twister to bardzo duży rejestr przesuwny z liniowym sprzężeniem zwrotnym. Algorytm działa na podstawie 19 937 bitów przechowywanych w 624-elementowej tablicy 32-bitowych liczb całkowitych bez znaku. Wartość 2 ^ 19937-1 jest liczbą pierwszą Mersenne'a; technika manipulowania ziarnem oparta jest na starszym algorytmie „skręcania” - stąd nazwa „Mersenne Twister”.

Atrakcyjnym aspektem Mersenne Twister jest wykorzystanie operacji binarnych - w przeciwieństwie do czasochłonnego mnożenia - do generowania liczb. Algorytm ma również bardzo długi okres i dobrą ziarnistość. Jest zarówno szybki, jak i skuteczny w aplikacjach niekryptograficznych.

int RandU(int nMin, int nMax)
{
    return nMin + (int)((double)rand() / (RAND_MAX+1) * (nMax-nMin+1));
}

Jest to mapowanie 32768 liczb całkowitych na (nMax-nMin + 1) liczb całkowitych. Mapowanie będzie całkiem dobre, jeśli (nMax-nMin + 1) jest małe (jak w twoim wymaganiu). Zauważ jednak, że jeśli (nMax-nMin + 1) jest duże, mapowanie nie zadziała (na przykład - nie możesz odwzorować 32768 wartości na 30000 wartości z równym prawdopodobieństwem). Jeśli takie zakresy są potrzebne - powinieneś użyć 32-bitowego lub 64-bitowego losowego źródła zamiast 15-bitowego rand () lub zignorować wyniki rand (), które są poza zakresem.

Oto bezstronna wersja, która generuje liczby w [low, high]:

int r;
do {
  r = rand();
} while (r < ((unsigned int)(RAND_MAX) + 1) % (high + 1 - low));
return r % (high + 1 - low) + low;

Jeśli twój zakres jest dość mały, nie ma powodu, aby buforować prawą stronę porównania w dopętli.

Polecam bibliotekę Boost.Random , jest bardzo szczegółowa i dobrze udokumentowana, pozwala wyraźnie określić, jaką dystrybucję chcesz, aw scenariuszach niekryptograficznych może faktycznie przewyższać typową implementację rand biblioteki C.

załóżmy, że min i max są wartościami int, [i] oznacza dołączenie tej wartości, (i) oznacza, że ​​nie należy uwzględniać tej wartości, używając powyższego do uzyskania właściwej wartości za pomocą c ++ rand ()

odniesienie: for () [] zdefiniuj, odwiedź:

https://en.wikipedia.org/wiki/Interval_(mathematics)

dla funkcji rand i srand lub zdefiniuj RAND_MAX odwiedź:

http://en.cppreference.com/w/cpp/numeric/random/rand

[minimum maksimum]

int randNum = rand() % (max - min + 1) + min

(minimum maksimum]

int randNum = rand() % (max - min) + min + 1

[minimum maksimum)

int randNum = rand() % (max - min) + min

(minimum maksimum)

int randNum = rand() % (max - min - 1) + min + 1

W tym wątku próbkowanie odrzucania było już omówione, ale chciałem zasugerować jedną optymalizację opartą na tym fakcie rand() % 2^something nie wprowadza żadnego błędu, jak już wspomniano powyżej.

Algorytm jest naprawdę prosty:

• obliczyć najmniejszą potęgę 2 większą niż długość interwału
• losuj jedną liczbę w tym „nowym” przedziale
• zwraca tę liczbę, jeśli jest mniejsza niż długość pierwotnego interwału
  • odrzucić inaczej

Oto mój przykładowy kod:

int randInInterval(int min, int max) {
    int intervalLen = max - min + 1;
    //now calculate the smallest power of 2 that is >= than `intervalLen`
    int ceilingPowerOf2 = pow(2, ceil(log2(intervalLen)));

    int randomNumber = rand() % ceilingPowerOf2; //this is "as uniform as rand()"

    if (randomNumber < intervalLen)
        return min + randomNumber;      //ok!
    return randInInterval(min, max);    //reject sample and try again
} 

Działa to dobrze zwłaszcza w przypadku małych interwałów, ponieważ potęga 2 będzie „bliżej” rzeczywistej długości interwału, a więc liczba chybionych będzie mniejsza.

PS
Oczywiście unikanie rekurencji byłoby bardziej wydajne (nie ma potrzeby obliczania ponad pułapem dziennika ..), ale pomyślałem, że w tym przykładzie będzie bardziej czytelny.

Zauważ, że w większości sugestii początkowa losowa wartość, którą otrzymałeś z funkcji rand (), która zwykle wynosi od 0 do RAND_MAX, jest po prostu marnowana. Tworzysz z tego tylko jedną liczbę losową, podczas gdy istnieje rozsądna procedura, która może dać ci więcej.

Załóżmy, że chcesz [min, max] region liczb całkowitych losowych. Zaczynamy od [0, max-min]

Weź podstawę b = max-min + 1

Zacznij od przedstawienia liczby otrzymanej z rand () w bazie b.

W ten sposób masz piętro (log (b, RAND_MAX)), ponieważ każda cyfra w bazie b, z wyjątkiem ostatniej, reprezentuje losową liczbę z zakresu [0, max-min].

Oczywiście końcowe przesunięcie do [min, max] jest proste dla każdej liczby losowej r + min.

int n = NUM_DIGIT-1;
while(n >= 0)
{
    r[n] = res % b;
    res -= r[n];
    res /= b;
    n--;
}

Jeśli NUM_DIGIT jest liczbą cyfr w podstawie b, którą można wyodrębnić, to znaczy

NUM_DIGIT = floor(log(b,RAND_MAX))

to powyższe jest prostą implementacją wyodrębnienia NUM_DIGIT liczb losowych od 0 do b-1 z jednej liczby losowej RAND_MAX, zapewniając b <RAND_MAX.

Wzór na to jest bardzo prosty, więc wypróbuj to wyrażenie,

 int num = (int) rand() % (max - min) + min;  
 //Where rand() returns a random number between 0.0 and 1.0

Poniższe wyrażenie powinno być bezstronne, jeśli się nie mylę:

std::floor( ( max - min + 1.0 ) * rand() ) + min;

Zakładam tutaj, że rand () daje losową wartość z zakresu od 0,0 do 1,0 NIE wliczając 1,0 i że max i min to liczby całkowite z warunkiem, że min <max.