Jakość liniowych generatorów kongruencjalnych dla liczb losowych

14

Wykonuję symulacje równania Langevina dla różnych sił zewnętrznych. Powiedziano mi, że C rand()z stdlib.hmoże wprowadzić błąd w moich wynikach, używam Twistera Mersenne.

Niemniej jednak chciałbym wiedzieć (i zobaczyć) dokładnie, jakie błędy liniowy kongruencjalny generator może wprowadzić w mojej symulacji. Oto rzeczy, które próbowałem:

  • Generowanie krotek losowych 3D, aby spróbować zobaczyć hiperpłaszczyzny. Nic nie widzę.
  • Wykonanie FFT dużego wektora liczb losowych. Jest prawie tak samo dla Mersenne Twister i rand().
  • Sprawdzanie zasady ekwipartycji dla cząstki w ruchu Browna. Oba integratorów zgodni w oczekiwanej wartości z taką samą liczbą cyfr znaczących.KE=12kBT
  • Widząc, jak dobrze kosztują w kilku pojemnikach, które nie są potęgą drugą. Oba dają takie same wyniki jakościowe, nikt nie jest lepszy.
  • Patrząc na Browna ścieżek, aby zobaczyć wyraźne rozbieżności z . Znowu nie ma szczęścia.x=0
  • Rozkład punktów w okręgu. Wypełnione i tylko na obwodzie. Pomiędzy nimi wszystkimi i między najbliższymi sąsiadami (odpowiedź Shora poniżej w komentarzach). Dostępne w tej liście, po prostu uruchom ją z Julią 0.5.0 po zainstalowaniu potrzebnych bibliotek (instrukcje znajdziesz w liście).

Chciałbym podkreślić, że szukam wprowadzonego uprzedzenia w kontekście symulacji fizycznych. Widziałem na przykład, jak rand()żałośnie nie udaje się zagorzałych testów, podczas gdy Mersenne Twister nie, ale w tej chwili nie znaczy to dla mnie zbyt wiele.

Czy masz jakieś fizyczne, konkretne przykłady tego, jak zły generator liczb losowych niszczy symulację Montecarlo?

Uwaga: Widziałem, jak PRNG RANDUmoże być okropny. Interesują mnie nieoczywiste przykłady generatorów, które wyglądają niewinnie, ale ostatecznie wprowadzają stronniczość.

RedPointyJackson
źródło
1
Nie mam żądanych przykładów, ale używam drand48 () / srand48 () zamiast rand () / srand () w moich własnych programach C. Ich odpowiednie strony podręcznika man dokumentują różne używane algorytmy prng (patrz algorytm rand losowy) i uważam, że drand48 jest ogólnie preferowany, chociaż moje szczegółowe zrozumienie jest znikomo małe. Kiedy chcę gwarantowaną przenośną odtwarzalność na różnych platformach, kodowałem ran1 () z Numerical Recipes in C, 2nd Edition, WHPress i in., Cambridge UP 1992, ISBN 0-521-43108-5, strona 280. Działa świetnie, o ile Mogę powiedzieć, ale nie testowałem ilościowo.
Użyj random () lub drand48 () / lrand48 () (zawsze używam tego drugiego do symulacji molekularnych i symulacji Monte Carlo i jest całkiem niezły). Spróbuj także użyć losowego ziarna. To powinno wystarczyć do symulacji równania Langevina dla pojedynczej cząstki.
valerio
Użyliśmy obwodu, a nie koła.
@PeterShor Dziękujemy za korektę. Zaktualizowałem odpowiedź, obawiam się, że nadal nie mam szczęścia.
RedPointyJackson
1
@DanielShapero random i urandom mają być kryptograficznie bezpiecznymi RNG, przeznaczonymi do celów kryptograficznych, takich jak generowanie kluczy. Sprzętu aspektem jest to, że na Linuksie, używają entropię środowiska, to nie jest taka sama jak z akceleracją sprzętową. W rzeczywistości nie są przeznaczone do niczego takiego jak symulacje Monte Carlo.
Kirill,

Odpowiedzi:

3

Ciekawym odniesieniem opisującym niepowodzenie symulacji Monte Carlo systemu fizycznego z powodu nieodpowiedniego RNG (chociaż nie używali LCG) jest:

A. Ferrenberg i DP Landau. Symulacje Monte Carlo: ukryte błędy w „dobrych” generatorach liczb losowych. Physical Review Letters 63 (23): 3382-3384, 1992.

Modele Isinga, które badali Ferrenberg i Landua, są dobrymi testami RNG, ponieważ można je porównać z dokładnym rozwiązaniem (dla problemu 2-D) i znaleźć błędy w cyfrach. Modele te powinny bezbłędnie pokazywać usterki w starej 32-bitowej arytmetyce PMMLCG.

