Generadores de números aleatorios

Georgina Flesia

FaMAF

26 de marzo, 2015

¿Que son los numeros aleatorios?

DefiniciónLos numeros aleatorios son observaciones independientes devariables aleatorias distribuidas uniformemente en el (0,1).Observaciones de variables aleatorias no uniformes se obtienengeneralmente a travez de transformaciones de muestras de numerosaleatorios.

¿Para qué se utilizan?

I Simulación.I Muestreo.I Análisis numérico.I Testeo de programas.I Juegos de azar.I Toma de decisiones.

Secuencias de números aleatorios

En simulación las secuencias con distribución uniforme en [0,1] seutilizan:

I en forma directa,I para generar distribuciones discretas y continuas,I para generar valores de variables aleatorias dependientes.

Un poco de historia

Antes de las computadoras, existieron diferentes métodos paragenerar secuencias de números aleatorios:

I Procedimientos físicos (monedas, dados, bolilleros, . . . ).I (1927) Tipett: tabla de 40000 dígitos aleatorios (no resultaron

con distribución uniforme).I (1939) Kendall y Babbington: dispositivo mecánico. Tabla de

100.000 números aleatorios.I (1955) Rand Corporation: ruido electrónico. Tabla de 1 millón de

números aleatorios.

Inconvenientes de métodos físicos

1. No puede repetirse una misma secuencia.2. No hay velocidad computacional.3. Incorporar una tabla a la computadora implica gran costo de

almacenamiento en relación a la cantidad de números.

Con la mejora en la rapidez de cálculo de las computadoras surgenmétodos de generación de secuencias de númerospseudo-aleatorios.

DefiniciónUna secuencia de números pseudo-aleatorios es una sucesión denumeros producidos en forma determinista que tienen la aparienciade ser variables aleatorias uniformes independientes en el (0,1)

Propiedades deseables de un generador de númerospseudo-aleatorios

Un generador de números aleatorios razonable debe cumplir:1. repetibilidad

Al repetir los parámetros del generador, se repite la secuencia.2. portabilidad,

La secuencia no debe depender del lenguaje computacional nide la computadora utilizada.

3. velocidad computacional.Puede ser conveniente utilizar lenguajes de bajo o medio nivel.

Principios generales de un buen generador

I La secuencia generada debe ser intuitivamente aleatoria.I Esa aleatoriedad debe ser establecida teóricamente o, al menos,

debe pasar ciertos tests de aleatoriedad.I Debe conocerse algo sobre las propiedades teóricas del

generador.

Un ejemplo

Secuencia de von Neumann (1946):1. X0: número de 5 dígitos. (03001)2. X 2

i : escrito con diez dígitos. (0009006001)3. Xi+1: 5 dígitos centrales. (09006)4. Volver a 2

3001,9006,81108,78507,63349,13095,71479,92474,51440

21000,41000,81000,61000,21000 . . .

con 4 dígitos: 3792,3792,3792,3792, . . .

Desventajas del Gen. de von Neuman

1. No se conocen propiedades teóricas del generador: (semillaconveniente, número de dígitos en cada término).

2. Forsythe: (con 4 dígitos) De 16 semillas, 12 degeneraban en6100, 2100, 4100, 6100,.... y 4 degeneraban en 0.

3. Metropolis: Secuencias de 20 bits degeneran en uno de 13ciclos, con longitud máxima 142.

4. En algunos casos el primer tramo es satisfactoriamente aleatorioy luego degenera.

Generador congruencial lineal

Genera una sucesion de enteros con la relacion recursiva

yi = ayi−1 + c mod M, 1 ≤ i ,

y los transforma a numeros del (0,1)

xi =yi

M, secuencia en el [0,1).

I y0: semillaI a: multiplicadorI c: incrementoI M: módulo

generador mixto: c 6= 0 generador multiplicativo: c = 0.

Ejemplo

¿Cuál es el siguiente número en la secuencia, entre 0 y 15?

0,1,6,15,12,13,2,11,8,9,14, ..

yi = 5yi−1 + 1 mod 16, y0 = 0.

Secuencia intuitivamente no aleatoria:

1,12,1,12,1,12,1,12 . . .

Propiedades

I El menor número K tal que yn+K = yn es el período de lasecuencia.

I Todo generador congruencial genera secuencias de períodofinito.

I El período de una secuencia está acotado por M.I Repetibilidad.I ¿Portabilidad?

