Generación de números aleatorios

Introducción

Existe una cierta variación en cuanto al término (números aleatorios o seudoaleato-rios) a utilizar para referirse a las secuencias de números, obtenidos mediante algúnmétodo automático de generación digital, que se ajustan a una distribución uniformeen el intervalo [0,1). A partir de ahora utilizaremos el término de números aleatoriospara referirnos a áquellos que cumplen lo anterior.

Aunque los diferentes métodos de generación de números aleatorios aparecen, casisiempre, descritos en libros de simulación, su utilización se extiende a muchas otrasaplicaciones de las computadoras, como experimentos estadísticos, juegos de orde-nador, 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:

f(x) ={

1 si 0 ≥ x ≥ 10 en el caso contrario

E(R) =∫ 10 xdx = x2

2 |10 = 1

2

V (R) =∫ 10 x2dx− (E(R))2 = x3

3 |10 − (1

2)2 = 13 −

14 = 1

12

Las consecuencias de la uniformidad e independencia referidas anteriormente son:

- Si (0,1) se divide en n clases o subintervalos de igual longitud, el número esper-ado de observaciones en cada subintervalo es N/n, donde N es el número totalde observaciones.

- La probabilidad de observar un valor en un determinado subintervalo es inde-pendiente de las anteriores.

Definición formal:

1

Un generador de números aleatorios (nomenclatura inglesa: RNG) puede ser defi-nido como una estructura G = (S, s0, T, U, G), donde S es un conjunto conjunto finitode estados, s0 ∈ S es el estado inicial, T : S → S es la función de transición, U es unconjunto finito de símbolos de salida, y G : S → U es la función de salida. El gener-ador comienza con el estado inicial s0 (llamado semilla) y evoluciona de acuerdo consi := T (si−1), para i ≥ 1 y en cada paso i, obtiene la observación ui := G(si). Se esperaque las observaciones se comporten como si fueran valores de variables aleatorias in-dependientes 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 valoresentre 0 y 1, representación aproximada de la distribución U[0,1). A partir de ahora seva a considerar este último caso. Puesto que S es finito, la secuencia de estados es, alfinal, periódica. El periodo es el entero positivo mas pequeño ρ tal que, dado algúnentero τ ≥ 0 y para todo n ≥ τ , sρ+n = sn. El τ mas pequeño con esa propiedad sedenomina 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 ungenerador de números aleatorios, también se consideran relevantes: la velocidad, lasnecesidades de memoria, la transportabilidad, la repetibilidad, la facilidad de imple-mentación y la disponibilidad de saltos hacia adelante.

A continuación se detallan algunos de los métodos de generación, atendiendo a suactualidad 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 = (a1xn−1 + . . . + akxn−k) mod m

donde m se llama módulo, a1, . . . , ak son enteros entre −m + 1 y m − 1 llamadosmultiplicadores con (ak 6= 0) y k es el orden de la recurrencia. Se puede definir elestado de la recurrencia en el paso n como sn = (xn, . . . , xn+k−1) ∈ Zk

m. La longituddel periodo máximo es mk − 1.

Método de congruencias aditivas.

Es un método rápido, puesto que no necesita realizar multiplicación. Se precisauna secuencia de números x1, x2 . . . , xn. El generador produce una extensión de lasecuencia xn+1, xn+2, . . . de la forma siguiente:

xi = (xi−1 + xi−n) mod m

2

Por definición a = bmodm si a−b es divisible por m (resto 0). Por ejemplo, en módulo4, los números 2, 6, 10, 14 son equivalentes porque (10− 2), (10− 6) . . . son todos di-visibles por 4. Hay que tener en cuenta que, cuando utilizamos módulo m, los valoresque 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 = 73x7 = (x6 + x2) mod 100 = 107 mod 100 = 7x8 = (x7 + x3) mod 100 = 96 mod 100 = 96

. . .

Los números aleatorios se obtienen a partir de la relación : Ui−n = Xim para i =

n + 1, n + 2, . . .. En el caso anterior los números serían:

U1 = X6100 = 0,73

U2 = X7100 = 0,07

U3 = X8100 = 0,96. . .

En general estos métodos parten de la idea de la secuencia de Fibonacci xn que esgenerada 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 peri-odo 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.

3

Generadores de congruencias lineales

Una gran mayoría de los generadores utilizados actualmente utilizan esta técnicaintroducida 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 comien-zo 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 =Zim . Además de ser no negativos se debe verificar que:

0 < m, a < m, c < m,Z0 < m

Objeciones: No es una secuencia aleatoria pura, ya que:

Z1 = (aZ0 + c) mod mZ2 = (aZ1 + c) mod m = (a2Z0 + ca + c) mod mZ3 = (aZ2 + c) mod m = (a3Z0 + ca2 + ca + c) mod m. . .

Por inducción:

Zi = (aiZ0 + cai−1 + cai−2 + . . . + c) mod mZi = (aiZ0 + c(ai−1 + ai−2 + . . . + 1)) mod m

Zi =[aiZ0 + c(ai−1)

a−1

]mod m

por lo que cualquier Zi está completamente determinado por m,a, c y Z0. A pesarde todo, eligiendo adecuadamente los parámetros, se consigue que los Ui obtenidosaparezcan como pertenecientes a una U(0,1) cuando se realizan distintos tests sobreellos.

Otra objeción a realizar es que los números Ui solo pueden tomar valores racionales0, 1

m , 2m , . . . (m−1)

m , por tanto la probabilidad de obtener un valor de Ui entre 0,1m y 0,9

m es0 cuando debiera ser 0,8

m > 0. Ahora bien, si se elige un m suficientemente grandelos puntos en el intervalo [0,1) serán muy densos, ya que si elegimos un m ≥ 109

4

hay como mínimo mil millones de valores posibles para Ui espaciados todos la mismadistancia. Se considera que esto es una aproximación lo suficentemente buena para lamayoría de los planteamientos.

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 alos obtenidos hasta ese momento. La longitud de cada uno de esos ciclos se conocecomo periodo del generador y se representa por p. Como Zi solo depende de Zi−1 yse verifica que 0 ≤ Zi ≤ m− 1 es claro que p ≤ m. Si p = m el generador se llama deperiodo 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

i Zi Ui i Zi Ui i Zi Ui i Zi Ui

0 7 - 5 10 0.625 10 9 0.563 15 4 0.2501 6 0.375 6 5 0.313 11 0 0.000 16 7 0.4382 1 0.063 7 12 0.750 12 3 0.188 17 6 0.3753 8 0.500 8 15 0.938 13 2 0.125 18 1 0.0634 11 0.688 9 14 0.875 14 13 0.813 19 8 0.500

Está claro que en simulación se necesitan generadores con periodos largos, fijandoun 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 solosi 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 distintopor lo que se distinguen generadores mixtos (c > 0) y generadores multiplicativos(c = 0).

Ejemplos:

IMSL utliza un generador multiplicativo con los siguientes parámetros:

5

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 semillafuera X0 = 123,457 se obtendría :

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

X1 = 2,074,941,799

U1 = X1231 = 0,9662

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

U2 = X2231 = 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 laclase java.util.Random de la biblioteca estándar Java1. Está basado en una recurrencialineal con periodo de longitud 248, pero cada valor de salida está construido tomandodos valores sucesivos de la recurrencia lineal, de la forma siguiente:

xi+1 = (25214903917xi + 11) mod 248

ui = (227bx2i/222c+ bx2i+1/ 221c)/ 253

Conviene indicar que el generador rand48 de la biblioteca estándar Unix utiliza lamisma recurrencia, pero produce su salida utilizando simplemente: ui = xi/248.

Visual Basic El generador utilizado en Microsoft Visual Basic2 es un generador decongruencias lineales (LCG) con un periodo de longitud 224, definido por:

xi = (1140671485xi−1 + 12820163) mod 224

ui = xi/224

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

6

Excel El generador implementado en Microsoft Excel3 es en esencia un LCG (gen-erador 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. Sulongitud de periodo depende normalmente de la precisión de los números en comaflotante utilizados para su implementación. En la documentación no está perfecta-mente definido. A efectos de pruebas de validez no se utiliza el algoritmo sino que segenera una secuencia muy grande que se guarda en un archivo para su comprobaciónposterior.

LCG16807 Este es un generador de congruencias lineales definido por:

xi = 16807xi−1 mod (231 − 1)ui = xi/(231 − 1)

con periodo de longitud 231 − 2, y propuesto originalmente por Lewis, Goodman,and Miller (1969). Este LCG se ha utilizado ampliamente en muchas bibliotecas desoftware para estadística, simulación, optimización, etc. así como en bibliotecas desistemas operativos. Se ha recomendado en varios libros, por ejemplo, Bratley, Fox, ySchrage (1987) y en artículos, Law y Kelton (1982). Como curiosidad, este generadores utilizado en Arena y uno similar se utilizó en AutoMod (con el mismo módulopero con el multiplicador 742938285) hasta hace poco, cuando los gestores de estosproductos tuvieron la buena idea de reemplazarlo por MRG32k3a. También se utilizaen otros productos de sofware de simulación.

MRG32k3a Este es el generador propuesto por L’Ecuyer (1999a). Combina dosgeneradores multiplicativos de orden 3 y su periodo tiene una longitud de aproxi-madamente 2191.

MT19937 Es el generador “trenzado” Mersene propuesto por Matsumoto y Nishimu-ra4(1998). Tiene el periodo de longitud más grande: 219937 − 1.

Generadores de Tausworthe.

Están relacionados con los métodos criptográficos, operan sobre los bits para for-mar números aleatorios. Se define una secuencia b1, b2, . . . de dígitos binarios por re-currencia.

3(http://support.microsoft.com/directory)4Matsumoto, 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

7

bi = (c1bi−1 + c2bi−2 + . . . + cqbi−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 deeste 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ónXOR.

bi ={

0 si bi−r = bi−q

1 si bi−r 6= bi−q

Para inicializar la secuencia hay que dar valor a los primeros bi hasta q. Esto es equiv-alente 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 2Los primeros 42 bits son :︷︸︸︷

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 dichasecuencia 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 ejemploanterior si se elige l = 4, se obtiene la secuencia:

1516

,816

,1316

,1316

,416

,216

,516

,916

,1516

,116

, . . .

l debe ser como máximo la longitud de palabra del ordenador, también se puedenelegir 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 yaún mayores incluso en micros de 16 bits.

8

Inconvenientes:

- Aunque globalmente suelen presentar propiedades de aleatoriedad buenas, aveces, localmente no las tienen.

- En general proporcionan malos resultados en las pruebas de rachas hacia arribay hacia abajo.

- Aunque la correlación de primer orden (un número con el siguiente) es casi cerose sospecha que algunos generadores pueden proporcionar valores elevados decorrelaciones de orden elevado.

- No todos los polinomios primitivos tiene las mismas cualidades.

Combinaciones de generadores.

Si se combinan dos secuencias de números aleatorios Xi e Yi se produce una nuevasecuencia Zi que reduce la no aleatoriedad. Si Xi e Yi están distribuidas sobre el rango0 → (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. Com-binando un generador de congruencias con otro de Tausworthe se pueden atenuarmuchas de las situaciones no deseadas que se producen. Wichmann y Hill5 informande buenos resultados utilizando el método 1), es decir, combinando generadores (eneste caso tres):

Xi+1 = 171 Xi mod 30269Yi+1 = 172 Yi mod 30307Zi+1 = 170 Zi mod 30323

Por tanto Ui+1 =Xi+1

30269+

Yi+1

30307+

Zi+1

30323

El periodo es del orden de 1013. Estos cálculos no precisan mas de 16 bits por loque se considera un buen generador para ordenadores de 16 bits 6.

5Wichman, B.A. and I.D. Hill, Algorithm AS 183: An Efficient and Portable Pseudo-Random Number Gen-erator, Applied Statistics, 31, 188-190, 1982.

6En una nota publicada por Microsoft, indican que la versión del RANDU de 2003 implementa estealgoritmo. http://support.microsoft.com/kb/828795

9

L’Ecuyer (1988)7 recomienda combinar:

Xi+1 = 40014×Xi mod 2147483563Yi+1 = 40692× Yi mod 2147483399

en la forma:

Xi+1 = (Xi − Yi) mod 2147483562

El periodo obtenido es (m1 − 1) ∗ (m2 − 1)/2 ≈ 2 1018, con m1 = 2147483563 ym2 = 2147483399. L’Ecuyer y otros autores han planteado, como se ha indicado ante-riormente, nuevas combinaciones de generadores para obtener periodos mucho máslargos.

Consideraciones sobre la selección de la semilla.

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

1. No utilizar el ceroAunque puede ser una semilla válida para los generadores de congruencias lin-eales mixtos, llevará a cero a los generadores de congruencias lineales multi-plicativos o a los de Tausworthe.

2. Prohibir valores paresEn teoría los valores pares son igual de buenos que los impares. De hecho paralos generadores de periodo total, todos los valores de semilla no cero son igual-mente buenos. Si un generador no es de periodo total, por ejemplo uno multi-plicativo con módulo m = 2k, la semilla debe ser impar. Para impedir este tipode situaciones por desconocimiento del tipo de generador es recomendable uti-lizar semillas impares.

3. No subdividir una serie (stream)Utilizar una única serie para todas las variables es un error muy común. porejemplo, si u1, u2, . . . es la secuencia generada por u0, se podría pensar en utilizar

7P. L’Ecuyer. Efficient and portable combined random number generators. Communications of the ACM,(31), 6, 742–751, 1988

10

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 entrelas variables.

4. Utilizar series (stream) no solapadasCada serie precisa una semilla separada. Si las semillas son tales que las seriesse solapan, existirá una correlación entre las series y las secuencias obtenidas noserá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 cuales el número siguiente :

xn = anx0 + c(an−1)a−1 mod m

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

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

6. No utilizar semillas aleatoriasNo se recomienda seleccionar semillas aleatorias. En particular no utilizar dosnú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.

11

Referencias

Las referencias básicas para este punto son:

Jerry Banks, John S. Carson, Barry L. Nelson Discrete-Event System SimulationPrentice-Hall, 1996.

James E. Gentle Random Number Generation and Monte Carlo Methods Springer-Verlag, 1998.

Raj Jain The Art of Computer Systems Performance Analysis Techniques for Experi-mental 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 muycompleto respecto a la selección de los parámetros de los generadores. Springer-Verlag, 2001.

Pierre L’Ecuyer Software for Uniform Random Number Generation: Distinguishingthe good and the bad Winter Simulation Conference, 2001

En la dirección http://csrc.nist.gov/rng/rng6_4.html hay documentación y distintostipos de pruebas para números aleatorios.

12