Jakość liniowych generatorów kongruencjalnych dla liczb losowych

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.$\langle \text{KE}\rangle=\frac{1}{2}k_BT$
• 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.$\langle x\rangle = 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ść.

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.

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$ $t$$d^t$$d^t$$I_1,\ldots$$\{I_{j+1}-I_j\}$ następuje w przybliżeniu znany rozkład.

• Zderzenie : generuj$n$ $t$wektory wymiarowe w $[0,1)^t$, dzieli je na $d^t$ równe hipersześciany i policz liczbę kolizji.

• MaxOft : generuj$n$ 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 $\mathbb{P}(X<x) = x^t$; wartości są$n=2\times 10^6$, $t=6$. Porównanie odbywa się za pomocą testu chi-kwadrat i testu Andersona-Darlinga:$\chi^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).