Elección de a, c y M

Las buenas propiedades dependen de una elección apropiada de a,c y M, y en algunos casos y0.

I La elección de M se relaciona con: longitud de la secuencia yvelocidad computacional.

I La elección de a y c, en función de M, se relacionan con laaleatoriedad.

Generadores mixtos

yi+1 = a yi + c mod M, c 6= 0

tiene período M si y sólo si

1. m.c.d .(c,M) = 12. a ≡ 1 mod p, para cualquier factor primo p de M.3. Si 4 | M, entonces a ≡ 1 mod 4.

Corolario: Si M es primo, el período máximo ocurre sólo si a = 1.

Generadores multiplicativos

I a es raíz primitiva de M si a(M−1)/p 6≡ 1 mod (M) para cualquierfactor primo p de M − 1.

Ejemplo: M = 7.

a = 2 a = 321 ≡ 2 mod 7 31 ≡ 3 mod 722 ≡ 4 mod 7 32 ≡ 2 mod 723 ≡ 1 mod 7 33 ≡ 6 mod 724 ≡ 2 mod 7 34 ≡ 4 mod 725 ≡ 4 mod 7 35 ≡ 5 mod 726 ≡ 1 mod 7 36 ≡ 1 mod 727 ≡ 2 mod 7 37 ≡ 3 mod 7

Generadores multiplicativos

Para un generador multiplicativo yi+1 = a yi mod M, la longitud K dela secuencia verifica

1. Si K = M − 1 entonces M es primo.2. K divide a M − 1.3. K = M − 1 si y sólo si a es raíz primitiva de M.

Problema: Encontrar raíces primitivas.

Propiedad útil: Si a es raíz primitiva y (k ,M − 1) = 1, entonces ak esraíz primitiva.

Un ejemplo con M primo

M = 231 − 1 a = 16807

M − 1 = 231 − 2 = 2 · 32 · 7 · 11 · 31 · 151 · 331,

I M: es un primo de Mersenne.I 7: raíz primitiva, (5,M − 1) = 1, implica que 75 = 16807 es raíz

primitiva.I secuencia de longitud máxima M − 1 = 2 147 483 646.

Módulo potencia de 2

Si M = 2k , c = 0, tomar módulo es computacionalmente sencillo.

yj = ajy0 mod (2k )

I secuencia de longitud máxima = 2k−2, para a raíz primitiva.I facilita cálculos (desplazamiento de bits).

I fenómeno de no aleatoriedad en bits menos significativos.I RANDU: M = 231, a = 216 + 3 = 65539.

Desventaja de un generador congruencial

En una secuencia y1, y2, . . . dada por un generador congruencialcualquiera, los puntos

(yj , yj+1, . . . , yj+k−1), j = 0,1,2, . . .

están ubicados en no más de (k !M)1/k hiperplanos paralelos.

I Cota máxima: (k !M)1/k : Estructura de redI Generador RANDU: Ternas ubicadas en 15 planos paralelos.

Generadores congruenciales, m=256

Generadores congruenciales

RANDU

Generadores portables

a = 75 = 16807

M = 231 − 1 = 2147483647

I "Buen generador".I Las multiplicaciones superan el rango de 32 bits.I Algoritmo de Schrage.

Algoritmo de Schrage

M = aq + r

q = [M/a], r = M mod a

Si r < q y 0 < z < M − 1 se puede probar que para todo z,0 < z < M:

I 0 ≤ a(z mod q) ≤ M − 1I 0 ≤ r [z/q] ≤ M − 1I

az mod M =

