Métodos congruenciales

Técnicas para generar números aleatorios

MÉTODOS CONGRUENCIALES Los principales generadores de números pseudo-aleatorios utilizados hoy en día son los llamados generadores congruenciales lineales, introducidos por Lehmer en 1951. Un método congruencial comienza con un valor inicial (semilla) x0, y los sucesivos valores xn,n > 1 se obtienen recursivamente con la siguiente fórmula:

xn = axn-1 + b modulo m, donde a, m y b son enteros positivos que se denominan, respectivamente, el multiplicador, el modulo y el incremento. Si b = 0, el generador se denomina multiplicativo; en caso contrario se llama mixto. La sucesión de números pseudo-aleatorios un,n > 1 se obtiene haciendo u, = Como el siguiente resultado demuestra, cada xi esta completamente caracterizado por a,b,m y x0.

Proposición 1.3. Los valores generados por un método congruencial verifican: an —

xn = anx0 + b mod m

1

a–1

Demostración. Para n =1, tenemos que xi = ax0 + b mod m, lo que implica que existe k E Z+ tal que km + x1 = ax0 + b. Ahora, para n = 2, se tiene que x2

=

ax1 + b mod m

=

a(ax0 + b — km) + b mod

=

a2x0 + b(a + 1) — akm mod m

=

a2x0 + b(a + 1) mod m.

Por recurrencia, se tiene que

x = a x + b(a + a + 1) mod m 3

3

2

0

x = a x + b(a - 1 + …. + a + 1) mod m n

n

n

0

= ax + b n

0

m


Por lo tanto, la primera objeción que se le puede hacer a este método, objeción común a todos los generadores de numeros pseudo-aleatorios, es que la sucesión de los valores xn no es en absoluto aleatoria.

Sin embargo, posteriormente

veremos que

si

elegimos las parámetros

iniciales

convenientemente, la sucesión {un} puede asemejarse a una sucesión de números aleatorios.

La segunda objeciín es que los valores u pueden tomar sílo los valores 0, mm, . . . , m-1, luego no hay posibilidad de generar un valor, por ejemplo entre m y m. Tomando m suficientemente grande m > 10 9, el conjunto de posibles valores es suficientemente denso en [0,1] como para que la sucesión asemeje a la de una variable continua uniforme en dicho intervalo. Ejemplo 1.4. Consideremos el generador congruencial xn = 5xn-1 + 1 mod 9 con x0 = 1.

Ejemplo Consideremos el generador congruencial xn = 5xn-1 + 3 mod 16, con xo = 7.

Se observa que el generador del ejemplo 1.4 no es adecuado, pues produce un ciclo de longitud 6 <


m = 9 (la longitud de ciclo de un operador se denomina periodo). Sin embargo, en el generador del ejemplo 1.5 la longitud de ciclo coincide exactamente con el módulo, lo cual es inevitable, pues es evidente que el periodo nunca puede exceder al módulo. En este caso, se dice que el generador es de ciclo completo. Obsérvese que un generador sea de ciclo completo es independiente de la semilla que se utilice; sea cual sea x0 E { 1 , . . . , m}, el ciclo siempre tendrá periodo m y se producirá en el mismo orden. Sin embargo, si un generador no es de ciclo completo, la longitud de ciclo puede depender de la semilla utilizada. Si en el ejemplo 1.4 utilizamos x0 = 5 o x0= 8, el periodo resultante es 2. De los ejemplos anteriores se desprende que una cuestión de interés es como elegir los parámetros del generador de forma que este tenga ciclo completo. El siguiente teorema, propuesto por Hull y Dobell (1962) proporciona un caracterización en este sentido.

Teorema 1.6.Un generador congruencial tiene periodo completo si y solo si se cumplen las siguientes condiciones: 1. m y b son primos entre si. 2. Si q es un numero primo que divide a m, entonces q divide a a - 1. 3. Si 4 divide a m, entonces 4 divide a a - 1. Controlario 1.7 Un generador congruencial multiplicativo no puede tener periodo completo Es evidente que es necesario que m sea grande con el fin de un periodo largo y una alta densidad en el intervalo [0,1]. Sin embargo, la operación de dividir por m y calcular el resto es relativamente lenta. Una elección de m adecuada computacional- mente es m = 2k, donde k-bits es el tamaño de palabra (unidad básica de trabajo) del microprocesador. El hecho de que esta elección sea ventajosa reside en que nos podemos aprovechar del desbordamiento de datos (integer overflow) para no tener que realizar la operación del generador explícitamente. La cuestión es que si tenemos una máquina de k-bits, entonces el mayor entero que puede ser representado es 2k — 1 y en cualquier intento de representar un entero mayor H, que ocuparía h > k dígitos binarios, se perderían los h — k dígitos binarios más a la izquierda, y los k dígitos que quedan se corresponden precisamente con H mod 2k. Ejemplo 1.8 Para ilustrar como aprovechar el desbordamiento de datos en el generador del ejemplo 1.5,


