Generación de números aleatorios Introducción Existe una cierta variación en cuanto al término (números aleatorios o seudoaleatorios) a utilizar para referirse a las secuencias de números, obtenidos mediante algún método automático de generación digital, que se ajustan a una distribución uniforme en el intervalo [0,1). A partir de ahora utilizaremos el término de números aleatorios para referirnos a áquellos que cumplen lo anterior. En 1951 D.H. Lehmer1 , uno de los pioneros en computación y, especialmente, en teoría computacional de números, proporcionó la siguiente definición para números aleatorios: Una secuencia aleatoria es una noción difusa... en la que cada término es impredecible para el no iniciado y cuyos dígitos deben pasar una serie de pruebas tradicionales para los estadísticos... Aunque los diferentes métodos de generación de números aleatorios aparecen, casi siempre, descritos en libros de simulación, su utilización se extiende a muchas otras aplicaciones de las computadoras, como experimentos estadísticos, juegos de ordenador, criptografía, protocolos de seguridad, etc. Una secuencia de números aleatorios debe tener dos propiedades importantes, uniformidad e independencia. La función densidad de probabilidad debe cumplir: 

1 si 0 ≥ x ≥ 1 0 en el caso contrario R1 2 E(R) = 0 xdx = x2 |10 = 12

f (x) =

V (R) =

R1 0

x2 dx − (E(R))2 =

x3 1 3 |0

− ( 12 )2 =

1 3

−

1 4

=

1 12

Las consecuencias de la uniformidad e independencia referidas anteriormente son: 1 D. H. Lehmer. Mathematical methods in large-scale computing units. In Proc. 2nd Sympos. on Large-Scale Digital Computing Machinery, pages 141– 146. Harvard University Press, 1951.

1

- Si (0,1) se divide en n clases o subintervalos de igual longitud, el número esperado de observaciones en cada subintervalo es N/n, donde N es el número total de observaciones. - La probabilidad de observar un valor en un determinado subintervalo es independiente de las anteriores. Definición formal: Un generador de números aleatorios (nomenclatura inglesa: RNG) puede ser definido como una estructura G = (S, s0 , T, U, G), donde S es un conjunto conjunto finito de estados, s0 ∈ S es el estado inicial, T : S → S es la función de transición, U es un conjunto finito de símbolos de salida, y G : S → U es la función de salida. El generador comienza con el estado inicial s0 (llamado semilla) y evoluciona de acuerdo con si := T (si−1 ), para i ≥ 1 y en cada paso i, obtiene la observación ui := G(si ). Se espera que las observaciones se comporten como si fueran valores de variables aleatorias independientes e idénticamente distribuidas, uniformemente sobre U. El conjunto U es, a veces, un conjunto de enteros de la forma 0, . . . m − 1, o un conjunto finito de valores entre 0 y 1, representación aproximada de la distribución U[0,1). A partir de ahora se va a considerar este último caso. Puesto que S es finito, la secuencia de estados es, al final, periódica. El periodo es el entero positivo mas pequeño ρ tal que, dado algún entero τ ≥ 0 y para todo n ≥ τ , sρ+n = sn . El τ mas pequeño con esa propiedad se denomina transitorio. Cuando τ = 0 la secuencia se dice que es puramente periódica. Independientemente de los criterios de calidad estadística, para la elección de un generador de números aleatorios, también se consideran relevantes: la velocidad, las necesidades de memoria, la transportabilidad, la repetibilidad, la facilidad de implementación y la disponibilidad de saltos hacia adelante. En Kahaner2 se definen cinco características sobre las que se debería valorar cualquier generador de números aleatorios: Calidad. Se deben satisfacer los tests estadísticos adecuados. Antes de la repetición de la secuencia debiera haber un periodo lo suficientemente largo. Eficiencia. El generador debiera ser rápido. Debe precisar la menor cantidad de almacenamiento posible. Repetibilidad. El método debe depender de una semilla, para permitir que un experimento se pueda repetir. Transportabilidad. Si un método dado se implementa sobre un sitema diferente debe producir los mismos resultados. 2 D. Kahaner. Numerical Methods and Software, Prentice Hall Series in Computational Mathematics (1977)

