2.1 Generadores congruenciales lineales

En los generadores congruenciales lineales se considera una combinación lineal de los últimos \(k\) enteros generados y se calcula su resto al dividir por un entero fijo \(m\). En el método congruencial simple (de orden \(k = 1\)), partiendo de una semilla inicial \(x_0\), el algoritmo secuencial es el siguiente: \[\begin{aligned} x_{i} & = (ax_{i-1}+c) \bmod m \\ u_{i} & = \dfrac{x_{i}}{m} \\ i & =1,2,\ldots \end{aligned}\] donde \(a\) (multiplicador), \(c\) (incremento) y \(m\) (módulo) son enteros positivos8 fijados de antemano (los parámetros de este generador). Si \(c=0\) el generador se denomina congruencial multiplicativo (Lehmer, 1951) y en caso contrario se dice que es mixto (Rotenberg, 1960).

Obviamente los parámetros y la semilla determinan los valores generados, que también se podrían obtener de forma no recursiva: \[x_{i}=\left( a^{i}x_0+c\frac{a^{i}-1}{a-1}\right) \bmod m\]

Este método está implementado (no de la forma más eficiente9) en la función rlcg() del paquete simres (fichero rng.R):

library(simres)
rlcg
## function(n, seed = as.numeric(Sys.time()), a = 7^5, c = 0, m = 2^31 - 1) {
##   u <- numeric(n)
##   for(i in 1:n) {
##     seed <- (a * seed + c) %% m
##     u[i] <- seed/m # (seed + 1)/(m + 1)
##   }
##   # Almacenar semilla y parámetros
##   assign(".rng", list(seed = seed, type = "lcg",
##           parameters = list(a = a, c = c, m = m)), envir = globalenv())
##   # .rng <<- list(seed = seed, type = "lcg", parameters = list(a = a, c = c, m = m))
##   # Para continuar con semilla y parámetros:
##   #   with(.rng, rlcg(n, seed, parameters$a, parameters$c, parameters$m))
##   # Devolver valores
##   return(u)
## }
## <bytecode: 0x0000015a8a846d40>
## <environment: namespace:simres>

Tratando de imitar el funcionamiento del generador uniforme de R, la semilla se almacena en el objeto .rng del entorno de trabajo y se genera a partir de la hora del sistema si no se establece. No obstante, en lugar de emplear esta función, puede resultar más cómodo utilizar set.rng() para establecer la semilla (y opcionalmente los parámetros de este generador, que ya se selecciona por defecto con type = "lcg") y posteriormente rng() para generar múltiples secuencias.

set.rng(543210)
u <- rng(1); u
## [1] 0.25136
u == rlcg(1, seed = 543210)
## [1] TRUE

A pesar de su simplicidad, una adecuada elección de los parámetros permite obtener de manera eficiente secuencias de números “aparentemente” i.i.d. (independientes e idénticamente distribuidos) \(\mathcal{U}(0,1)\). Algunos ejemplos de parámetros, empleados en la actualidad o en el pasado, son:

  • \(c=0\), \(a=7^{5}=16807\) y \(m=2^{31}-1\) (primo de Mersenne), Park y Miller (1988) minimal standar, empleado por las librerías IMSL y NAG.

  • \(c=0\), \(a=48271\) y \(m=2^{31}-1\) actualización del minimal standar propuesta por Park et al. (1993).

  • \(c=0\), \(a=2^{16}+3=65539\) y \(m=2^{31}\), generador RANDU de IBM (no recomendable como veremos más adelante).

Durante los primeros años, el procedimiento habitual consistía en escoger \(m\) de forma que se pudiera realizar eficientemente la operación del módulo, aprovechando la arquitectura del ordenador (por ejemplo \(m = 2^{31}\) si se emplean enteros con signo de 32 bits). Posteriormente se seleccionaban \(c\) y \(a\) de forma que el período \(p\) fuese lo más largo posible (o suficientemente largo), empleando los resultados mostrados a continuación.