veamos cómo generar x10 = 9 a partir de x9 = 14. Observemos que 5x14 + 3 = 73, que en binario es 1001001. Como nuestra capacidad es de 4 bits, el desbordamiento de datos produce que se pierdan los tres primeros dígitos, quedando 1001 que es la representación binaria de x10 = 9.

Corolario 1.9. Un generador congruencial con m = 2k > 4 tiene periodo completo si y solo si b es impar y 1 = a mod 4.

Los generadores multiplicativos tienen la ventaja computacional de que no es necesario la suma de la constante b, pero, como hemos visto, no pueden tener periodo completo. Afortunadamente, es posible obtener periodo m — 1 si m y a se eligen adecuadamente. Los generadores multiplicativos se introdujeron antes que los mixtos, han sido objeto de más estudio y actualmente se utilizan mas que los mixtos. Al igual que los generadores mixtos, la mejor forma de elegir el modulo desde el punto de vista computacional es tomar m = 2k. Sin embargo, en este caso, el periodo máximo posible será una cuarta parte del módulo.

Observemos que 5x14 + 3 = 73, que en binario es 1001001. Como nuestra capacidad es de 4 bits, el desbordamiento de datos produce que se pierdan los tres primeros dígitos, quedando 1001 que es la representación binaria de x10 = 9.

Corolario 1.9. Un generador congruencial con m = 2k > 4 tiene periodo completo si y solo si b es impar y 1 = a mod 4.

Los generadores multiplicativos tienen la ventaja computacional de que no es necesario la suma de la constante b, pero, como hemos visto, no pueden tener periodo completo. Afortunadamente, es posible obtener periodo m — 1 si m y a se eligen adecuadamente. Los generadores multiplicativos se introdujeron antes que los mixtos, han sido objeto de más estudio y actualmente se utilizan más que los mixtos.


Al igual que los generadores mixtos, la mejor forma de elegir el modulo desde el punto de vista computacional es tomar m = 2k. Sin embargo, en este caso, el periodo máximo posible será una cuarta parte del módulo. Proposición 1.10. El periodo máximo de un generador congruencial multiplicativo con m = 2k > 16 es 2k_2. Dicho periodo maximal se alcance si y solo si x0 es impar y 3 = a mod 8 o 5 = a mod 8

Un generador multiplicativo muy utilizado, conocido como RANDU, tomaba m = 231 y a = 216 + 3. Sin embargo, se ha demostrado que tiene propiedades estadísticas bastante malas. Debido a los problemas que surgen de elegir m = 2k en los generadores multiplicativos, se han investigado otras alternativas para el valor de m. La más satisfactoria es la propuesta de Hutchinson en 1966, que consiste en tomar como m el mayor primo menor que 2 k. Por ejemplo, si k = 31, el mayor primo menor que 231 es 231 — 1. A este respecto se tiene el siguiente resultado.

Teorema 1.11. Sea t la longitud de un ciclo maximal de un generador congruencial multiplicativo. Se verifica que:

1.

Si t = m — 1, entonces m es primo

2.

Si m es primo, entonces t divide a m — 1.

3.