2

Simplicidad. El generador debe ser suficientemente sencillo de implementar y utilizar. A continuación se detallan algunos de los métodos de generación, atendiendo a su actualidad más que a su desarrollo histórico.

Secuencias lineales recursivas. La mayoría de los RNG están basados en congruencias lineales de la forma: xn = (a1 xn−1 + . . . + ak xn−k ) mod m donde m se llama módulo, a1 , . . . , ak son enteros entre −m + 1 y m − 1 llamados multiplicadores con (ak 6= 0) y k es el orden de la recurrencia. Se puede definir el k . La longitud estado de la recurrencia en el paso n como sn = (xn , . . . , xn+k−1 ) ∈ Zm k del periodo máximo es m − 1.

Método de congruencias aditivas. Es un método rápido, puesto que no necesita realizar multiplicación. Se precisa una secuencia de números x1 , x2 . . . , xn . El generador produce una extensión de la secuencia xn+1 , xn+2 , . . . de la forma siguiente: xi = (xi−1 + xi−n ) mod m Por definición a = b mod m si a − b es divisible por m (resto 0). Por ejemplo, en módulo 4, los números 2, 6, 10, 14 son equivalentes porque (10 − 2), (10 − 6) . . . son todos divisibles por 4. Hay que tener en cuenta que, cuando utilizamos módulo m, los valores que resultarán estarán comprendidos entre 0 y m-1. Ejemplo: Sea una secuencia de enteros dados: x1 , x2 , x3 , x4 y x5 (57, 34, 89, 92 y 16), por tanto n = 5. Consideremos m = 100. La secuencia anterior puede ser ampliada por este método: x6 = (x5 + x1 ) mod 100 = 73 mod 100 = 73 x7 = (x6 + x2 ) mod 100 = 107 mod 100 = 7 x8 = (x7 + x3 ) mod 100 = 96 mod 100 = 96 ... Los números aleatorios se obtienen a partir de la relación : Ui−n = n + 1, n + 2, . . .. En el caso anterior los números serían:

Xi m

para i =

3

U1 = U2 = U3 =

X6 100 X7 100 X8 100

= 0,73 = 0,07 = 0,96

... En general estos métodos parten de la idea de la secuencia de Fibonacci xn que es generada de la forma siguiente: xn = xn−1 + xn−2 Los números obtenidos por la modificación de la secuencia : xn = (xn−1 + xn−2 ) mod m no tienen buenas propiedades aleatorias. Existe una extensión de esta idea de la forma: xn = (xn−5 + xn−17 ) mod 2k Marsaglia(1983) considera que este generador pasa la mayoría de los tests estadísticos. Utiliza, sin embargo, 17 posiciones de memoria con 17 enteros no todos impares. La suma ya es módulo 2k en las máquinas de k bits con complemento a 2. El periodo de este generador es sk (217 − 1). Para k=8, 16, 32 el periodo es 1,6 × 107 , 4,3 × 109 y 2,8 × 1014 respectivamente.

Generadores de congruencias lineales Una gran mayoría de los generadores utilizados actualmente utilizan esta técnica introducida por Lehmer en 1951. Una secuencia de números enteros Z1 , Z2 , . . . está definida por la fórmula recursiva: Zi = (aZi−1 + c) mod m donde el módulo m, el multiplicador a, el incremento c y la semilla o valor de comienzo Z0 son enteros no negativos. Evidentemente, por definición, se verifica que: 0 ≤ Zi < m Para obtener un número aleatorio de la distribución uniforme [0,1) se debe hacer Ui = Zi m . Además de ser no negativos se debe verificar que: 4