Innym interesującym odniesieniem jest:

H. Bauke i Stephan Mertens. Pseudolosowe monety pokazują więcej głów niż ogonów. arXiv: cond-mat / 0307138 [cond-mat.stat-mech]

Bauke i Mertens zdecydowanie popierają generatory liczb losowych w stylu rejestru przesuwnego z binarnym sprzężeniem zwrotnym. Bauke i Mertens mają na ten temat kilka innych dokumentów.

Znalezienie płaszczyzn Marsaglii na wykresie punktowym 3D może być trudne. Możesz spróbować obrócić fabułę, aby uzyskać lepszy widok, a czasem wyskakują na ciebie. Możesz również przeprowadzić testy 3D o jednolitości statystycznej - w przypadku starszych 32-bitowych LCG zawiodą przy dość małej liczbie pojemników. np. test jednorodności z siatką pojemników 20 x 20 x 20 w 3 wymiarach jest wystarczający do wykrycia braku jednolitości dla powszechnie stosowanego (i wcześniej dobrze ocenianego) PMMLCG o m = 2 ^ 31-1, a = 7 ^ 5.

Brian Borchers
źródło
1

Aby to sprawdzić, można użyć zestawu testów PRNG TestU01 który z tych testów randkończy się niepowodzeniem. (Patrz TestU01: Biblioteka AC do testowania empirycznego generatorów liczb losowych aby się z zestawem testów). To łatwiejsze niż własnych symulacji Monte Carlo. W pewnym sensie jest to również kwestia kompozycyjności oprogramowania (i poprawności oprogramowania): biorąc pod uwagę PRNG, który wydaje się działać dobrze na małych, prostych testach, skąd wiesz, że jego patologiczne zachowania nie będą wyzwalane przez większy program?

Oto kod:

#include "TestU01.h"

int main() {
  // Same as rand() on my machine
  unif01_Gen* gen = ulcg_CreateLCG(2147483647, 16807, 0, 12345);

  bbattery_SmallCrush(gen);
  bbattery_Crush(gen);

  return 0;
}

W przypadku pakietu SmallCrush istnieją 3 testy zakończone niepowodzeniem na 15 (patrz guidelongtestu01.pdf w TestU01 dla długich opisów i wszystkich odniesień; są to 15 statystyk z 10 testów).

  • n tretretja1,{jajot+1-jajot} następuje w przybliżeniu znany rozkład.

  • Zderzenie : generujn twektory wymiarowe w [0,1)t, dzieli je na ret równe hipersześciany i policz liczbę kolizji.

  • MaxOft : generujn grupy t wartości w [0,1), obliczyć maksimum X w każdej grupie i porównaj rozkład n maksima z rozkładem teoretycznym P.(X<x)=xt; wartości sąn=2)×106, t=6. Porównanie odbywa się za pomocą testu chi-kwadrat i testu Andersona-Darlinga:χ2) test kończy się niepowodzeniem z wartością p <10-300, podczas gdy test AD przechodzi pomyślnie (ten test AD kończy się niepowodzeniem w niektórych większych podobnych testach w Crush).

Zakładając, że są to wszystkie „typowe” symulacje Monte Carlo (choć mogą nie być podobne do problemów, o których myślałeś), wniosek jest taki, że randzawodzi jakiś nieznany ich podzbiór. Nie wiem jednak, dlaczego jest to konkretnie ten podzbiór, więc nie mogę powiedzieć, czy zadziała na twoim problemie, czy nie.

MaxOft wydaje się szczególnie podejrzliwy, biorąc pod uwagę, jak prosty jest opis.

Wśród testów w pakiecie Crushrand nie udaje się 51 ze 140 (140 statystyk w 96 testach). Niektóre z nieudanych testów (jak Fourier3 ) są wykonywane na ciągach bitów, więc być może możliwe, że nie będą one dla ciebie odpowiednie. Innym ciekawym testem, który się nie powiedzie, jest GCD , który testuje rozkład GCD dwóch losowych liczb całkowitych. (Znów nie wiem, dlaczego ten konkretny test kończy się niepowodzeniem ani czy cierpi z tego powodu twoja symulacja).

PS : Kolejną rzeczą, na którą należy zwrócić uwagę, jest to, że rand()faktycznie jest wolniejsze niż niektóre PRNG, które pomyślnie przeszły wszystkie testy SmallCrush , Crush , BigCrush , takie jak MRG32k3a (patrz artykuł L'Ecuyer i Simard powyżej).

Cyryl
źródło