Si m es primo, entonces t = m — 1 si y sólo si a es una raóz primitiva de m (a es una raíz

Primitiva tal que 1 =

Los generadores multiplicativos mas famosos utilizados por IBM tomaban m = 2 31 — 1 y a1 = 75 o a2 = 630360016.


MÉTODO FIBONACCI EXTENDIDOS Historia La sucesión fue descrita por Fibonacci como la solución a un problema de la cría de conejos: " Cierto

hombre tenía una pareja de conejos juntos en un lugar cerrado y uno desea saber cuántos son creados a partir de este par en un año cuando es su naturaleza parir otro par en un simple mes, y en el segundo mes los nacidos parir también".[1]

Número de Mes

Explicación de la genealogía

Comienzo del mes 1

Nace una pareja de conejos (pareja A).

Fin del mes 1

La pareja A tiene un mes de edad. Se cruza la pareja A.

Fin del mes 2

La pareja A da a luz a la pareja B. Se vuelve a cruzar la pareja A.

Fin del mes 3

La pareja A da a luz a la pareja C. La pareja B cumple 1 mes. Se cruzan las parejas A y B

Fin del mes 4

Las parejas A y B dan a luz a D y E. La pareja C cumple 1 mes. Se cruzan las parejas A,

Fin del mes 5

A, B y C dan a luz a F, G y H. D y E cumplen un mes. Se cruzan A, B, C, D y E.

Fin del mes 6

A, B, C, D y E dan a luz a I, J, K, L y M. F, G y H cumplen un mes. Se cruzan A, B, C, D,

Nota: al contar la cantidad de letras distintas en cada mes, se puede saber la cantidad de parejas totales que hay hasta ese mes. De esta manera Fibonacci presentó la sucesión en su libro Liber Abaci, publicado en 1202. Muchas propiedades de la sucesión de Fibonacci fueron descubiertas por Édouard Lucas, responsable de


haberla denominado como se la conoce en la actualidad.[2] También Kepler describió los números de Fibonacci, y el matemático escocés Robert Simson descubrió en 1753 que la relación entre dos números de Fibonacci sucesivos

se acerca a la

relaciónáurea fi ( ) cuanto más se acerque a infinito; es más: el cociente de dos términos sucesivos de toda sucesión recurrente de orden dos tiende al mismo límite. Esta sucesión ha tenido popularidad en el siglo XX especialmente en el ámbito musical, en el que compositores con tanto renombre como Béla Bartók, Olivier Messiaen, la banda Tool y Delia Derbyshire la han utilizado para la creación de acordes y de nuevas estructuras de frases musicale

Las grandes ventajas de estos generadores son su rapidez y que tienen un periodo muy largo. La fomentación teórica en la que se basan es diferente a la de los generadores congruenciales lineales. Los generadores de Fibonacci se basan en una recurrencia del tipo

Donde r s son enteros dados y

denota alguna de las operaciones

TRANSFORMACIÓN INVERSA La técnica de la transformación inversa es probablemente la más fácil de entender y la más comúnmente utilizada, este método consiste en utilizar la función de distribución acumulada F , y se basa en el siguiente Teorema Teorema

Sea

U

una variable aleatoria uniforme en (0,1) y

función de distribución

F , si F

X

una variable aleatoria continua con

es invertible entonces la variable aleatoria X =

una distribución F .

Demostración Sea F X denota la función de distribución de X = F _1(U). Entonces

F _1(U) tiene


Ahora como F es una función de distribución invertible, se tiene que creciente de x, así la desigualdad a <

b

F (x) es una función monótona no-

es equivalente a la desigualdad F (a) < F (b). Por lo tanto

Sea X una variable aleatoria con soporte (a,b) y función de densidad muestra aleatoria de tamaño

F (x),

n

f

(x), se desea generar una

de esta distribución. Denótese a la función de distribución acumulada

entonces utilizando el resultado del teorema anterior, y considerando el hecho de que F (x) es

monótona creciente en el soporte de la variable aleatoria en donde no es constante, entonces en este intervalo

la función de distribución tiene inversa, sea ésta F- 1, entonces se propone el

siguiente algoritmo para esta técnica:

a)

Generar n valores u i ~ U (0,1)

b)

Regresar x ¡ = F- l (u i )

Es decir, Tomando una muestra distribuida uniformemente a lo largo del eje y entre 0 y 1, y posteriormente evaluando la distribución inversa en este valor, se está generando una muestra de la variable aleatoria X. Puede notarse que en donde la función de distribución acumulativa F(x) tiene una pendiente relativamente mayor, dará como resultado una alta densidad de puntos en el eje x (dando un mayor valor de f (x)), y por otro lado, donde F(x), tiene una pendiente relativamente menor (crece más despacio), dará como resultado una densidad menor de puntos a lo largo del eje x (con un menor valor de función de densidad f (x)).


EJEMPLO 1

Esta técnica puede ser ilustrada con la distribución Weibull. En este caso, se tiene que F(x) = 1 , luego si U ~U(0,1) y U F(x) entonces se encuentra que x = b[-ln(1-u)]1/c . Realizando el programa en el Proyecto R, resulta después de Generar 50,000 números aleatorios.