0 < m, a < m, c < m, Z0 < m Objeciones: No es una secuencia aleatoria pura, ya que: Z1 = (aZ0 + c) mod m Z2 = (aZ1 + c) mod m = (a2 Z0 + ca + c) mod m Z3 = (aZ2 + c) mod m = (a3 Z0 + ca2 + ca + c) mod m ... Por inducción: Zi = (ai Z0 + cai−1 + cai−2 + . . . + c) mod m i−2 + . . . + 1)) mod m Zi = h(ai Z0 + c(ai−1 + i a i −1) Zi = ai Z0 + c(aa−1 mod m

por lo que cualquier Zi está completamente determinado por m, a, c y Z0 . A pesar de todo, eligiendo adecuadamente los parámetros, se consigue que los Ui obtenidos aparezcan como pertenecientes a una U(0,1) cuando se realizan distintos tests sobre ellos. Otra objeción a realizar es que los números Ui solo pueden tomar valores racionales 0,9 0, por tanto la probabilidad de obtener un valor de Ui entre 0,1 m y m es 0,8 0 cuando debiera ser m > 0. Ahora bien, si se elige un m suficientemente grande los puntos en el intervalo [0,1) serán muy densos, ya que si elegimos un m ≥ 109 hay como mínimo mil millones de valores posibles para Ui espaciados todos la misma distancia. Se considera que esto es una aproximación lo suficentemente buena para la mayoría de los planteamientos. (m−1) 1 2 m, m, . . . m ,

Dada la forma de la expresión Zi es inevitable el comportamiento como un bucle, es decir, que en el momento que se repita un Zi todos los siguentes serán iguales a los obtenidos hasta ese momento. La longitud de cada uno de esos ciclos se conoce como periodo del generador y se representa por p. Como Zi solo depende de Zi−1 y se verifica que 0 ≤ Zi ≤ m − 1 es claro que p ≤ m. Si p = m el generador se llama de periodo total. Ejemplo: Sea m = 16, a = 5, c = 3 y Z0 = 7, la secuencia de los Zi obtenidos será: Zi = (5Zi−1 + 3) mod 16 5

i 0 1 2 3 4

Zi 7 6 1 8 11

Ui 0.375 0.063 0.500 0.688

i 5 6 7 8 9

Zi 10 5 12 15 14

Ui 0.625 0.313 0.750 0.938 0.875

i 10 11 12 13 14

Zi 9 0 3 2 13

Ui 0.563 0.000 0.188 0.125 0.813

i 15 16 17 18 19

Zi 4 7 6 1 8

Ui 0.250 0.438 0.375 0.063 0.500

Está claro que en simulación se necesitan generadores con periodos largos, fijando un m grande sería conveniente conseguir que tuviera periodo total. Teorema: El generador definido de la forma Zi = (aZi−1 + c) mod m tiene periodo total si y solo si se cumplen las siguientes condiciones: c y m son primos entre sí. Si q es un número primo que divide a m entonces q divide a (a-1) Si 4 divide a m entonces 4 divide a (a-1). Dependiendo del valor de c el comportamiento de los generadores puede ser distinto por lo que se distinguen generadores mixtos (c > 0) y generadores multiplicativos (c = 0). Ejemplos: IMSL utliza un generador multiplicativo con los siguientes parámetros:

a = 75 = 16,807 m = 231 − 1 = 2,147,483,647 c=0 Como m es un número primo, el periodo p, es igual a (m-1). Suponiendo que la semilla fuera X0 = 123,457 se obtendría :

X1 = 75 (123457) mod (231 − 1) = 2,074,941,799 X1 = 2,074,941,799 U1 =

X1 231

= 0,9662

X2 = 75 (2,074,941,799) mod (231 − 1) = 559,872,160 6

U2 =

X2 231

= 0,2607

NOTA.- La rutina de IMSL divide por (m+1) en lugar de hacerlo por m; sin embargo, para un valor tan grande como el de m, el efecto es despreciable. SIMSCRIPT II.5 utiliza un generador semejante con a = 630,360,016 Java Este es el generador utilizado para implementar el método nextDouble en la clase java.util.Random de la biblioteca estándar Java3 . Está basado en una recurrencia lineal con periodo de longitud 248 , pero cada valor de salida está construido tomando dos valores sucesivos de la recurrencia lineal, de la forma siguiente: xi+1 = (25214903917xi + 11) mod 248 ui = (227 bx2i /222 c + bx2i+1 / 221 c)/ 253 Conviene indicar que el generador rand48 de la biblioteca estándar Unix utiliza la misma recurrencia, pero produce su salida utilizando simplemente: ui = xi /248 . Visual Basic El generador utilizado en Microsoft Visual Basic4 es un generador de congruencias lineales (LCG) con un periodo de longitud 224 , definido por: xi = (1140671485xi−1 + 12820163) mod 224 ui = xi /224 Excel El generador implementado en Microsoft Excel5 es en esencia un LCG (generador de congruencias lineales), excepto que su recurrencia está dada por: ui = (9821,0ui−1 + 0,211327) mod 1 Está implementado directamente para los ui en la aritmética de punto flotante. Su longitud de periodo depende normalmente de la precisión de los números en coma flotante utilizados para su implementación. En la documentación no está perfectamente definido. A efectos de pruebas de validez no se utiliza el algoritmo sino que se genera una secuencia muy grande que se guarda en un archivo para su comprobación posterior. LCG16807 Este es un generador de congruencias lineales definido por: xi = 16807xi−1 mod (231 − 1) ui = xi /(231 − 1) 3

(http://java.sun.com/j2se/1.3/docs/api/java/util/Random.html) (http://support.microsoft.com/support/kb/articles/Q231/8/47.ASP) 5 (http://support.microsoft.com/directory) 4

7

con periodo de longitud 231 − 2, y propuesto originalmente por Lewis, Goodman, and Miller (1969). Este LCG se ha utilizado ampliamente en muchas bibliotecas de software para estadística, simulación, optimización, etc. así como en bibliotecas de sistemas operativos. Se ha recomendado en varios libros, por ejemplo, Bratley, Fox, y Schrage (1987) y en artículos, Law y Kelton (1982). Como curiosidad, este generador es utilizado en Arena y uno similar se utilizó en AutoMod (con el mismo módulo pero con el multiplicador 742938285) hasta hace poco, cuando los gestores de estos productos tuvieron la buena idea de reemplazarlo por MRG32k3a. También se utiliza en otros productos de sofware de simulación. MRG32k3a Este es el generador propuesto por L’Ecuyer (1999a). Combina dos generadores multiplicativos de orden 3 y su periodo tiene una longitud de aproximadamente 2191 . MT19937 Es el generador “trenzado” Mersene propuesto por Matsumoto y Nishimura6 (1998). Tiene el periodo de longitud más grande: 219937 − 1 y unas cualidades estadísticas muy buenas. Es uno de los generadores que oferta actualmente Matlab (versión 7.1).

Generadores de Tausworthe. Están relacionados con los métodos criptográficos, operan sobre los bits para formar números aleatorios. Se define una secuencia b1 , b2 , . . . de dígitos binarios por recurrencia. bi = (c1 bi−1 + c2 bi−2 + . . . + cq bi−q ) mod 2 donde c1 , c2 , . . . cq son contantes que pueden tomar el valor 0 o 1. El periodo máximo es 2q − 1. En casi todas las aplicaciones de los generadores de este tipo se suelen utilizar solo dos coeficientes ck distintos de 0, por tanto bi = (bi−r + bi−q ) mod 2 para enteros r y q que cumplen 0 < r < q. Hay que tener en cuenta que la suma de bits módulo 2 corresponde a la operación XOR.  bi =

0 si bi−r = bi−q 1 si bi−r 6= bi−q

6 Matsumoto, M. y T. Nishimura. Mersene twister: a 623- dimensionally equidistributed uniform pseudo-random number generator. ACM Transactions on Modeling and Computer Simulation 8 (1):3-30

8

Para inicializar la secuencia hay que dar valor a los primeros bi hasta q. Esto es equivalente a proporcionar la semilla Z0 en los geneneradores de congruencias. Ejemplo (Lewis y Payne 1973): Si r = 3 y q = 5 consideremos b1 , b2 , . . . b5 = 1 Para i ≥ 6 bi = (bi−3 + bi−5 ) mod 2 Los primeros 42 bits son : z }|{ z }| { z }| { z }| { z }|{ z }| { z }| { z }| { 1111 1 000 11 01 110 1 0100 0010 0 101 10 01 111 10001 10 . . . El periodo de los bits es 31 (2q − 1) La cuestión es: una vez obtenida la secuencia de bits bi , como se transforma dicha secuencia en números aleatorios U(0,1). Una posibilidad es juntar l consecutivos bi para formar un número de l bits entre 0 y 2l − 1 que será dividido por 2l . En el ejemplo anterior si se elige l = 4, se obtiene la secuencia: 15 8 13 13 4 2 5 9 15 1 , , , , , , , , , ,... 16 16 16 16 16 16 16 16 16 16 l debe ser como máximo la longitud de palabra del ordenador, también se pueden elegir secuencias de elementos dejando algunos bits sin utilizar cada vez. Ventajas: - Independientes del ordenador y de su tamaño de palabra. - Se pueden obtener secuencias de longitud considerable como 2521 − 1 > 10156 y aún mayores incluso en micros de 16 bits. Inconvenientes: - Aunque globalmente suelen presentar propiedades de aleatoriedad buenas, a veces, localmente no las tienen. - En general proporcionan malos resultados en las pruebas de rachas hacia arriba y hacia abajo. - Aunque la correlación de primer orden (un número con el siguiente) es casi cero se sospecha que algunos generadores pueden proporcionar valores elevados de correlaciones de orden elevado. - No todos los polinomios primitivos tiene las mismas cualidades. 9

Combinaciones de generadores. Si se combinan dos secuencias de números aleatorios Xi e Yi se produce una nueva secuencia Zi que reduce la no aleatoriedad. Si Xi e Yi están distribuidas sobre el rango 0 → (m − 1) se podría elegir una de estas opciones: 1.

Zi = (Xi + Yi ) mod m.

2.

Utilizar Yi para barajar Xi y hacer Zi igual a la secuencia barajada.

3.

Zi = Xi XOR Yi .

No hay garantía que una combinación del tipo 1 o 3 elimine la no aleatoriedad. Combinando un generador de congruencias con otro de Tausworthe se pueden atenuar muchas de las situaciones no deseadas que se producen. Wichmann y Hill7 informan de buenos resultados utilizando el método 1), es decir, combinando generadores (en este caso tres): Xi+1 = 171 Xi mod 30269 Yi+1 = 172 Yi mod 30307 Zi+1 = 170 Zi mod 30323 Por tanto Ui+1 =

Xi+1 Yi+1 Zi+1 + + 30269 30307 30323

El periodo es del orden de 1013 . Estos cálculos no precisan mas de 16 bits por lo que se considera un buen generador para ordenadores de 16 bits 8 . L’Ecuyer (1988)9 recomienda combinar: Xi+1 = 40014 × Xi mod 2147483563 Yi+1 = 40692 × Yi mod 2147483399 en la forma: Zi+1 = (Xi − Yi ) mod 2147483562 (

Zi+1 2,147,483,563 2,147,483,562 2,147,483,563

si Zi+1 ≥ 0 El periodo obtenido si Zi+1 = 0 es (m1 − 1) ∗ (m2 − 1)/2 ≈ 2 1018 , con m1 = 2147483563 y m2 = 2147483399. L’Ecuyer El número aleatorio devuelto es: Ui+1 =

7 Wichman, B.A. and I.D. Hill, Algorithm AS 183: An Efficient and Portable Pseudo-Random Number Generator, Applied Statistics, 31, 188-190, 1982. 8 En una nota publicada por Microsoft, indican que la versión del RANDU de 2003 implementa este algoritmo. http://support.microsoft.com/kb/828795 9 P. L’Ecuyer. Efficient and portable combined random number generators. Communications of the ACM, (31), 6, 742–751, 1988

10

y otros autores han planteado, como se ha indicado anteriormente, nuevas combinaciones de generadores para obtener periodos mucho más largos.

Consideraciones sobre la selección de la semilla. Normalmente la selección de la semilla para el comienzo de una secuencia de un generador no debe afectar a los resultados de la simulación. Si el generador es de periodo total y solo se precisa una variable, cualquier semilla es válida. Cuando se precisan distintas semillas por necesitar varias variables aleatorias diferentes en la misma simulación, se debe tener bastante cuidado. Casi todas las simulaciones usuales pertenecen a este tipo, por ejemplo, para simular una cola se necesita una llegada aleatoria y un tiempo de servicio aleatorio. Las normas que se dan a continuación son para este tipo de casos. Las dos primeras valen para el caso de una sola variable. 1.

No utilizar el cero Aunque puede ser una semilla válida para los generadores de congruencias lineales mixtos, llevará a cero a los generadores de congruencias lineales multiplicativos o a los de Tausworthe.

2.

Prohibir valores pares En teoría los valores pares son igual de buenos que los impares. De hecho para los generadores de periodo total, todos los valores de semilla no cero son igualmente buenos. Si un generador no es de periodo total, por ejemplo uno multiplicativo con módulo m = 2k , la semilla debe ser impar. Para impedir este tipo de situaciones por desconocimiento del tipo de generador es recomendable utilizar semillas impares.

3.

No subdividir una serie (stream) Utilizar una única serie para todas las variables es un error muy común. por ejemplo, si u1 , u2 , . . . es la secuencia generada por u0 , se podría pensar en utilizar u1 para generar los tiempos entre llegadas y u2 los tiempos de servicio y así sucesivamente u3 , u4 . . . Esto puede producir una correlación muy fuerte entre las variables.

4.

Utilizar series (stream) no solapadas Cada serie precisa una semilla separada. Si las semillas son tales que las series se solapan, existirá una correlación entre las series y las secuencias obtenidas no serán independientes. La idea es que si precisamos 10.000 números disponemos de una semilla para u0 y otra para u10,000 . En función de la recursividad de los números aleatorios se puede conocer cual es el número siguiente : xn = an x0 +

c(an −1) a−1

mod m 11

En algunos libros (Jain 91, pag 455) se proporcionan relaciones de semillas para secuencias de 100.000 en 100.000, para el generador xn = 75 xn−1 mod (231 − 1) 5.

Reutilizar semillas para sucesivas réplicas de la misma simulación

6.

No utilizar semillas aleatorias No se recomienda seleccionar semillas aleatorias. En particular no utilizar dos números aleatorios sucesivos, obtenidos del generador, como semillas. Utilizar, por ejemplo, la hora del día produce dos problemas: La simulación no puede ser reproducida. No es posible garantizar que series múltiples no se solapen entre si.

12

Referencias Las referencias básicas para este punto son: Jerry Banks, John S. Carson, Barry L. Nelson Discrete-Event System Simulation Prentice-Hall, 1996. James E. Gentle Random Number Generation and Monte Carlo Methods SpringerVerlag, 1998. Raj Jain The Art of Computer Systems Performance Analysis Techniques for Experimental Design, Measurement, Simulation, and Modeling John Wiley & Sons, 1992 Averill M. Law, W. David Kelton Simulation Modeling and Analysis McGraw-Hill, 1991. Otras referencias son: P. Bratley, B.L. Fox y L. Schrage A Guide to Simulation. 2a edición Springer-Verlag, 1987. George S. Fishman Discrete-Event Simulation. Tiene un tratamiento teórico muy completo respecto a la selección de los parámetros de los generadores. SpringerVerlag, 2001. Pierre L’Ecuyer Software for Uniform Random Number Generation: Distinguishing the good and the bad Winter Simulation Conference, 2001 En la dirección http://csrc.nist.gov/rng/rng6_4.html hay documentación y distintos tipos de pruebas para números aleatorios.

13