Teorema 2.1 (Hull y Dobell, 1962)
Un generador congruencial tiene período máximo \(p = m\) si y solo si:

  1. \(c\) y \(m\) son primos relativos (i.e. \(m.c.d.(c, m) = 1\)).

  2. \(a-1\) es múltiplo de todos los factores primos de \(m\) (i.e. \(a \equiv 1 \bmod q\), para todo \(q\) factor primo de \(m\)).

  3. Si \(m\) es múltiplo de \(4\), entonces \(a-1\) también lo ha de ser (i.e. \(m \equiv 0 \bmod 4\Rightarrow a \equiv 1 \bmod 4\)).

Algunas consecuencias:

  • Un generador multiplicativo no cumple la condición 1 (\(m.c.d.(0, m) = m\)).

  • Si \(m\) primo y \(c \neq 0\), \(p = m\) si y solo si \(a = 1\).

Teorema 2.2 (Gauss, 1801)
Un generador multiplicativo tiene período máximo \(p = m-1\) si:

  1. \(m\) es primo.

  2. \(a\) es una raíz primitiva de \(m\) (i.e. \(a \neq 0\) y \(a^{(m-1)/k} \not\equiv 1 \bmod m\) para todo factor primo \(k\) de \(m - 1\)).

Sin embargo, además de preocuparse de la longitud del ciclo, sería mucho más importante que las secuencias generadas se comporten de forma similar a muestras aleatorias simple de una \(\mathcal{U}(0,1)\). Uno de los principales problemas con este tipo de generadores (y con muchos otros) es que los valores generados pueden mostrar una estructura reticular. Este es el caso del generador RANDU de IBM muy empleado en la década de los 70 y que muestra una clara estructura reticular cuando se consideran más de dos dimensiones (aunque ya se detectó este problema en 1963 se siguió utilizando hasta los 90). Para ilustrar este problema podríamos emplear el conjunto de datos randu del paquete base datasets que contiene 400 tripletas de números sucesivos generados con la implementación de VAX/VMS 1.5 (de 1977, y que no se corrigió hasta la versión 2.0 de 1980). Aunque también podemos emplear el siguiente código10 (ver Figura 2.1):

u <- rlcg(n = 9999, seed = 543210, a = 2^16 + 3, c = 0, m = 2^31)
library(plot3D)
xyz <- stats::embed(u, 3)
points3D(xyz[,3], xyz[,2], xyz[,1], colvar = NULL, phi = 60, 
         theta = -50, pch = 21, cex = 0.2)
Grafico de dispersión de tripletas del generador RANDU de IBM (contenidas en 15 planos).

Figura 2.1: Grafico de dispersión de tripletas del generador RANDU de IBM (contenidas en 15 planos).

En general todos los generadores de este tipo van a presentar estructuras reticulares. Marsaglia (1968) demostró que las k-uplas de un generadores multiplicativo están contenidas en a lo sumo \(\left(k!m\right)^{1/k}\) hiperplanos paralelos (2344 como máximo con \(k=3\) y \(m=2^{31}\)). Por tanto habría que seleccionar adecuadamente los parámetros del generador congruencial de forma que la estructura reticular sea imperceptible, teniendo en cuenta el número de datos que se pretende generar (por ejemplo de forma que la distancia mínima entre los puntos sea próxima a la esperada en teoría). Para más detalles sobre la estructura reticular ver por ejemplo Ripley (1987, sec. 2.7).

Se han propuesto diversas pruebas (ver Sección 2.3) para determinar si un generador tiene problemas de este tipo y se han realizado numerosos estudios para determinadas familias (e.g. Park y Miller, 1988, estudiaron los multiplicadores adecuados para \(m=2^{31}-1\), analizando el comportamiento de sus raíces primitivas.). En ciertos contextos muy exigentes (por ejemplo en criptografía), se recomienda considerar un “periodo de seguridad” \(\approx \sqrt{p}\) para evitar este tipo de problemas.

Aunque estos generadores tienen limitaciones en su capacidad para producir secuencias muy largas de números i.i.d. \(\mathcal{U}(0,1)\), son un elemento básico en generadores más avanzados (incluyendo el empleado por defecto en R) como veremos en la siguiente sección.

Ejemplo 2.1

