TEMA 2: M´ etodos eto dos de generaci´ on de variables aleatoon rias discretas y continuas 1.
Intr Introd odu ucci´ ci´ on on
Este tema es una revisi´ revisi´on on de los generadores generadores cl´ asicos de n´ umeros pseudoaleartorios umeros uniformes uniformes y de las t´ ecnicas ecnicas usuales para la generaci´ on de variables aleatorias discretas y continuas a partir de una secuencia de n´ umeros pseudoaleatorios. Puesto que ambos umeros aspectos son b´ asicos y ampliamente conocidos, en este tema, se proporcionar´ asicos a s´olo olo una introducci´on on esquem´ atica atica y algor´ algor´ıtmica de los mismos para facilitar f acilitar su implementaci´ on y la resoluci´on on de las actvidades propuestas. En relaci´ on on con la formulaci´ on on de generadores de n´ umeros pseudoaleatorios uniformes, umeros en la Secci´ on 2, se proporcionan los algoritmos y c´ on odigos para su implementaci´ odigos on. o n. En relaci´on on con las t´ecnicas ecnicas de generaci´on on de variables aleatorias, en las Secciones Se cciones 3 y 4, se describen brevemente brevemente los m´etodos etodos y se formulan formulan los algoritmos, algoritmos, as´ as´ı como se muestran muestran algunos ejemplos implementados en c´ odigo odigo MatLab.
2.
Gen Generad erador ores es de n´ umeros pseudoaleatorios umeros De los generadores que se introducen a continuaci´ on on M´etodo etodo de los cuadrados medios M´etodo eto do de Lehmer Leh mer M´etodo eto do Cong C ongrue ruenci ncial, al,
el m´etodo etodo congruencial es el m´ as efectivo, dado que podemos determinar de antemano el as ciclo de repetici´ on on de d e la secuencia secue ncia generada, gener ada, as´ as´ı como com o presenta pr esenta propiedade propi edadess m´as as pr´ proximas ´ a la aleatoriedad.
2.1.
M´ etodo etodo de los cuadrados cuadrados medios medios
Se formula seguidamente un algoritmo que describe la implementaci´ on on del generador de n´ umeros umeros pseudoaleatorios pseudoaleatorios dado por el m´etodo etodo de los cuadrados medios junto con un ejemplo implementado en c´ odigo odigo MatLab. Algoritmo
Paso 1. Elegir una semilla semilla o valor inicial x0 de 2n cifras. Paso 2. Elevar Elevar al cuadrado la semilla x20 . 1
Paso 3. Seleccionar las 2n cifras centrales de x20 . Paso 4. Como siguiente elemento de la secuencia de n´ umeros pseudoaleatorios, considerar las 2n cifras centrales seleccionadas en el paso anterior. Es decir, considerar ahora como semilla, las 2n cifras obtenidas en el paso anterior. Paso 5. Repetir los Pasos 2, 3, 4 y 5 tantas veces como n´ umeros pseudoaleatorios se deseen generar. Un inconveniente de este m´etodo es que el ciclo maximal de repetici´ on no puede ser determinado previamente y depende de la elecci´ on de la semilla, present´andose con frecuencia secuencias degeneradas. Ejemplo 1.
Seguidamente se muestra un ejemplo, implementado en MatLab, considerando los siguientes par´ ametros Valor inicial X = 4122. N´umero de elementos de la secuencia m = 6. clear all X = 4122; m = 6; for I = 1 : m Y = X ∗ X ; Z = fix(Y /100); V = fix(Z/(104 )); T = V ∗ (104 ); W = Z − T ; X = W ; end
2.2.
M´ etodo de Lehmer
Se formula ahora un algoritmo que describe la implementaci´ on del generador de n´umeros pseudoaleatorios dado por el m´etodo de Lehmer junto con un ejemplo implementado en c´ odigo MatLab.
2
Algoritmo
Paso 1. Elegir una semilla o valor inicial x0 de n cifras y multiplicarla por un n´ umero fijo de k cifras, obteni´endose un n´ umero l de n + k cifras. Paso 2. Se separan las k cifras de la izquierda del n´ umero l del paso anterior, obteni´endose un n´ umero h de n cifras. Paso 3. Se calcula la diferencia d = h − v, siendo v el n´ umero dado por las k cifras de la izquierda del n´ umero l, separadas en el Paso 2. Paso 4. Como nuevo elemento de la secuencia de n´ umeros pseudoaleatorios se considera el valor d, es decir, se considera como nuevo valor de la semilla, el valor de la secuencia generado en el paso anterior. Paso 5. Se repiten los Pasos 2, 3, 4 y 5 tantas veces como n´ umeros pseudoaleatorios se deseen generar. Ejemplo 2.
Seguidamente se muestra un ejemplo, implementado en MatLab, considerando los siguientes par´ ametros Valores iniciales X = 4122 y K = 76. N´umero de elementos de la secuencia m = 16. clear all X = 4122; K = 76; m = 16; for I = 1 : m X 1 = X ; X 2 = K ; Z = X 1 ∗ X 2; W = Z/(104 ); T = fix(W ); S = (W − T ); V = S ∗ (104); F = fix(V ); R = F − T X = R; end 3
2.3.
M´ etodo Congruencial
Se formula ahora un algoritmo que describe la implementaci´ on del generador de n´umeros pseudoaleatorios dado por el m´etodo congruencial junto con un ejemplo en c´ odigo MatLab. Algoritmo
Se considera el m´etodo congruencial mixto con par´ ametros a y b. Paso 1. Elegir una semilla o valor inicial x0. Paso 2. Transformar linealmente x0, obteniendo y = ax 0 + b. Paso 3. Calcular
y y x = m. − m
m
Paso 4. Como nuevo elemento de la secuencia de n´ umeros pseudoaleatorios se considera el valor x del paso anterior. Es decir, como nuevo valor de la semilla se considera el valor de la secuencia generado en el paso anterior. Paso 5. Se repiten los Pasos 2, 3, 4 y 5 tantas veces como n´ umeros pseudoaleatorios se deseen generar. En el algoritmo anterior se ha denotado por [·] la funci´on parte entera. Ejemplo 3.
Seguidamente se muestra un ejemplo, implementado en MatLab, considerando los siguientes par´ ametros para un m´etodo congruencial mixo Semilla s = 20. Escala a = 13. Or´ıgen b = 1. M´odulo de congruencia m = 16.
4
clear all s = 20; a = 13; b = 1; m = 16; for I = 1 : m n = a ∗ s + b; u = mod(n, m); v = u/m s = u; end
3.
Generaci´ on de variables aleatorias discretas
Se describen brevemente los m´etodos usuales de generaci´ on de variables aleatorias discretas a partir de una secuencia de n´ umeros pseudoaleatorios uniformes.
3.1.
M´ etodo de la transformaci´ on cuantil
Para la implementaci´ on de este m´etodo se considera una variable aleatoria discreta X con valores X 1 , . . . Xk , con probabilidades p1 , . . . pk , se tiene entonces el siguiente algoritmo de generaci´on de sus valores a partir de una secuencia de n´ umeros pseudoaleatorios uniformes. Algoritmo
Paso 1. Generar un n´ umero pseudoaleatorio uniforme u. Paso 2. Si el valor u generado en el paso anterior verifica la desigualdad j −1
j
p < u ≤ p , i
i=1
i
j = 1, . . . , k ,
i=1
entonces se considera el valor X j de la variable aleatoria X. Paso 3. Se repiten los Pasos 1, 2 y 3 tantas veces como valores independientes de la variable X se deseen generar. Ejemplo 4.
Seguidamente se muestra un ejemplo, implementado en MatLab, sobre la generaci´ on de una variable aleatoria Binomial a partir de una distribuci´ on de Bernoulli, o equivalentemente, como suma de variables aleatorias independientes de Bernoulli. 5
Semilla s = 20. Escala a = 13. Or´ıgen b = 1. M´odulo de congruencia m = 16. Probabilidad de ´exito de la distribuci´ on Bernoulli p = 1/2. Generaci´on de 10 valores de la Binomial. Par´ ametros de la Binomial n = 16 y p = 1/2. clear all s = 20; a = 13; b = 1; m = 16; p = 0,5; l = 10; for J = 1 : l Bi = 0; for I = 1 : m n = a ∗ s + b; u = mod(n, m); v = u/m; if v < p B = 1; else B = 0; end Bi = Bi + B; s = u; end Binomial(J ) = Bi s = exp(J ); end Binomial = Binomial ;
6
3.2.
M´ etodo de Aceptaci´ on-Rechazo
En los casos en los que la aplicaci´on del m´etodo de la Transformada Inversa requiera elevados tiempos computacionales, o bien no sea viable, se puede considerar, como alternativa, el m´etodo de aceptaci´ on-rechazo, que permite generar la variable de inter´es X, con distribuci´ on de probabilidad { p1 , . . . pk }, a partir de la generaci´ on de otra variable Y, con distribuci´ on de probabilidad { q 1 , . . . qk }, f´acilmente generable (mediante el m´etodo de la transformda inversa, u otro m´etodo). Se supone que existe una constante c tal que p j ≤ c, q j
∀ j = 1, . . . , k .
Algoritmo
Paso 1. Generar un valor de la variable Y, y0 . Paso 2. Generar un n´ umero uniforme u. Paso 3. Si u < py /cq y , entonces generar el valor x0 = y0 de X, en caso contrario no se genera ning´ un valor de X y se va al Paso 1. 0
0
Paso 4 Repetir los Pasos 1, 2 y 3 tantas veces como valores de X se desee generar. Ejemplo 5.
Seguidamente se muestra un ejemplo, implementado en MatLab, sobre la generaci´ on de una variable aleatoria discreta X con diez valores X 1 , . . . , X10 , con probabilidades p1 = 0,11, p2 = 0,12, p3 = 0,09, p4 = 0,08, p5 = 0,12, p6 = 0,10, p7 = 0,09, p8 = 0,09, p9 = 0,10, p10 = 0,10. Semilla s = 20. Escala a = 13. Or´ıgen b = 1. M´odulo de congruencia m = 32. Comparaci´on con una variable aleatoria y que sigue una uniforme discreta con diez valores clear all
7
s = 20; a = 13; b = 1; m = 32; p = 0,5; for I = 1 : m n1 = a ∗ s + b; u1 = mod(n1, m); v1 = u1/m; D1 = fix(10 ∗ v1) + 1; if D1 == 1 P D1 = 0,11; elseif D1 == 2 P D1 = 0,12; elseif D1 == 3 P D1 = 0,09; elseif D1 == 4 P D1 = 0,08; elseif D1 == 5 P D1 = 0,12; elseif D1 == 6 P D1 = 0,10; elseif D1 == 7 P D1 = 0,09; elseif D1 == 8 P D1 = 0,09; elseif D1 == 9 P D1 = 0,10; elseif D1 == 10 P D1 = 0,10; end Q = P D1/0,12; s = u1; n2 = a ∗ s + b; u2 = mod(n1, m); v2 = u2/m; if v2 < Q D2 = D1; s = u2; else I == I + 1; end end
8
3.3.
M´ etodo de Composici´ on
Se aplica a la generaci´ on de variables aleatorias discretas cuya funci´ on masa de probabilidad viene dada por la siguiente expresi´ on: m
m
P [X = j] = α p , i i j
j = 1, . . . , k ,
i=1
α = 1, i
i=1
siendo p ji = P [X i = j ], j = 1, . . . , k , para i = 1, . . . , m , las funciones masa de probabilidad de las variables aleatorias X 1 , . . . , Xm , que intervienen en la definici´on de X. Para este tipo de variables, cuya distribuci´ on es una mixtura o composici´ on de otras, la transformada inversa y t´ecnica de aceptaci´ on-rechazo no suelen ser apropiadas. Se utiliza entonces la t´ecnica de composici´on dada por el siguiente algoritmo: Algoritmo
Paso 1. Generar un valor uniforme u. l−1
l
α < u ≤ α , l = 1, . . . , m , se genera la variable aleatoria X con distribuPaso 2. Si i
i=1
i
l
i=1
ci´on masa de probabilidad p jl , j = 1, . . . , k . Paso 3. Se genera un segundo uniforme para generar la variable aleatoria X l por el m´etodo de la Transformada Inversa (resp. aceptaci´ on-rechazo). Paso 4 Repetir los Pasos 1, 2 y 3 tantas veces como valores de X se desee generar.
4.
Generaci´ on de variables aleatorias continuas
Se introduce brevemente la versi´ on continua de las t´ecnicas de simulaci´ on anteriores para la generaci´on de variables aleatorias con densidad de probabilidad.
4.1.
Transformada Inversa
El m´etodo de la transformada inversa para generaci´ on de variables aleatorias continuas es f´acilmente implementable, seg´ un se aprecia en el siguiente algoritmo, para una variable aleatoria continua X con funci´on de distribuci´on F X . Algoritmo
Paso 1. Generar un n´ umero pseudoaleatorio uniforme u. Paso 2. Calcular X = F X 1 (u). −
9
Paso 3. Repitir los Pasos 1 y 2 tantas veces como valores independientes de la variable X se desee generar. Ejemplo 6.
Seguidamente se muestra un ejemplo, implementado en MatLab, sobre la generaci´ on de una variable aleatoria con distribuci´on exponencial con par´ ametro λ = 1/5. clear all s = 20; a = 13; b = 1; m = 16; lambda = 1/5; for I = 1 : m n = a ∗ s + b; u = mod(n, m); v = u/m; exp = (1/lambda) ∗ log(v); s = u; end
4.2.
T´ ecnica de Aceptaci´ on-Rechazo
Se pretende generar una variable aleatoria X con densidad de probabilidad f X a partir de la generaci´ on de una variable Y con densidad de probabilidad gY , f´acilmente generable. La condici´ on que debe verificar la densidad gY es la siguiente: Existe una constante c tal que f (x) ≤ c, ∀x ∈ Rg(X ). g(x) Se tiene entonces el siguiente algoritmo de generaci´ on mediante la t´ecnica de Aceptaci´ onRechazo: Algoritmo
Paso 1. Generar un valor y0 de la variable Y con funci´ on de densidad gY , mediante el m´etodo de la Transformada Inversa. Paso 2. Generar un n´ umero aleatorio uniforme u. Paso 3. Si se verifica la desigualdad u ≤
f X (y0 ) , cgY (y0 )
10
entonces, X = y 0 . En caso, contrario, volver al Paso 1. Paso 4. Se repiten los Pasos 1, 2 y 3 tantas veces como valores independientes de la variable X se deseen generar. Ejemplo 7.
Seguidamente se muestra un ejemplo, implementado en MatLab, sobre la generaci´ on de una beta con par´ ametros 2 y 4, es decir, sobre la generaci´ on de una variable aleatoria X con densidad de probabilidad f X (x) = 20x(1 − x)3 ,
0 < x < 1,
mediante la t´ecnica de aceptaci´ on-rechazo. La variable aleatoria Y, con la que se compara en esta ocasi´ on, se distribuye seg´ u n una uniforme en el intervalo (0, 1). Derivando el cociente f X (y) gY (y) e igualando a cero se obtiene el valor de la constante c que en este caso viene dada por c = Se tiene entonces que
f X (1/4) 135 = . 64 gY (1/4)
f X (y) 256 = y(1 − y)3 . cgY (y) 27
Se muestra ahora en c´ odigo MatLab la generaci´ on de dicha beta a partir de los c´ alculos anteriores. Semilla s = 20. Escala a = 13. Or´ıgen b = 1. M´odulo de congruencia m = 64. Comparaci´on con una variable aleatoria Y que se distribuye seg´ un una uniforme en el intervalo (0, 1).
11
clear all s = 20; a = 13; b = 1; m = 64; for I = 1 : m n = a ∗ s + b; u = mod(n, m); v = u/m; H = (256/27) ∗ v ∗ (1 − v)3 ; if v <= H s = u; n = a ∗ s + b; u = mod(n, m); v = u/m; Beta = v else I = I + 1; s = u; end end
4.3.
T´ ecnica de Composici´ on
Esta t´ecnica se aplica cuando la densidad de probabilidad de la variable aleatoria X de inter´es se define como la convoluci´ on de dos funciones, la densidad de pesos y la densidades condicionadas que conforman la mixtura continua que define la densidad de la variable aleatoria X. Algoritmo
Paso 1. Generar un valor de la variable aleatoria Y asociada a la densidad de probabilidad de pesos (mediante el m´etodo de la Transformada Inversa). Paso 2. Generar el valor correspondiente de la variable aleatoria X a partir de la densidad condicionada al valor de Y generado (mediante el m´etodo de la Transformada Inversa). Paso 3. Repetir los Pasos 1 y 2 tantas veces como valores de la variable aleatoria X se desee generar.
12
Ejemplo 8.
Seguidamente se muestra un ejemplo, implementado en MatLab, sobre la generaci´ on de una variable aleatoria X con funci´on de densidad ∞
f X (x) = n
y
n
−
exp(−yx)dy.
1
Se considera como densidad de pesos la funci´ on p(y) =
n y n+1
,
1 ≤ y, n ≤ 1,
y como densidades condicionadas h(x/y) = y exp(−yx),
x ≥ 0.
Semilla S = 15. Escala A = 9. Or´ıgen B = 13. M´odulo de congruencia m = 64. Densidad de inter´es, una mixtura de exponenciales con par´ametros dados por los valores de la variable aleatoria generada a partir de la densidad de pesos.
13
clear all S = 15; A = 9; B = 13; m = 64; K = 3; S 2 = 12; A2 = 5; B2 = 15; for I = 1 : m n = A ∗ S + B; u = mod(n, m); V = u/m; if V == 0 S = u; I = I + 1; else DPESO = (K 1/K +1)/(V 1/K +1); end if DPESO == 0 S = u; I = I + 1; else S = u; n = A ∗ S + B; u = mod(n, m); V = u/m; if V == 0 S = u; I = I + 1; else CONDICIONAL = (−1/DPESO) ∗ log(V ) end S = u; end S = u; end
14
5.
Simulaci´ on de vectores aleatorios
En esta secci´on se considera la aplicaci´ on del m´etodo de la transformada inversa para la generaci´ on de vectores aleatorios. Por tanto, el objetivo es generar la variable multidimensional X = (X 1 , . . . , Xn ) con distribuci´ on F X . Se distinguen dos casos: Las variables aleatorias X 1 , . . . , Xn son independientes. Por tanto, su densidad de probabilidad viene dada por n
f X ,...,X
n
1
(x , . . . , x ) = f (x ), n
1
i
i
i=1
siendo f i la densidad marginal de X i , para i = 1, . . . , n . Este caso se reduce al caso unidimensional, pues el vector X se genera a partir de la ecuaci´ on X i = F i 1(U i ), −
siendo
i = 1, . . . , n ,
xi
F i (xi ) =
f i (yi )dyi ,
i = 1, . . . , n ,
−∞
la funciones de distribuci´ on de las componentes X i , i = 1, . . . , n , del vector X, respectivamente. La variable U i se distribuye seg´ un una uniforme sobre el intervalo [0, 1], para i = 1, . . . , n , como es usual en la aplicaci´on del m´etodo de la transformada inversa. Las variables aleatorias X 1 , . . . , Xn son dependientes. Se considera pues el vector (U 1 , . . . , Un ) de componentes aleatorias independientes uniformemente distribuidas sobre el intervalo [0, 1]. El vector X = (X 1 , . . . , Xn ) se genera entonces a partir de la resoluci´ on del siguiente sistema de ecuaciones: F 1 (X 1 ) = U 1 F 2 (X 2 | X 1 ) = U 2 F 3 (X 3 | X 1 , X 2) = U 3 .. . F n (X n | X 1 , . . . , Xn 1 ) = U n . −
Por tanto, el algoritmo que nos permite la programaci´ on del m´etodo se basa u´nicamente en dos pasos:
• Paso 1. Generar n variables independientes uniformemente distrbuidas U (0, 1). on X = (X 1 , . . . , Xn ) del sistema de ecuaciones • Paso 2. Encontrar la soluci´ F 1 (X 1 ) = U 1 F 2 (X 2 | X 1 ) = U 2 F 3 (X 3 | X 1 , X 2 ) = U 3 .. . F n (X n | X 1 , . . . , Xn 1 ) = U n. −
15
Hay n! posibles ordenaciones de las componentes aleatorias X 1 , . . . , Xn del vector on del sistema X. Por tanto, hay n! posibilidades de generar X a partir de la resoluci´ de ecuaciones anterior. Por ejemplo, en el caso n = 2, se tienen 2! = 2 posibilidades dadas por f X ,X (x1 , x2 ) = f 1 (x1 )f 2 (x2 | x 1 ) 1
2
f X ,X (x1 , x2 ) = f 2 (x2 )f 1 (x1 | x 2 ). 1
2
(1) La eficiencia de las simulaciones generalmente depender´a del orden en el cual se seleccionan las variables X 1 , . . . , Xn para conformar el vector X = (X 1, . . . , Xn ).
16
