Qualidade de geradores congruenciais lineares para números aleatórios

Estou fazendo algumas simulações da equação de Langevin, para várias forças externas. Sendo dito que C's rand()destdlib.h podem introduzir viés nos meus resultados, eu estou usando um Mersenne Twister.

No entanto, eu gostaria de saber (e ver) exatamente que tipo de erros um gerador congruencial linear pode introduzir em minha simulação. Estas são as coisas que eu tentei:

• Gerando tuplas 3D de randoms para tentar ver hiperplanos. Eu não consigo ver nada.
• Fazendo a FFT de um grande vetor de números aleatórios. É quase o mesmo para o Mersenne Twister erand() .
• Verificando o princípio da equipartição para uma partícula em movimento browniano. Ambos os integradores de concordar no valor esperado de $\langle \text{KE}\rangle=\frac{1}{2}k_BT$ com o mesmo número de dígitos significativos.
• Vendo o quão bem eles se encaixam em um número de caixas que não é um poder dois. Ambos dão os mesmos resultados qualitativos, ninguém sendo melhor.
• Olhando para caminhos browniano para ver divergências claras de $\langle x\rangle = 0$ . Mais uma vez, sem sorte.
• Distribuição de pontos em um círculo. Preenchido, e apenas no perímetro. Entre todos eles e entre os vizinhos mais próximos (resposta de Shor, abaixo nos comentários). Disponível nesta lista , basta executá-la com Julia 0.5.0 após instalar as bibliotecas necessárias (consulte a lista para obter instruções).

Gostaria de enfatizar que estou procurando viés introduzido no contexto de simulações físicas. Por exemplo, eu vi comorand() falha miseravelmente os testes do dieharder, enquanto o Mersenne Twister não, mas no momento isso não significa muito para mim.

Você tem exemplos físicos, concretos, de como um gerador de números aleatórios ruim destrói uma simulação de Montecarlo?

Nota: Vi como o PRNG RANDUpode ser horrível. Estou interessado em exemplos não óbvios, de geradores que parecem inocentes, mas que acabam apresentando viés.

Uma referência interessante que descreve uma falha de uma simulação de Monte Carlo de um sistema físico devido ao RNG inadequado (embora eles não usassem um LCG) é:

A. Ferrenberg e DP Landau. Simulações de Monte Carlo: erros ocultos de geradores de números aleatórios "bons". Physical Review Letters 63 (23): 3382-3384, 1992.

Os modelos de Ising que Ferrenberg e Landua estudaram são bons testes de RNGs, porque você pode comparar com uma solução exata (para o problema 2-D) e encontrar erros nos dígitos. Esses modelos devem mostrar as falhas em um PMMLCG aritmético de 32 bits à moda antiga sem muita dificuldade.

Outra referência interessante é:

H. Bauke e Stephan Mertens. Moedas pseudo-aleatórias mostram mais cabeças do que caudas. arXiv: cond-mat / 0307138 [cond-mat.stat-mech]

Bauke e Mertens defendem fortemente os geradores de números aleatórios no estilo do registro de deslocamento de feedback linear binário. Bauke e Mertens têm alguns outros documentos relacionados a isso.

Pode ser difícil encontrar os aviões Marsaglia em um gráfico de dispersão 3D. Você pode tentar girar a plotagem para ter uma visão melhor e, às vezes, elas aparecem em você. Você também pode fazer testes 3D de uniformidade estatística - para LCGs de 32 bits mais antigos, eles falharão em um número bastante pequeno de posições. por exemplo, um teste de uniformidade com uma grade de caixas de 20x20x20 em 3 dimensões é suficiente para detectar falta de uniformidade para o PMMLCG amplamente utilizado (e anteriormente bem considerado) com m = 2 ^ 31-1, a = 7 ^ 5.

É possível usar o conjunto TestU01 de testes PRNG para descobrir qual desses testesrand falha. (Consulte TestU01: Biblioteca AC para testes empíricos de geradores de números aleatórios para obter uma visão geral da suíte de testes.) Isso é mais fácil do que criar simulações de Monte Carlo. De certa forma, também é uma questão de composição de software (e correção de software): dado um PRNG que parece funcionar bem em testes pequenos e simples, como você sabe que seus comportamentos patológicos não serão acionados por um programa maior?

Aqui está o código:

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

Para o pacote SmallCrush , há três testes com falha em 15 (consulte guidelongtestu01.pdf no TestU01 para obter descrições longas e todas as referências; essas são 15 estatísticas de 10 testes).

• $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$ máximos com a distribuição teórica $\mathbb{P}(X<x) = x^t$; os valores são$n=2\times 10^6$, $t=6$. A comparação é feita através do teste do qui-quadrado e do teste de Anderson-Darling: o$\chi^2$ teste falha com valor-p $<10^{-300}$, enquanto o teste do AD passa (esse teste do AD falha em alguns dos testes similares maiores no Crush).

Supondo que sejam todas simulações "típicas" de Monte Carlo (embora possam não ser como os problemas que você tinha em mente), a conclusão é que randfalha algum subconjunto desconhecido deles. Não sei por que é especificamente esse subconjunto, portanto, é impossível para mim dizer se funcionará no seu próprio problema ou não.

O MaxOft parece particularmente suspeito, dado o quão direta é a descrição.

Entre os testes no pacote Crush , randfalha 51 de 140 (140 estatísticas em 96 testes). Alguns dos testes com falha (como Fourier3 ) são feitos em cadeias de bits, então talvez seja possível que não sejam relevantes para você. Outro teste curioso que falha é o GCD , que testa a distribuição do GCD de dois números inteiros aleatórios. (Novamente, não sei por que esse teste específico falha ou se a sua simulação sofrerá com isso.)

PS : Outro aspecto a ser observado é que, rand()na verdade, é mais lento do que alguns PRNGs que são aprovados com êxito em todos os testes de SmallCrush , Crush , BigCrush , como o MRG32k3a (consulte o artigo da L'Ecuyer & Simard acima).