R: Problema con runif: el número generado se repite (más a menudo de lo esperado) después de menos de 100 000 pasos

Question

Después de ejecutar el código

RNGkind(kind="Mersenne-Twister")  # the default anyway
set.seed(123)
n = 10^5
x = runif(n)
print(x[22662] == x[97974])

TRUE ¡es la salida!

Si utilizo, por ejemplo, RNGkind(kind="Knuth-TAOCP-2002") ocurre lo mismo: Obtengo "sólo" 99 995 valores diferentes en x . Dados los periodos de ambos generadores aleatorios, los resultados parecen muy poco probables.

¿Estoy haciendo algo mal? Necesito generar al menos un millón de números aleatorios.

Estoy usando Windows 8.1 con la versión 3.6.2 de R; Plataforma: x86_64-w64-mingw32/x64 (64-bit) y RStudio 1.2.5033.

---

Resultados adicionales:

1. Tener una bolsa con $n$ diferentes bolas, elegimos una bola $m$ veces y lo devuelve cada vez. La probabilidad $p_{n, m}$ que todas las bolas elegidas son diferentes es igual a ${n\choose m} / (n^m m!)$ .
2. La documentación de R apunta a un enlace donde está disponible la implementación de Mersenne-Twister para máquinas de 64 bits: http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt64.html

El muestreo uniforme de $[0, 1]$ se obtiene eligiendo primero un entero aleatorio de 64 bits, por lo que calculé las probabilidades anteriores para los 64 bits y (cuando $p_{64, 10^5}$ resultó ser bastante bajo) caso de 32 bits: $$ p_{2^{64}, 10^5}\doteq 0.9999999999972... \qquad p_{2^{32}, 10^5} \doteq 0.3121... $$

A continuación, he probado con 1000 semillas aleatorias y he calculado la proporción de casos en los que todos los números generados son diferentes: 0,303.

Así que, actualmente, asumo que por alguna razón, se utilizan enteros de 32 bits.

Answer 1

La documentación de R sobre la generación de números aleatorios tiene algunas frases al final, que confirman tu expectativa de que se usen enteros de 32 bits y podrían explicar lo que estás observando:

No confíe en la aleatoriedad de los bits de orden inferior de los RNG. La mayoría de los generadores uniformes suministrados devuelven valores enteros de 32 bits que se convierten en dobles, por lo que toman como máximo 2^32 valores distintos y las ejecuciones largas devolverán valores duplicados (Wichmann-Hill es la excepción, y todos dan al menos 30 bits variables).

Así que la implementación en R parece ser diferente a lo que se explica en la página web de los autores del Mersenne Twister. Posiblemente combinando esto con la paradoja del cumpleaños, se esperarían duplicados con sólo 2^16 números con una probabilidad de 0,5, y 10^5 > 2^16. Probando el algoritmo de Wichmann-Hill como se sugiere en la documentación:

RNGkind(kind="Wichmann-Hill") 
set.seed(123)
n = 10^8
x = runif(n)
length(unique(x))
# 1e8

Tenga en cuenta que el generador de números aleatorios original de Wichmann-Hill tiene la propiedad de que su siguiente número puede predecirse por el anterior, y por lo tanto no cumple los requisitos de no predictibilidad de un PRNG válido. Véase este documento de Dutang y Wuertz, 2009 (sección 3)

Answer 2

Sólo para enfatizar la aritmética del $2^{32}$ en términos del número de valores potenciales distintos: si se muestrea $10^5$ tiempos de $2^{32}$ valores con reemplazo, se esperaría una media de $2^{32}\left(1-\left(1-\frac{1}{2^{32}}\right)^{10^5}\right) \approx 10^5 - 1.1634$ valores distintos, observando que $\frac{(10^5)^2}{2 \times 2^{32}} \approx 1.1642$ se acerca a este déficit

Por lo tanto, es de esperar que haya muchos ejemplos anteriores. Hay dos pares con set.seed(1) :

n = 10^5
set.seed(1)
x = runif(n)
x[21101] == x[56190]
x[33322] == x[50637]

Si se hace algo similar a la primera $2000$ semillas en R para el valor por defecto runif se obtiene una media de $10^5 - 1.169$ valores únicos, lo que se acerca a la expectativa calculada. Sólo $30.8\%$ de estas semillas no producen duplicados de la muestra de $10^5$

Muestra $10^6$ veces y se esperaría que el problema fuera unas cien veces peor y, de hecho, el número medio de valores únicos para la primera $2000$ semillas es $10^6 - 116.602$ y ninguna de estas semillas no produjo duplicados

Hay otra forma de reducir la probabilidad de solapamientos sin dejar de tener una distribución uniforme: pruebe con pnorm(rnorm(n))

  set.seed(123)
  n = 10^8
  x = runif(n) 
  length(unique(x))
  # 98845390
  y = pnorm(rnorm(n))
  length(unique(y))
  # 100000000

Answer 3

Aunque sea contraintuitivo, hay buenas razones que explican este fenómeno, esencialmente porque un ordenador utiliza una precisión finita. Se acaba de publicar un preprint (marzo de 2020) en ArXiv (como ya se ha mencionado en la discusión, por cierto) y trata esta cuestión a fondo. Ha sido escrito por un investigador experimentado en estadística computacional (no soy yo ni un amigo mío) y utiliza R. Todos los códigos son reproducibles y puedes comprobar fácilmente los códigos y las afirmaciones por ti mismo. Sólo citaré algunas líneas (primeras líneas de la Conclusión) de la conclusión que parecen responder a tu pregunta:

De forma poco intuitiva (pero, como resulta, no inesperada), la generación de números aleatorios puede dar lugar a empates. Para generar $n$ números aleatorios en un $k$ -de bits, demostramos que el número esperado de empates es $n-2^{k}(1-(1-2^{-k})^{n})$ . Además, derivamos una fórmula numéricamente robusta para calcular este número. Para una arquitectura de 32 bits como la que se sigue utilizando en los generadores de números aleatorios (ya sea por razones históricas, de reproducibilidad o de tiempo de ejecución), el número esperado de empates al generar un millón de números aleatorios es de 116.

La versión citada es la publicada el 18 de marzo de 2020.

https://arxiv.org/abs/2003.08009

Answer 4

Aquí hay dos problemas. El primero ha sido bien cubierto en las otras respuestas, a saber: por qué aparecen duplicados para ciertas configuraciones de los argumentos de entrada.

La otra es muy importante: hay una gran diferencia entre "aleatorio con reemplazo" y "aleatorio permutación de un conjunto conocido. "Matemáticamente, es completamente válido que la secuencia de enteros aleatorios contenga, por ejemplo, 6,6,6,6 . La mayoría de los PRNG no hacen un "reemplazo" completo en su algoritmo, por lo que lo que acabamos obteniendo está mucho más cerca (pero no exactamente, como muestra el ejemplo de la pregunta publicada) de una permutación aleatoria del conjunto de valores. De hecho, como la mayoría de los PRNG generan el siguiente valor basándose en el valor actual (y posiblemente en algunos anteriores), son casi procesos de Markov. Los llamamos "aleatorios" porque estamos de acuerdo en que un observador externo no puede determinar el algoritmo del generador, por lo que el siguiente número que aparece es imprevisible para ese observador.

Considere, entonces, la diferencia entre runif y sample donde este último tiene un argumento que indica explícitamente si se selecciona con o sin reemplazo.