Consideramos el generador congruencial, de ciclo máximo, definido por: \[\begin{aligned} x_{n+1} & =(5x_{n}+1)\ \bmod\ 512,\nonumber\\ u_{n+1} & =\frac{x_{n+1}}{512},\ n=0,1,\dots\nonumber \end{aligned}\]

Generamos 500 valores de este generador, obteniendo de paso el tiempo de CPU requerido:

set.rng(321, "lcg", a = 5, c = 1, m = 512)  # Establecer semilla y parámetros
nsim <- 500
system.time(u <- rng(nsim))
##    user  system elapsed 
##       0       0       0

Representamos la distribución de los valores generados mediante un histograma, en escala de densidades, de forma que podemos compararla con la densidad teórica:

hist(u, freq = FALSE)
abline(h = 1, col = "blue") # Densidad uniforme
Histograma de los valores generados.

Figura 2.2: Histograma de los valores generados.

En este caso concreto la distribución de los valores generados es aparentemente más uniforme de lo que cabría esperar, lo que induciría a sospechar de la calidad de este generador (ver Ejemplo 2.2 en Sección 2.3).

A partir de los valores generados podríamos aproximar características de la distribución teórica. Por ejemplo, la aproximación por simulación la media teórica sería:

mean(u)
## [1] 0.49996

Como se puede observar, esta aproximación es muy cercana al valor teórico 0.5 (el error absoluto es \(3.90625\times 10^{-5}\)).

Como ejemplo adicional, la aproximación (mediante simulación) de la probabilidad del intervalo \((0.4, 0.8)\) sería:

sum((0.4 < u) & (u < 0.8))/nsim
## [1] 0.402
mean((0.4 < u) & (u < 0.8))     # Alternativa
## [1] 0.402

Si la comparamos con la probabilidad teórica (\(0.8 - 0.4 = 0.4\)), también observamos que es muy próxima, lo que podría aumentar las sospechas de que este generador no reproduce adecuadamente la variabilidad de una distribución \(\mathcal{U}(0,1)\) (en el Ejemplo 2.2 se realizará un análisis más riguroso de la calidad de este generador).

Bibliografía

Bratley, P., Fox, B. L., y Schrage, L. E. (1983). A Guide to Simulation. Springer. https://link.springer.com/book/10.1007/978-1-4684-0167-7
Lehmer, D. H. (1951). Mathematical models in large-scale computing units. The Annals of the Computation Laboratory of Harvard University, 26, 141-146.
Marsaglia, G. (1968). Random numbers fall mainly in the planes. Proceedings of the National Academy of Sciences, 61(1), 25-28.
Park, S. K., y Miller, K. W. (1988). Random number generators: good ones are hard to find. Communications of the ACM, 31(10), 1192-1201. https://doi.org/10.1145/63039.63042
Park, S. K., Miller, K. W., y Stockmeyer, P. K. (1993). Technical correspondence: Response. Communications of the ACM, 36(7), 108-110.
Ripley, B. D. (1987). Stochastic Simulation. Wiley. https://www.wiley.com/en-us/Stochastic+Simulation-p-9780470009604
Rotenberg, A. (1960). A new pseudo-random number generator. Journal of the ACM, 7(1), 75-77.
Schrage, L. (1979). A more portable Fortran random number generator. ACM Transactions on Mathematical Software (TOMS), 5(2), 132-138.

  1. Se supone además que \(a\), \(c\) y \(x_0\) son menores que \(m\), ya que, dadas las propiedades algebraicas de la suma y el producto en el conjunto de clases de resto módulo \(m\) (que es un anillo), cualquier otra elección de valores mayores o iguales que \(m\) tiene un equivalente verificando esta restricción.↩︎

  2. Para evitar problemas computacionales se recomienda realizar el cálculo de los valores empleando el método de Schrage (1979; ver e.g. Bratley et al., 1983, sec. 6.5.2).↩︎

  3. La función stats::embed(u, 3) devuelve una matriz en la que la fila i-ésima contiene los valores u[i+2], u[i+1] y u[i]. Además, en lugar de la función plot3D::points3D() se podría utilizar la función plot3d() del paquete rgl, y rotar la figura (pulsando con el ratón) para ver los hiperplanos rgl::plot3d(xyz).↩︎