En las gráficas de la figura 1.1ª) muestran en el primer cuadro los números aleatorios sin agrupación, en el segundo se presenta el histograma de dichos puntos en esta se observa que siguen un patrón parecida a la densidad Weibull, en la segunda imagen figura1.1b) se muestran las dos densidades, de los puntos generados con la técnica y la densidad teórica (línea punteada) en donde se observa que son muy similares de esta forma se puede asegurar que la técnica funciona, además entre mayor sea el número de valores aleatorios mejor es la aproximación. La transformación inversa es una técnica simple, eficiente para obtener la función de densidad de probabilidades, pero tiene un inconveniente muy fuerte, requiere que el decisor sea capaz de encontrar la inversa de la función de distribución F(x), como esto no siempre es posible, es necesario considerar otras técnicas. EJEMPLO 2 Si X es una variable aleatoria exponencial con parámetro 1, entonces la función de distribución está dada por.

F(x) 1ex Denotando x = F- 1 (u ), se tiene F(x) 1e

-x

ó F(x)=1e-x

Tomando logaritmos, x = log(1u) Por lo tanto, se puede generar una exponencial con parámetro 1, generando Una variable aleatoria U, estableciendo.

x F 1(U) log(1u) 


Se puede obtener por notación que 1U es también uniforme sobre (0,1) y por lo tanto log(1u) tiene la misma distribución como - log(U). Esto es, el logaritmo negativo de un número aleatorio distribuido exponencialmente con parámetro 1. Además, note que si X es exponencial con media 1, entonces para alguna constante positiva c , cX es una exponencial con media c . Por lo tanto una variable aleatoria exponencial X con parámetro λ (media 1/λ) puede ser generada por número aleatorio U y estableciendo X=-

A continuación se muestran las gráficas hechas en el Proyecto estadístico R para generar 50,000 números aleatorios con una distribución exponencial con parámetro 1.


Generar Números Aleatorios con Métodos Congruenciales En los experimentos de simulación es necesario el uso de números aleatorios que representen una distribución de probabilidad uniforme. Para tal propósito tenemos que generar números a través de reglas determinísticas, por lo que es más preciso llamarles números pseudoaleatorios. De hecho, podemos decir que “una sucesión cumple la facultad de aleatoria si satisface las pruebas estadísticas de aleatoriedad” Métodos Congruenciales Los generadores de números aleatorios que más se usan son los generadores congruenciales lineales (LCG) ideados por Lehmer. El objetivo de un LCG es generar un valor aleatorio a partir de otro anterior. En este documento estudiaremos los tres métodos congruenciales lineales más populares: Congruencial Mixto La fórmula (o relación de recurrencia) es sencilla: X a X c m n n ( )mod 1 = × + + Donde: · X0 es la semilla · a el multiplicador · c la constante aditiva y · m el módulo A tener en cuenta: Los valores a, X0 y c tienen que ser mayores que cero y la variable m debe ser mayor que las tres anteriores. Para entrar en acción vamos a darle valores arbitrarios a cada uno de estos parámetros y estudiar que reacción tienen en la relación de recurrencia. Supongamos que a = 5, c = 7, X0 = 7 y m = 8. Entonces los resultados son:


Nótese que después de 8 pasadas el valor inicial de X se repite. Decimos entonces que el periodo del generador es 8… igualito al valor del módulo… Eso no siempre es así. Veamos un caso donde el periodo es menor a m. El valor de los parámetros es a = c = X0 = 4 y m = 6. Ahora lo resultados son:

Para que lo anterior quede mas claro, a continuación se muestra un programa en Matlab que genera números aleatorios implementando el generador congruencial lineal mixto. Cabe hacer la aclaración que los números deben estar entre 0 y 1, por lo que los residuos obtenidos con la función mod se deben dividir entre m-1. % Genera un conjunto de numeros pseudo aleatorios por el metodo congruencial lineal a partir de: % Sintaxis: genelin(a, c, n, s), donde: % a es el factor de multiplicacion igual 1 + 4k, % c es el factor aditivo y es primo relativo a m, % n es el numero de datos a generar y % s es la semilla o valor de inicio, que se sugiere sea impar. % g y k son numeros enteros. function [x, r] = genelin(a, c, n, s) i=1; while n > 2^i, i = i + 1; end % m es el modulo que se obtiene m = 2^g > n, m=2^i; for I=1:n u=a*s+c;

x(I)=mod(u,m); r(I)=x(I)/(m-1); format long, r; s=x(I); end


Métodos congruenciales

Recommend Documents