Qualité des générateurs congruentiels linéaires pour les nombres aléatoires

Je fais quelques simulations de l'équation de Langevin, pour diverses forces externes. Être dit que C est rand()de stdlib.hpeut introduire un biais dans mes résultats, je suis sur un Mersenne Twister.

Néanmoins, je voudrais savoir (et voir) exactement quel genre d'erreurs un générateur congruentiel linéaire peut introduire dans ma simulation. Ce sont des choses que j'ai essayées:

• Génération de tuples 3D de randoms pour essayer de voir les hyperplans. Je ne vois rien.
• Faire la FFT d'un grand vecteur de nombres aléatoires. C'est presque la même chose pour le Mersenne Twister et rand().
• Vérification du principe d'équipartition pour une particule en mouvement brownien. Les deux intégrateurs sont d' accord sur la valeur attendue de avec le même nombre de chiffres significatifs.$\langle \text{KE}\rangle=\frac{1}{2}k_BT$
• Voyant à quel point ils se rangent bien dans un certain nombre de bacs qui ne sont pas deux. Les deux donnent les mêmes résultats qualitatifs, personne n'est meilleur.
• En regardant les chemins browniens pour voir des divergences claires de . Encore une fois, pas de chance.$\langle x\rangle = 0$
• Répartition des points dans un cercle. Rempli, et uniquement dans le périmètre. Entre tous et entre voisins les plus proches (réponse de Shor, ci-dessous dans les commentaires). Disponible dans cet essentiel , exécutez-le simplement avec Julia 0.5.0 après avoir installé les bibliothèques nécessaires (voir l'essentiel pour les instructions).

Je voudrais souligner que je recherche des biais introduits dans le cadre de simulations physiques. Par exemple, j'ai vu comment rand()échoue lamentablement les tests de la perceuse alors que le Mersenne Twister ne le fait pas, mais pour le moment cela ne signifie pas trop pour moi.

Avez-vous des exemples physiques et concrets sur la façon dont un mauvais générateur de nombres aléatoires détruit une simulation Montecarlo?

Remarque: J'ai vu à quel point les PRNG RANDUpeuvent être horribles. Je m'intéresse à des exemples pas évidents, de générateurs qui semblent innocents mais qui introduisent finalement des biais.

Une référence intéressante qui décrit un échec d'une simulation Monte Carlo d'un système physique en raison d'un RNG inadéquat (bien qu'ils n'aient pas utilisé de LCG) est:

A. Ferrenberg et DP Landau. Monte Carlo Simulations: Erreurs cachées des "bons" générateurs de nombres aléatoires. Lettres d'examen physique 63 (23): 3382-3384, 1992.

Les modèles Ising que Ferrenberg et Landua ont étudiés sont de bons tests des RNG car vous pouvez comparer avec une solution exacte (pour le problème 2D) et trouver des erreurs dans les chiffres. Ces modèles devraient montrer les défauts d'un PMMLCG arithmétique 32 bits à l'ancienne sans trop de difficulté.

Une autre référence intéressante est:

H. Bauke et Stephan Mertens. Les pièces pseudo-aléatoires montrent plus de têtes que de queues. arXiv: cond-mat / 0307138 [cond-mat.stat-mech]

Bauke et Mertens plaident fortement contre les générateurs de nombres aléatoires de style à registre à décalage à rétroaction linéaire binaire. Bauke et Mertens ont d'autres articles à ce sujet.

Il peut être difficile de trouver les avions Marsaglia dans un nuage de points 3D. Vous pouvez essayer de faire pivoter l'intrigue pour obtenir une meilleure vue et parfois ils vous apparaîtront simplement. Vous pouvez également effectuer des tests 3D d'uniformité statistique - pour les LCG 32 bits plus anciens, ceux-ci échoueront sur un nombre de casiers assez petit. Par exemple, un test d'uniformité avec une grille de bacs de 20x20x20 en 3 dimensions est suffisant pour détecter le manque d'uniformité pour le PMMLCG largement utilisé (et précédemment bien considéré) avec m = 2 ^ 31-1, a = 7 ^ 5.

Il est possible d'utiliser la suite TestU01 de tests PRNG afin de savoir lequel de ces tests randéchoue. (Voir TestU01: Bibliothèque AC pour le test empirique des générateurs de nombres aléatoires pour un aperçu de la suite de tests.) C'est plus facile que de créer ses propres simulations Monte Carlo. D'une certaine manière, c'est aussi une question de composabilité logicielle (et d'exactitude logicielle): étant donné un PRNG qui semble fonctionner correctement sur de petits tests simples, comment savez-vous que ses comportements pathologiques ne seront pas déclenchés par un programme plus vaste?

Voici le code:

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

Pour la suite SmallCrush , il y a 3 tests qui échouent sur 15 (voir guidelongtestu01.pdf dans TestU01 pour de longues descriptions et toutes les références; ce sont 15 statistiques de 10 tests).

• $n$ $t$$d^t$$d^t$$I_1,\ldots$$\{I_{j+1}-I_j\}$

• $n$ $t$$[0,1)^t$$d^t$

• $n$$t$$[0,1)$$X$$n$$\mathbb{P}(X<x) = x^t$$n=2\times 10^6$$t=6$$\chi^2$$<10^{-300}$

En supposant que ce sont toutes des simulations Monte Carlo "typiques" (même si elles ne ressemblent pas aux problèmes que vous aviez en tête), la conclusion est que randcertains sous-ensembles inconnus échouent. Je ne sais pas pourquoi c'est spécifiquement ce sous-ensemble, cependant, il m'est impossible de dire si cela fonctionnera sur votre propre problème ou non.

MaxOft semble particulièrement suspect, étant donné la simplicité de la description.

Parmi les tests de la suite Crush , randéchoue 51 sur 140 (140 statistiques en 96 tests). Certains des tests ayant échoué (comme Fourier3 ) sont effectués sur des chaînes de bits, il est donc possible qu'ils ne soient pas pertinents pour vous. Un autre test curieux qui échoue est GCD , qui teste la distribution du GCD de deux entiers aléatoires. (Encore une fois, je ne sais pas pourquoi ce test particulier échoue ou si votre simulation en souffrira.)

PS : Encore une autre chose à noter est que rand()c'est en fait plus lent que certains PRNG qui réussissent tous les tests SmallCrush , Crush , BigCrush , tels que MRG32k3a (voir le document L'Ecuyer & Simard ci-dessus).