{a(z mod q)− r [z/q] si es ≥ 0a(z mod q)− r [z/q] + M c.c.

El generador ran0

yj+1 = ayj mod M

a = 75 = 16807, M = 231 − 1 = 2147483647

Schrage: se utiliza q = 127773 y r = 2836Desventajas:

I Sucesiones de números muy pequeños.I Inconvenientes en el plano: (yi , yi+1): el test χ2 falla para

N ≈ O(107) ≤ M − 2

Shuffling

Se ”almacenan” los últimos valores generados en una tabla, y lasalida se obtiene eligiendo aleatoriamente un elemento de dichatabla y reponiéndolo por un nuevo valor generado.

I El generador ran1.I a = 75 = 16807, M = 231 − 1 = 2147483647I Tabla de 32 posiciones.

El generador ran1

RAN

SALIDA

1

y

v0

v31

2

3

Combinación de generadores

TeoremaSean W1,W2, . . . ,Wn variables aleatorias discretas, tales queW1 ∼ U([0,d − 1]). Entonces

W = (n∑

j=1

Wj) mod d

es una v.a. uniforme discreta en [0,d − 1].Ejemplo: tirar 2 dados, y sumar módulo 6.

Combinación de congruenciales

I Combinar secuencias de generadores congruenciales.Sugerencia: restar.

I Se obtiene un generador de v.a. uniformes.I La longitud de la secuencia es mayor = mínimo común múltiplo

de los generadores.

xn = 40014xn−1 mod 231 − 85

yn = 40692yn−1 mod 231 − 249

El 2 es el único factor común.Período ≈ 2.3× 1018

Combinación de congruenciales

El generador ran2

I Utiliza 2 generadores congruenciales de enteros de 32 bits.I Tabla de 32 posiciones: shuffling.I Salida = combinación de x e y .

G1 7→ xn = 40014xn−1 mod 231 − 85

G2 7→ yn = 40692yn−1 mod 231 − 249

ran2

y

G1

G2

2

4

v31

v0

SALIDA=vj−y

1

3

vj

El generador ran3

Ver en Knuth, Seminumerical Algorithms, sección 3.2.1.

I Se construye inicialmente una tabla de 55 posiciones, connúmeros aleatorios.

I i y j recorren cíclicamente la tabla, separados por 24 posiciones.I Salida: diferencia entre los elementos i y j , que reemplaza a su

vez el lugar i .I El ciclo que se obtiene es de longitud 232.

ran3

SALIDA=vi−vj

j vj

vii

1

2

3

j−i= 24 (55)

j−>j−1i−>i−1

Marsaglia - Zaman 1994

Some portable very-long-period random number generators, GeorgeMarsaglia and Arif Zaman. Proponen mejoras sobre la definicion deran2

I La combinación de generadores produce mejoras. Sin embargo,conviene combinar estructuras diferentes.

II Utilizar módulo 32 bits en lugar de 31 bits.I Utilizar ”primos seguros” (safeprimes), con a = 2k + α.

231 − 69 231 − 535

232 − 1409, 232 − 209

Otras consideraciones sobre ran2

I ¿Por qué no un rango de enteros?I La subrutina iran2( ).

UNI( )=.5+.232830E-9*iran2()

VNI( )=.4656613E-9*iran2()

Algunos ejemplos de candidatos para combinar

Módulo Secuencia Período232 xn = 69069 xn−1 + impar 232

232 xn = xn−1 ∗ xn−2 231

232 xn = xn−1 + xn−2 + C 258

231 − 69 xn = xn−3 − xn−1 262

232 − 18 xn = xn−2 − xn−3 − C 295

Ejemplos: generador mzran( )

Combina los generadores:

xn = 69069xn−1 + 1013904243 mod 232

Período: ≈ 232

xn = xn−3 − xn−1 mod 231 − 69

Período: ≈ 294

El período es mayor a 294, o 1028.

Ejemplos: generador mzran13( )

Combina los generadores:

xn = 69069xn−1 + 1013904243 mod 232

Período: ≈ 232

xn = xn−2 − xn−3 − c′ mod 232 − 18

Período: ≈ 295

El período es del orden de 2125, o 1036.

Rutina propuesta por Marsaglia y Zaman

typedef unsigned long int unlong;unlong x=521288629, y=362436069, z=16163801,

c=1, n=1131199209;unlong mzran13(){ long int s;

if (y>x+c) (s=y-(x+c); c=0;}else { s=y-(x+c)-18; c=1; }x=y; y=z;return (z=s) + (n=69069*n+1013904243);

};void ran13set(unlong xx, unlong yy, unlong zz, long nn)

{ x=xx; y=yy; z=zz; n=nn; c=y>z; }

Generadores de Fibonacci demorados (LaggedFibonacci)

Se basan en la ecuacion

xi = (xj−1 � xi−k ) mod m 0 < j < k

donde � representa un operador binario como suma, multiplicación o(XOR).

I Típicamente, m = 2M con M = 32 o 64.I Necesita un bloque semilla de tamaño k para ser inicializado,

que usualmente se genera usando la rutina rand2().

Generador GFSR

I Cuando el operador es el (XOR), el generador de Fibonacci sedenomina two-tap generalized feedback tap register (GFSR).

I La teoría es bastante compleja, y no hay pruebas rigurosas deldesempeño de estos generadores.

I Su calidad se basa en test estadísticos.I Mas aun, los valores j y k tienen que ser elegidos

cuidadosamente.I Para que el generador alcance su máximo periodo el polinomio

y = xk + x j + 1

debe ser primitivo sobre los enteros modulo 2.I El algoritmo r1279(), multiplicative lagged Fibonacci, con

k = 1279 pasa todos los test conocidos y tieneimplementaciones rápidas en librerías como GSL (C++ y C).

Mersenne Twister

I Mersenne Twister fue desarrollado en 1997 por Matsumoto yNishimura y es una version de un two-tap generalised feedbacktap register (GFSR).

I La implementación mt19937() es parte de muchas lenguajes ylibrerías científicas como Matlab, R, Python, Boost, GSL.

I Tiene un periodo de ρ = 219937 − 1. Para palabras de k -bits,mt19937() genera numeros aleatorios con distribuciónpseudo-uniforme en el rango [0,2k − 1].

I El algoritmo es complejo, no es recomendable hacer unaimplementación sucia, sino utilizar una implementaciónreconocida, con desempeño publicado.

good C random number generator

13/05/2003Hi everybody, Does anyone know of a good random numbergenerator in C? I know of the random generators in NumericalRecipes. Are there other good sources? Thank you in advance.David

George Marsaglia (March 12, 1924 – February 15, 2011)

Most RNGs work by keeping a certain number, say k, of the mostrecently generated integers, then return the next integer as a functionof those k. The initial k integers, the seeds, are assumed to berandomly chosen, usually 32-bits. The period of the RNG is related tothe number of choices for seeds, usually 2(32k), so to get longerperiods you need to increase k.

Rutina propuesta

Probably the most common type has k=1, and needs a single seed,with each new integer a function of the previous one. An example isthis congruential RNG, a form of which was the system RNG in VAXsfor many years:

static unsigned long x=123456789;/* a random initial x to be */

/* assigned by the calling program */unsigned long cong(void ){ return (x=69069*x+362437);}

Simple, k=1, RNGs can perform fairly well in tests of randomnesssuch as those in the new version of Diehard,

csis.hku.hk/~diehard

but experience has shown that better performances come from RNGswith k’s ranging from 4 or 5 to as much as 4097.

Rutina propuesta

Here is an example with k=5, period about 2160, one of the fastestlong period RNGs, returns more than 120 million random 32-bitintegers/second (1.8MHz CPU), seems to pass all tests:

static unsigned longx=123456789,y=362436069,z=521288629,w=88675123,v=886756453;/* replace defaults with five random seed valuesin calling program */

unsigned long xorshift(void){unsigned long t;t=(x^(x>>7)); x=y; y=z; z=w; w=v;v=(v^(v<<6))^(t^(t<<13)); return (y+y+1)*v;}

Rutina propuesta

Another example has k=257, period about 28222. Uses a static arrayQ[256] and an initial carry ’c’, the Q array filled with 256 random32-bit integers in the calling program and an initial carryc<809430660 for the multiply-with-carry operation. It is very fast andseems to pass all tests.

static unsigned long Q[256],c=362436;/* choose random initial c<809430660

and *//* 256 random 32-bit integers for Q[]*/unsigned long MWC256(void){

unsigned long long t,a=809430660LL;static unsigned char i=255;

t=a*Q[++i]+c; c=(t>>32);return(Q[i]=t); }

The Mersenne Twister (check Google) is an excellent RNG, withk=624. But it requires an elaborate C program and is slower thanmany RNGs that do as well in tests, have comparable or longerperiods and require only a few lines of code.

Rutina propuesta

Here is a complimentary-multiply-with-carry RNG with k=4097 and anear-record period, more than 1033000 times as long as that of theTwister. (2131104vs.219937)

static unsigned long Q[4096],c=362436;/* choose random initial c<809430660 and *//* 4096 random 32-bit integers for Q[] */

unsigned long CMWC4096(void){unsigned long long t, a=18782LL;static unsigned long i=4095;unsigned long x,r=0xfffffffe;

i=(i+1)&4095;t=a*Q[i]+c;c=(t>>32); x=t+c; if(x<c){x++;c++;}return(Q[i]=r-x); }

You will find several more CMWC RNGs and comments on choice ofseeds in the May 2003 Communications of the ACM.