Simulaci´ on estad´ıstica

26 de enero de 2009

2

Cap´ıtulo 1 Introducci´ on 1.1.

Conceptos b´ asicos

La simulaci´on es la t´ecnica que consiste en realizar experimentos de muestreo sobre el modelo de un sistema. Un modelo no es m´as que un conjunto de variables junto con ecuaciones matem´aticas que las relacionan y restricciones sobre dichas variables. La modelizaci´on es una etapa presente en la mayor parte de los trabajos de investigaci´on (especialmente en las ciencias experimentales). En muchas ocasiones, la realidad es bastante compleja como para ser estudiada directamente y es preferible la formulaci´on de un modelo que contenga las variables m´as relevantes que aparececen en el fen´omeno en estudio y las relaciones m´as importantes entre ellas. Frecuentemente, la resoluci´on de los problemas que se pretenden abordar puede realizarse por procedimientos anal´ıticos sobre el modelo construido (normalmente mediante el uso de herramientas matem´aticas como las de resoluci´on de ecuaciones ordinarias o de ecuaciones diferenciales, el c´alculo de probabilidades, etc.). En otras circunstancias dicha resoluci´on anal´ıtica no es posible (o es tremendamente complicada o costosa) y es preferible una aproximaci´on de la soluci´on mediante simulaci´on.

1.2.

Experimentaci´ on real y simulaci´ on

La experimentaci´on directa sobre la realidad puede tener muchos inconvenientes: un coste muy alto gran lentitud en ocasiones las pruebas son destructivas a veces no es ´etica (experimentaci´on sobre seres humanos) puede resultar imposible (un acontecimiento futuro) 3

4

´ CAP´ITULO 1. INTRODUCCION

Razones como esas (y algunas otras) pueden indicar la ventaja de trabajar con un modelo del sistema real. La estad´ıstica es la ciencia que se preocupa de c´omo estimar los par´ametros y contrastar la validez de un modelo a partir de los datos observados del sistema real que se pretende modelizar. Una vez se ha construido un modelo, la primera tentativa debe ser siempre tratar de resolver anal´ıticamente el problema que nos ocupa. En caso de ser esto posible, la soluci´on es exacta (a menudo la resoluci´on tambi´en es r´apida). En caso contrario puede recurrirse a la simulaci´on que involucrar´a mucha labor de procesado. Gracias a la gran potencia de c´alculo de los computadores actuales los programas de simulaci´on pueden ofrecer una soluci´on aproximada r´apida en la mayor parte de los problemas susceptibles de ser modelados.

Ejemplo 1.2.1 Sup´ongase que se quiere calcular la probabilidad de aparici´ on de exactamente dos caras en tres lanzamientos de una moneda. La experimentaci´on sobre la situaci´ on real consistir´ıa en repetir numerosas veces los tres lanzamientos y observar con qu´e frecuencia se obtienen exactamente dos caras. El sistema real es el mecanismo por el cual se realizan los lanzamientos. Un modelo razonable para este sistema es el de una variable aleatoria X ∈ B (3, 0. 5) (supuesto que la moneda tiene la misma probabilidad de cara que de cruz). Bajo este modelo, se tratar´ıa de calcular P (X = 2). En este caso la resoluci´ on anal´ıtica es factible y muy sencilla    2   1 1 3 3 1 = = 0. 375 P (X = 2) = 2 2 8 2 La simulaci´on consistir´ıa en obtener n´ umeros aleatorios (en ordenador) para replicar artificialmente los tres lanzamientos en gran cantidad de ocasiones, observando la frecuencia relativa con la que aparecen exactamente dos caras. Ejemplo 1.2.2 Sup´ongase el siguiente juego: un jugador lanza una moneda (abonando un euro por cada lanzamiento) hasta que el n´ umero de caras supere en tres al n´ umero de cruces obtenidas. En ese momento el jugador recibe 10 unidades monetarias. ¿Resulta rentable jugar? De nuevo aqu´ı la experimentaci´on real es muy lenta. La modelizaci´on puede realizarse de nuevo gracias a la teor´ıa de la probabilidad. En esta ocasi´ on, sin embargo, la resoluci´on anal´ıtica ser´ıa complicada (salvo que se tengan conocimientos de cadenas de Markov). Parece, por tanto, conveniente una aproximaci´ on mediante simulaci´ on. Se tratar´ıa de ir replicando artificialmente el lanzamiento de la moneda simulando un gran n´ umero de posibles partidas y examinando la p´erdida o ganancia media.

´ 1.3. VENTAJAS E INCONVENIENTES DE LA SIMULACION

1.3.

5

Ventajas e inconvenientes de la simulaci´ on

Ventajas: 1. En casos en los que la resoluci´on anal´ıtica no puede llevarse a cabo. 2. Cuando existen medios de resolver anal´ıticamente el problema pero dicha resoluci´on es complicada y costosa. 3. Si se desea experimentar antes de que exista el sistema. 4. Cuando es imposible experimentar sobre el sistema real por ser dicha experimentaci´on destructiva. 5. En ocasiones en las que la experimentaci´on sobre el sistema es posible pero no ´etica. 6. Es de utilidad en sistemas que evolucionan muy lentamente en el tiempo.

Inconvenientes: 1. La construcci´on de un buen modelo puede ser una tarea muy laboriosa. 2. Frecuentemente el modelo omite variables o relaciones importantes entre ellas. 3. Resulta dif´ıcil conocer la precisi´on de la simulaci´on, especialmente en lo relativo a la precisi´on del modelo formulado.

6

´ CAP´ITULO 1. INTRODUCCION

Cap´ıtulo 2 Generaci´ on de n´ umeros pseudoaleatorios uniformes en (0. 1) 2.1.

Introducci´ on

Casi todos los m´etodos de simulaci´on se basan en la posibilidad de generar n´ umeros aleatorios con distribuci´on U (0. 1). Hasta el gran desarrollo de los ordenadores los n´ umeros aleatorios se obten´ıan por procedimientos experimentales (loter´ıas, ruletas) y se almacenaban en tablas. En la actualidad estos n´ umeros son generados por ordenador y se denominan pseudoaleatorios ya que, en realidad, todos los n´ umeros de la sucesi´on que se genera son predecibles a partir del primero, llamado semilla. En cualquier caso, todo generador de n´ umeros pseudoaleatorios m´ınimamente aceptable debe comportarse como si se tratase de una muestra genuina de datos independientes de una U (0. 1).

2.1.1.

Propiedades deseables de un generador de n´ umeros pseudoaleatorios

Para poder utilizar sin reservas un generador de n´ umeros pseudoaleatorio ´este debe satisfacer los contrastes estad´ısticos m´as habituales en este contexto: los de aleatoriedad (los contrastes de rachas o los de saltos), los de independencia (como los basados en autocorrelaciones, el test de Ljung-Box, el contraste de pares seriados, etc) y los de bondad de ajuste a una U (0. 1) (entre ellos es test chi-cuadrado y el de KolmogorovSmirnov). Tambi´en existen otros contrastes espec´ıficos que tratan de indagar a la vez sobre varios de los aspectos anteriores. Entre ellos destacamos el contraste del poker y el del coleccionista.

7

´ DE NUMEROS ´ CAP´ITULO 2. GENERACION UNIFORMES EN (0. 1)

8

Adem´as de estas propiedades de tipo estad´ıstico existen otros requisitos computacionales. Unos y otros pueden resumirse en la siguiente lista. Requisitos deseables para un generador 1. Producir muestras seg´ un una distribuci´on U (0. 1). 2. Pasar los contrastes de aleatoriedad e independencia m´as habituales. 3. Que la sucesi´on generada sea reproducible a partir de la semilla. 4. Tener una longitud de ciclo tan grande como se desee. 5. Generar valores a alta velocidad. 6. Ocupar poca memoria.

2.2.

M´ etodo de los cuadrados medios

Es debido a von Neumann y tiene fundamentalmente s´olo inter´es hist´orico. 1. Se toma un n´ umero entero inicial, x0 , llamado semilla, de 2n cifras. 2. Se eleva al cuadrado obteniendo un n´ umero de 4n cifras (completando, quiz´a, con ceros a la izquierda). 3. Se considera x1 el n´ umero entero formado por las 2n cifras centrales. 4. Se eleva al cuadrado x1 y se repite el proceso anterior tantas veces como sea preciso. 5. Finalmente se consideran los n´ umeros ui =

xi , 102n

ya en el intervalo (0. 1).

Ejemplo 2.2.1 T´omese n = 2 y x0 = 4122. Resulta: x0 = 4122 x20 = 16|9908|84 x1 = 9908 x21 = 98|1684|64 x2 = 1684 x22 = 02|8358|56 x3 = 8358 x23 = 69|8561|64 x4 = 8561 x24 = 73|2907|21 x5 = 2907 x25 = 08|4506|49 De esta forma, los n´ umeros pseudoaleatorios en (0. 1) son u0 = 0. 4122 u1 = 0. 9908 u2 = 0. 1684 u3 = 0. 8385 u4 = 0. 8561 u5 = 0. 2907 Siguiendo, de nuevo, con n = 2, pero tomando como semilla x0 = 3708, se obtiene x0 x2 x4 x6 x8

= 3708 = 1300 = 6100 = 4100 = 6100

x20 x22 x24 x26

= 13|7492|64 = 01|6900|00 = 37|2100|00 = 16|8100|00

x1 x3 x5 x7

= 7292 = 6900 = 2100 = 8100

x21 x23 x25 x27

= 56|1300|64 = 47|6100|00 = 04|4100|00 = 65|6100|00

As´ı pues, como x8 = x4 , los valores u4, u5 , u6 , u7 se repetir´ an c´ıclicamente de forma indefinida. Este tipo de fen´omenos de ciclo corto son su mayor inconveniente.

´ 2.3. METODO DE LEHMER

2.3.

9

M´ etodo de Lehmer

El m´etodo consiste en los siguientes pasos: 1. Se toma como semilla un n´ umero entero, x0 , de n cifras. 2. Se elige otro entero, c, de k cifras. Suele tomarse k < n. 3. Se calcula x0 · c, n´ umero de, a lo sumo, n + k cifras. 4. Se separan las k cifras de la izquierda de x0 · c y al n´ umero formado por las n cifras restantes se le resta el que forman esas k cifras de la izquierda, dando lugar a x1 . 5. Se repite este proceso tantas veces como sea necesario. 6. Se devuelven los valores ui =

xi . 102n

Ejemplo 2.3.1 Tomando n = 4, k = 2, x0 = 4122 y c = 76, se obtiene x0 x1 x2 x3 x4 x5

= 4122 = 3241 = 6292 = 8145 = 8959 = 0816

x0 · c = 31|3272 x1 · c = 24|6316 x2 · c = 47|8192 x3 · c = 61|9020 x4 · c = 68|0884 x5 · c = 06|2016

3272 − 31 = 3241 6316 − 24 = 6292 8192 − 47 = 8145 9020 − 61 = 8959 0884 − 68 = 0816 2016 − 06 = 2010

De esta forma u0 = 0. 4122 u1 = 0. 3241 u2 = 0. 6292 u3 = 0. 8145 u4 = 0. 8959 u5 = 0. 0816 Todav´ıa en el caso de que n = 4 y k = 2, pero con x0 = 2000 y c = 50, se tiene x0 · c = 10|0000 y as´ı x1 = 0000 − 10 = −10 < 0. Este es precisamente uno de los peores inconvenientes de este m´etodo: la aparici´ on de iterantes negativos. Tambi´en aparecen, con frecuencia, ciclos cortos (en particular, el cero es un valor absorbente de este generador).

´ DE NUMEROS ´ CAP´ITULO 2. GENERACION UNIFORMES EN (0. 1)

10

2.4.

M´ etodos congruenciales

Se basan en la idea de considerar una combinaci´on lineal de los u ´ltimos k enteros generados y calcular su resto al dividir por un entero fijo m. El m´etodo congruencial lineal simple que procede como sigue: 1. Elegir un n´ umero entero positivo m (normalmente en relaci´on con el tipo de enteros que se va a usar) y otros dos n´ umeros enteros, a y c, tales que 0 < a < m y 0 ≤ c < m. 2. Fijar la semilla x0 , un valor entero inicial que cumpla 0 ≤ x0 < m. 3. Obtener de forma recurrente xn = (axn−1 + c) mod m para n = 1, 2, . . . 4. Devolver los valores un =

xn , m

n = 0. 1, . . .

Cuando tomamos c = 0 el generador se dice congruencial multiplicativo. Si c > 0, se dice congruencial mixto. Ejemplo 2.4.1 Consid´erese un generador congruencial con m = 8, a = 5, c = 4: xn = (5xn−1 + 4) mod 8 Tomando como semilla los valores 5 ´o 2 se obtiene: x0 = 5 x1 = 5 x2 = 5 · · · x0 = 2 x1 = 6 x2 = 2 x3 = 6 · · · que presentan ciclos de longitud 1 y 2 respectivamente. Cambiando el valor de c a 2 se tiene xn = (5xn−1 + 2) mod 8 y as´ı, x0 = 5 x1 = 3 x2 = 1 x3 = 7 x4 = 5 · · · x0 = 2 x1 = 4 x2 = 6 x3 = 0 x4 = 2 · · · donde ambos ciclos son de longitud cuatro. Finalmente dejando el mismo valor de m pero eligiendo a = 5 y c = 5, se tiene xn = (5xn−1 + 5) mod 8, que conduce a x0 x5 x0 x5

=5 =2 =2 =3

x1 x6 x1 x6

=6 =7 =7 =4

x2 x7 x2 x7

=3 =0 =0 =1

x3 x8 x3 x8

=4 =5 =5 =2

con ciclo de longitud 8, que es el m´aximo valor posible.

x4 = 1 ··· x4 = 6 ···

´ 2.4. METODOS CONGRUENCIALES

2.4.1.

11

Generadores congruenciales de ciclo m´ aximo

Adem´as de las propiedades estad´ısticas deseables para cualquier generador, una cuesti´on importante para los generadores congruenciales (como se ha visto en los ejemplos previos) es la de garantizar que el ciclo del generador sea m´aximo (o, cuando menos, muy elevado). En la pr´actica tratar´a de tomarse el ciclo igual o muy pr´oximo al n´ umero de enteros de tipo largo del lenguaje en cuesti´on. En general, se define la longitud del ciclo (o per´ıodo) de un generador de n´ umeros pseudoaleatorios uniformes, como el menor n´ umero entero positivo, p, que cumple que existe un n0 natural tal que xi+p = xi para todo i = n0 . n0 + 1, . . . En el caso de un generador congruencial mixto el m´aximo valor para el per´ıodo es m. Tambi´en puede demostrarse que si un generador congruencial tiene ciclo m´aximo para cierta elecci´on de la semilla, entonces lo tiene para cualquier otra. Un resultado que sirve para caracterizar qu´e propiedades deben cumplir los par´ametros de un generador congruencial para que tenga per´ıodo m´aximo es el teorema de Knuth (1969). Teorema 2.4.1 (Knuth) Las siguientes condiciones son necesarias y suficientes para que un generador congruencial con par´ ametros m, a y c, tenga per´ıodo m´ aximo (i.e. p = m). 1. c y m son primos entre s´ı (i.e. m.c.d. (c, m) = 1). 2. a − 1 es m´ ultiplo de todos los factores primos de m (i.e. a ≡ 1mod g, para todo g factor primo de m). 3. Si m es m´ ultiplo de 4, entonces a − 1 tambi´en lo ha de ser (i.e. m ≡ 0 mod 4 ⇒ a ≡ 1 mod 4). A la luz del teorema de Knuth, es f´acil darse cuenta porqu´e s´olo el tercero de los generadores del ejemplo anterior ten´ıa per´ıodo ´optimo.

´ DE NUMEROS ´ CAP´ITULO 2. GENERACION UNIFORMES EN (0. 1)

12

2.4.2.

Generadores congruenciales de algunos lenguajes y bibliotecas de rutinas

En ordenadores binarios es muy com´ un elegir m = 2β o m = 2β − 1, donde β depende del tama˜ no de palabra (t´ıpicamente m ser´a el mayor entero representable en el ordenador o una unidad mayor que ´el). En los generadores con m = 2β resulta especialmente f´acil expresar las condiciones del teorema de Knuth de forma mucho m´as sencilla: As´ı, al ser m una potencia de 2, su u ´nico factor primo es el 2 y, por tanto la primera condici´on equivale a que c sea impar. Para β ≥ 2 se tiene que m es m´ ultiplo de 4 y, por tanto la tercera condici´on impone que a − 1 tambi´en lo sea. Por u ´ltimo, de nuevo por ser el 2 el u ´nico factor primo de m, la segunda condici´on pedir´ıa que a − 1 fuese par, lo cual ya es consecuencia de que sea m´ ultiplo de 4. β En resumen, si m = 2 , con β ≥ 2, el generador congruencial tiene per´ıodo m´aximo si y s´olo s´ı c es impar y a = 4k + 1, siendo k un n´ umero natural. Algunos generadores habituales en lenguajes con enteros (con signo) de 36 bits corresponden con las elecciones m = 235 m = 235

a = 27 + 1 c = 1 a = 515 c=1

Todos ellos tienen tienen per´ıodo m´aximo (e igual a 235 ' 3,44 × 1010 ). Otros generadores congruenciales para enteros (con o sin signo) de 32 bits y algunos lenguaje o bibliotecas que los usan o los han usado en el pasado son m = 231 m = 231 − 1 m = 231 − 1 m = 232

a = 314159269 a = 16807 a = 630360016 a = 663608941

c = 453805245 c=0 APL, IMSL y SIMPL/I c=0 Algunos FORTRAN c=0 Ahrens y Dieter (1974)

Aunque s´olo el primero tiene per´ıodo m´aximo, los dem´as lo tienen muy elevado. El lenguaje C (bajo UNIX) posee un generador congruencial de n´ umeros pseudoa48 leatorios de 48 bits: el drand48. Sus par´ametros son m = 2 , a = 25214903917 y c = 11. La semilla se establece mediante la sentencia srand48() introduciendo como argumento un entero de 32 bits que corresponde a los 32 bits m´as significativos de x0 (entero de 48 bits). Los 16 bits de menor orden se toman siempre coincidentes con el n´ umero (decimal) 13070. Los par´ametros a y c se pueden modificar por el usuario desde otras rutinas del C. Los valores por defecto para estas cantidades ofrecen un generador de per´ıodo m´aximo ya que m = 2β , con β = 48, c es impar y a = 6 303 725 979 · 4 + 1.

2.5. MEDIDAS ESTAD´ISTICAS DE LA CALIDAD DE UN GENERADOR

2.4.3.

13

Otros generadores congruenciales

Generador congruencial lineal m´ ultiple xn = (a1 xn−1 + a2 xn−2 + · · · + ak xn−k + c) mod m Generadores congruenciales no lineales xn = g(xn−1 ) mod m Knuth(1981): xn = g(xn−1 ) mod m, con g(x) = ax2 + bx + c. Generadores congruenciales matriciales xn = (Axn−1 + C) mod m con xn un vector d−dimensional y A y C matrices d × d. Los elementos de los vectores y de las matrices son enteros entre 1 y m − 1.

2.5.

Medidas estad´ısticas de la calidad de un generador de n´ umeros pseudoaleatorios

La mayor´ıa de los contrastes estad´ısticos para estudiar la calidad de un generador de n´ umeros aleatorios se basan en medir posibles discrepancias (en alg´ un sentido) de las muestras generadas por el m´etodo en cuesti´on con respecto a las hip´otesis de aleatoriedad, independencia y ajuste a una distribuci´on U (0. 1). En casi todos los casos se acaba recurriendo a un contraste de tipo chi-cuadrado en el que se comparan las frecuencias esperadas, ei , de ciertas modalidades, con las observadas, oi , mediante el estad´ıstico D=

k X (oi − ei )2 i=1

ei

,

que seguir´a, aproximadamente, una distribuci´on χ2k−1 , bajo la hip´otesis nula que se contrasta. A continuaci´on detallamos algunos de los contrastes m´as habituales.

2.5.1.

Contraste chi-cuadrado de Pearson

El test chi-cuadrado est´a basado en un estad´ıstico que, para una variable discreta o de tipo cualitativo, compara la frecuencia de aparici´on de cada una de las modalidades observadas (ni ) con las frecuencias esperadas, en base a la distribuci´on de probabilidad especificada (ei ). Concretamente, para una variable discreta con k modalidades, el contraste se basa en el estad´ıstico sugerido por Pearson en el a˜ no 1900: Q=

k X (ni − ei )2 i=1

ei

14

´ DE NUMEROS ´ CAP´ITULO 2. GENERACION UNIFORMES EN (0. 1)

cuya distribuci´on aproximada es la de una χ2k−1 , siempre que la distribuci´on especificada sea la correcta. Comentarios: 1. Es muy corriente aplicar el test chi-cuadrado a´ un en casos en los que la distribuci´on de la variable no est´a totalmente especificada sino que depende de alg´ un par´ametro que, por tanto, habr´a de ser estimado (por ejemplo el caso en que se suponga que la variable en concreto sigue una distribuci´on de Poisson y resta por especificar su par´ametro λ). En estos casos la distribuci´on aproximada del test ha de ser corregida umero de para incorporar esta informaci´on pasando a ser una χ2k−r−1 , siendo r el n´ par´ametros estimados por m´axima verosimilitud. 2. El contraste chi-cuadrado se utiliza habitualmente incluso cuando la variable objeto de estudio es continua. En este caso, dicha variable ha de ser agrupada en intervalos de clase. 3. Una limitaci´on, bastante recomendable en la pr´actica, es la de no llevar a cabo el contraste chi-cuadrado cuando la frecuencia esperada de alguna clase sea menor que 5. Aquellos casos en que esta condici´on falla pueden tratarse agrupando varios valores distintos hasta que se cumpla esta restricci´on. Ejemplo 2.5.1 Consid´erese los siguientes datos sobre “espacio de disco (en MBytes) ocupado por usuario en una estaci´on de trabajo”: 35, 45, 47, 50. 31, 30. 25, 33, 35, 40. 45, 47, 49, 42, 40. 50. 46, 55, 42, 46. Est´ udiese la bondad de ajuste de dichos datos a una distribuci´ on U [a, b], realizando el contraste chi-cuadrado. En primer lugar se estiman por m´axima verosimilitud los par´ ametros de la distribuci´ on uniforme, obteni´endose respectivamente el m´ınimo y el m´ aximo muestral: a ˆ = 25 ˆ y b = 55. Dividiendo el intervalo [25, 55] en cuatro clases de igual longitud se obtiene la siguiente tabla: clases ni [25, 320 5) 3 [320 5, 40) 5 [40. 470 5) 8 [470 5, 55] 4 Total 20

ei (ni − ei )2 /ei 5 00 8 5 0 5 10 8 5 00 2 20 20 8

En este caso el valor del estad´ıstico es Q = 20 8, que corresponde a un nivel cr´ıtico de p = 00 09426, para una distribuci´on chi-cuadrado con 4 − 2 − 1 = 1 grado de libertad, con lo cual aunque puede aceptarse la hip´ otesis de que la distribuci´ on poblacional es 0 0 0 uniforme con α = 0 01 o α = 0 05, no se aceptar´ıa con α = 0 1.

2.5. MEDIDAS ESTAD´ISTICAS DE LA CALIDAD DE UN GENERADOR

2.5.2.

15

Contraste de Kolmogorov-Smirnov

A diferencia del procedimiento anterior, el test de Kolmogorov-Smirnov est´a especialmente dise˜ nado para el contraste de ajuste a distribuciones continuas. El contraste de Kolmogorov-Smirnov est´a basado en la distribuci´on del estad´ıstico Dn : Dn = sup |Fn (x) − F (x)| x∈R

que representa la m´axima discrepancia, en vertical, entre la funci´on de distribuci´on emp´ırica N´ umero de observaciones Xi ≤ x Fn (x) = n y la te´orica (F ). Dn = m´ax Dn,i = m´ax{|Fn (x(i) ) − F (x(i) )|, |Fn− (x(i) ) − F (x(i) )|} i=1,2,...,n

Ejemplo 2.5.2 Contrastar que la muestra de datos de ocupaci´ on de disco provenga de una distribuci´on N (40. 3). x(i) 25 30 31 33 35 40 42 45 46 47 49 50 55

F (x(i) ) Fn (x(i) ) Fn− (x(i) ) 0 00 05 0 0 0 0 0 00043 0 10 0 05 0 0 0 00135 0 15 00 10 00 00982 00 20 00 15 00 04779 00 30 00 20 0 0 05 0 40 00 30 00 74751 00 50 00 40 00 95221 00 60 00 50 0 0 0 97725 0 70 00 60 00 99019 00 80 00 70 0 0 0 99865 0 85 00 80 00 99957 00 95 00 85 1 1 00 95

Dn,i 00 05 00 09957 00 14865 00 19018 00 25221 00 2 00 34751 00 45221 00 37725 00 29019 00 19865 00 14957 00 05

Se obtiene Dn = 00 45221, cuyo nivel cr´ıtico es p < 00 01, que indica un claro rechazo de la hip´otesis establecida.

2.5.3.

Contrastes de independencia

Una de las hip´otesis estad´ısticas m´as importantes asumidas en toda la inferencia param´etrica cl´asica es que la muestra no es m´as que una realizaci´on de un conjunto de variables independientes. M´ etodos gr´ aficos Como paso previo a la utilizaci´on de m´etodos cuantitativos siempre ser´a ilustrativo representar gr´aficamente los datos frente al tiempo, es decir, secuencialmente seg´ un su orden de recogida (Xi frente al ´ındice i).

16

´ DE NUMEROS ´ CAP´ITULO 2. GENERACION UNIFORMES EN (0. 1)

Las situaciones de dependencia se corresponden con gr´aficas demasiado estables o demasiado cambiantes, es decir, o bien muy planas, o bien en continuos dientes de sierra o zigzagueantes. Por el contrario, los casos de ausencia de dependencia se identificar´an mediante gr´aficas moderadamente cambiantes. Otro tipo de gr´aficas que tambi´en permite detectar situaciones de dependencia son las gr´aficas en las que se representa cada valor de la muestra frente al valor anteriormente observado (Xi+1 frente a Xi ). Otra forma de estudiar si hay dependencia en la muestra es mediante el c´alculo de los coeficientes de autocorrelaci´on. Del mismo modo que el coeficiente de correlaci´on lineal mide, de alguna forma, la presencia o ausencia de dependencia entre las variables en cuesti´on, pueden definirse toda una serie de coeficientes de autocorrelaci´ on que, en ese mismo sentido, miden la dependencia entre los datos observados con cierto n´ umero de instantes de diferencia. Dada una muestra X1 , X2 , . . . , Xn , se define el coeficiente de autocorrelaci´on muestral de orden uno como el valor Pn (Xi − X)(Xi−1 − X) r(1) = i=2 Pn 2 i=1 (Xi − X) que expresa, de alguna manera, la correlaci´on entre cada dato y el observado un instante antes. En general, el coeficiente de autocorrelaci´on de orden k (o a k retardos), se define por Pn r(k) =

i=k+1 (Xi − X)(Xi−k Pn 2 i=1 (Xi − X)

− X)

Una forma gr´afica de visualizar los coeficientes de autocorrelaci´on es el llamado correlograma. En ´el se representan, para los primeros retardos, unas barras con altura igual al coeficiente de autocorrelaci´on correspondiente al retardo en cuesti´on. As´ı pues, las barras del correlograma oscilan entre −1 y 1. Es obvio que la dependencia positiva se caracterizar´a por la presencia de muchos coeficientes de correlaci´on sustancialmente mayores que cero, mientras que la negativa vendr´a dada por autocorrelaciones de signo alternado, siendo la del primer retardo negativa. El caso de independencia se corresponde con autocorrelaciones cercanas a cero. Bajo la hip´otesis de independencia, cada coeficiente de autocorrelaci´on muestral, r(k), tiene distribuci´on l´ımite normal de media cero y varianza 1/n. Este hecho permite, por s´ı mismo, el contrastar la hip´otesis sobre si el coeficiente de autocorrelaci´on te´orico es cero. √De hecho, es habitual incluir en el correlograma bandas de confianza a una distancia 2/ n del eje horizontal, de manera que se considerar´an siginificativamente distintos de cero aquellas autocorrelaciones que sobresalgan de los l´ımites. Otra posibilidad es incluir bandas de amplitud variable, teniendo en cuenta que, en el supuesto de que sean no nulas las primeras k − 1 autocorrelaciones (ρ1 , . . . , ρk−1 ), la P 2 1 + k−1 j=1 ρj varianza de r(k) es . n

2.5. MEDIDAS ESTAD´ISTICAS DE LA CALIDAD DE UN GENERADOR

17

Contrastes basados en rachas Pensemos en una muestra de una variable con dos posibles resultados (por ejemplo ESSSEESEESES, con E: “error en un dispositivo” y S: “dispositivo sin error”). Definici´ on 2.5.1 Una racha es una sucesi´ on de valores consecutivos repetidos que est´e flanqueada por valores adyacentes distintos. En el ejemplo anterior las rachas ser´ıan E, SSS, EE, S, EE, S, E, S. Independientemente de lo probable que sea observar los valores concretos de la variable, es obvio que el n´ umero total de rachas (o las longitudes de las mismas) constituye una medida de lo aleatoriamente que est´an repartidos los posibles valores en cuesti´on a lo largo de la muestra observada. Demasiadas rachas implican excesiva alternancia de valores (dependencia negativa), mientras que pocas rachas indican largas sucesiones de valores contiguos repetidos (dependencia positiva). El contraste del n´ umero total de rachas Consid´erese una muestra de tama˜ no n correspondiente a una variable con dos posibles resultados, de manera que existen n1 observaciones de un tipo y n2 iguales al otro valor de la variable en cuesti´on (n1 + n2 = n). Si llamamos R al n´ umero total de rachas observadas en la muestra, pueden obtenerse la media y la varianza de esta variable aleatoria (condicionadas a que hay n1 y n2 de elementos de cada tipo): E(R) = 1 + V ar(R) =

2n1 n2 n

2n1 n2 (2n1 n2 − n) n2 (n − 1)

Cuando el tama˜ no muestral n tiende a infinito, de forma que adem´as n1 /n tienda a una constante, la distribuci´ on de la variable R se puede aproximar por la distribuci´on   p N E(R), V ar(R) . Aunque originalmente el test de las rachas est´a dise˜ nado para una distribuci´on con s´olo dos posibles valores, suele aplicarse a aquellos casos en los que la distribuci´on en cuesti´on es continua codificando las observaciones con los valores + o −, seg´ un el dato en cuesti´on quede por arriba o por abajo de la mediana muestral.

18

´ DE NUMEROS ´ CAP´ITULO 2. GENERACION UNIFORMES EN (0. 1)

El contraste de rachas ascendentes y descendentes Cuando la variable en cuesti´on presenta una distribuci´on de tipo continuo, a pesar de que el test de las rachas por encima o por debajo de la mediana se puede usar, existe una contraste, basado en el n´ umero de cierto tipo de rachas, que trata de hacer un uso m´as intensivo de la continuidad de la variable. Se trata del contraste basado en el n´ umero total de rachas ascendentes o descendentes: Para cada par de datos consecutivos se anota un signo + si est´an en orden ascendente y − si est´an en orden descendente. Con el conjunto de datos (n en total) se consigue formar una tira de n − 1 signos + o −, sobre los cuales se cuenta el n´ umero total de rachas R. La distribuci´on del estad´ıstico R est´a tabulada para tama˜ nos muestrales pequen ˜os (normalmente para n < 25), mientras que, cuando el n´ umero de datos, n, excede de los valores tabulados, puede usarse una aproximaci´ on por la distribuci´on   p N (2n − 1)/3, (16n − 29)/90 . El contraste de Ljung-Box Existen distintos procedimientos para contrastar la independencia mediante la realizaci´on de un contraste sobre si los primeros m coeficientes de autocorrelaci´on son cero. Uno de los m´as utilizados a tal efecto es el de Ljung-Box, que utiliza para ello el estad´ıstico m X r(k)2 Q = n(n + 2) n−k k=1 que se distribuye, aproximadamente, seg´ un una χ2m−1 . Ejemplo 2.5.3 Consider´ense los datos de “espacio de disco duro ocupado por usuario”: 35, 45, 47, 50. 31, 30. 25, 33, 35, 40. 45, 47, 49, 42, 40. 50. 46, 55, 42, 46. Contrastar la independencia usando los test de las rachas y el de Ljung-Box. Test del n´ umero de rachas Dada la muestra original, se clasifican ahora los datos de la misma en funci´ on de 0 que se hallen por encima o por debajo de la mediana, Me = (42 + 45)/2 = 43 5 : − + + + − − − − − − + + + − − + + + −+ El n´ umero de rachas es R = 8 y n1 = n2 = 10. En las tablas puede encontrarse el p-valor: p = 00 256. Por tanto, se acepta la independencia.

2.5. MEDIDAS ESTAD´ISTICAS DE LA CALIDAD DE UN GENERADOR

19

Test de las rachas ascendentes y descendentes La secuncia de signos + y - correspondiente a los datos de la muestra es: + + + − − − + + + + + + − − + − + − + El n´ umero de rachas es R = 9. Bajo independencia, la probabilidad de que el n´ umero 0 de rachas ascendentes y descendentes sea menor o igual que 9 es 0 0255. De aqu´ı se deduce que el nivel cr´ıtico del contraste es p = 00 051. Este valor plantea serias dudas sobre la hip´otesis de independencia: tendr´ıamos una aceptaci´ on muy justa de dicha 0 0 hip´ otesis con nivel α = 0 05 pero ya no con α = 0 06. Test de Ljung-Box Tomemos m = 4. Los cuatro primeros coeficientes de autocorrelaci´ on son r(1) = 0 0 0 0 0 52550. r(2) = 0 33034, r(3) = −0 09887 y r(4) = −0 06706. A partir de estos valores, el estad´ıstico de Ljung-Box resulta Q = 90 44, cuyo pvalor est´a comprendido entre 00 01 y 00 025, como se puede ver en las tablas de una chi-cuadrado con tres grados de libertad. En resumen, parece que existen razones para rechazar la independencia.

2.5.4.

El contraste del coleccionista

Este procedimiento requiere fijar un entero positivo, M , y discretizar los valores generados, X1 , X2 , . . . , Xn , de la forma dM · Xi e+1, donde dxe denota la parte entera de x. As´ı se consigue una sucesi´on de enteros aleatorios cuyos valores est´an comprendidos entre 1 y M . Ahora se procede (como un coleccionista) a contabilizar cu´al es el n´ umero, Q, (aleatorio) de valores a generar hasta que se completa la colecci´on de todos los enteros entre 1 y M . Obviamente, bajo las hip´otesis de aleatoriedad y distribuci´on U (0. 1), cada posible entero entre 1 y M tiene la misma probabilidad de aparecer en cada generaci´on y, por tanto, resulta posible calcular la distribuci´on de probabilidad de Q. De esta forma podemos utilizar los valores calculados de las probabilidades P (Q = M ) , P (Q = M + 1) , . . . para calcular las frecuencias esperadas de cada clase y confrontarlas con las observadas v´ıa el estad´ıstico chi-cuadrado.

´ DE NUMEROS ´ CAP´ITULO 2. GENERACION UNIFORMES EN (0. 1)

20

2.5.5.

Contrastes de salto

Dados dos n´ umeros α y β tales que 0 ≤ α < β ≤ 1, los contrastes de saltos tratan de examinar, para cada valor generado, Xi , si se cumple α ≤ Xi ≤ β, anotando, en ese caso, un 1 (0 en caso contrario). En estas condiciones, la probabilidad de que aparezca un 1 es p = β − α y la de que aparezcan j ceros desde la aparici´on de un uno hasta la del siguiente uno es pj = p (1 − p)j , j = 0. 1, 2, . . . (que corresponde a una distribuci´on geom´etrica). De nuevo puede aplicarse el test chi-cuadrado a las clases resultantes. Las elecciones m´as habituales de α y β dan lugar a los siguientes contrastes: El test de rachas bajo la mediana Consiste en tomar α = 0 y β = 1/2. El test de rachas sobre la mediana Corresponde al caso α = 1/2 y β = 1. El test del tercio medio Que no es m´as que la elecci´on α = 1/3 y β = 2/3.

2.5.6.

El contraste de permutaciones

Dada la sucesi´on de n´ umeros pseudoaleatorios generada se consideran bloques de T valores consecutivos. Cada uno de los bloques puede presentar una cualquiera de las T ! posibles ordenaciones de esos T valores. Adem´as, de ser el generador adecuado, cada posible ordenaci´on ocurrir´a con igual probabilidad: T1! . El test consiste en observar una gran n´ umero de bloques y comparar las frecuencias observadas de cada posible ordenaci´on con las esperadas mediante el test chi-cuadrado. Las elecciones m´as comunes son T = 3, 4, ´o 5.

2.5.7.

El contraste del poker

En un primer momento se procede como en el contraste del coleccionista con M = 10. A partir de aqu´ı hay varias formas de actuar: Poker 1 Se toman conjuntos sucesivos de cinco enteros y, para cada uno, se determina cu´al de las siguientes posibilidades se da: 1. Un mismo entero se repite cinco veces (abreviadamente, AAAAA). 2. Un mismo entero se repite cuatro veces y otro distinto aparece una vez (AAAAB). 3. Un entero se repite tres veces y otro distinto se repite dos (AAABB). 4. Un entero se repite tres veces y otros dos distintos aparecen una vez cada uno (AAABC).

2.5. MEDIDAS ESTAD´ISTICAS DE LA CALIDAD DE UN GENERADOR

21

5. Un entero se repite dos veces, otro distinto se repite tambi´en dos veces y un tercer entero diferente aparece una s´ola vez (AABBC). 6. Un entero se repite dos veces y otros tres distintos aparecen una vez cada uno (AABCD). 7. Los cinco enteros que aparecen son todos distintos (ABCDE). Bajo la hip´otesis de aleatoriedad y ajuste a una U (0. 1), pueden calcularse las probabilidades de estas modalidades: P (AAAAA) = 0. 0001, P (AAAAB) = 0. 0045, P (AAABB) = 0. 0090. P (AAABC) = 0. 0720. P (AABBC) = 0. 1080. P (AABCD) = 0. 5040. P (ABCDE) = 0. 3024. Es frecuente que las clases AAAAA y AAAAB se agrupen a la hora de aplicar el test chi-cuadrado, ya que, en caso contrario, la restricci´on habitual ei ≥ 5 llevar´ıa a que 0. 0001 · n5 ≥ 5, es decir, n ≥ 250 000. Poker 2 Algo tambi´en bastante habitual es usar conjuntos de cinco enteros (como en el caso anterior) pero definiendo las categor´ıas seg´ un el n´ umero de enteros distintos de entre los cinco observados. As´ı P (1 entero diferente) = 0. 0001, P (2 enteros diferentes) = 0. 0135, P (3 enteros diferentes) = 0. 1800. P (4 enteros diferentes) = 0. 5040. P (5 enteros diferentes) = 0. 3024, procediendo frecuentemente a agrupar las dos primeras modalidades. Poker 3 A menudo se consideran conjuntos de cuatro enteros. En tal caso, P (AAAA) = 0. 001, P (AAAB) = 0. 036, P (AABB) = 0. 027, P (AABC) = 0. 432, P (ABCD) = 0. 504, siendo tambi´en bastante habitual el agrupar las dos primeras categor´ıas.

22

2.5.8.

´ DE NUMEROS ´ CAP´ITULO 2. GENERACION UNIFORMES EN (0. 1)

El contraste de pares seriados

La idea consiste en fijar un entero M ≥ 2 y considerar los enteros dM · Xi e + 1, tomar ahora estos valores apareados y utilizar el contrate chi-cuadrado considerando como categor´ıas los posibles pares (i, j) tales que i, j ∈ {1, 2, . . . , M }. As´ı se medir´a la discrepancia entre la frecuencias observadas en estas categor´ıas y las esperadas, iguales todas a n2 M12 . La elecciones m´as frecuentes son M = 3, 10 ´o 20.

2.5.9.

Chi-cuadrado sobre chi-cuadrado

Todos los contrastes anteriores se han planteado desde la perspectiva de la realizaci´on de una u ´nica prueba. Es decir, se toma un n´ umero, n (normalmente grande), de valores obtenidos por el generador y se realiza el contraste evaluando el estad´ıstico y compar´andolo con el punto cr´ıtico de una chi-cuadrado para decidir si se acepta o rechaza la hip´otesis (independencia, ajuste, aleatoriedad). En realidad tiene mucho m´as sentido la realizaci´on de un gran n´ umero de pruebas, evaluando en cada una el valor del estad´ıstico y, o bien observar que la proporci´on de rechazos del test se aproxima al valor nominal fijado (normalmente α = 0. 01 ´o α = 0. 05), o m´as precisamente aplicando, de nuevo, el contraste chi cuadrado para comprobar el ajuste de la distribuci´on del estad´ıstico a la chi-cuadrado especificada bajo la hip´otesis nula.

Cap´ıtulo 3 M´ etodos universales para la simulaci´ on de variables continuas En lo que sigue se expondr´an dos de los m´etodos generales para simular distribuciones continuas: el m´etodo de inversi´on y el de aceptaci´on/rechazo. Ambos son aplicables a gran n´ umero de contextos siempre que la distribuci´on que se desea simular tenga ciertas caracter´ısticas. En ambos casos la herramienta indispensable es alg´ un m´etodo de generaci´on de n´ umeros pseudoaleatorios uniformes en (0,1).

3.1.

El m´ etodo de inversi´ on

Es el m´etodo universal por antonomasia para simular distribuciones continuas. Tambi´en a veces se denota m´etodo de Montecarlo. Est´a basado en el siguiente resultado te´orico. Teorema 3.1.1 (de inversi´ on) Sea X una variable aleatoria con funci´ on de distribuci´ on F , continua e invertible. Entonces, la variable aleatoria U = F (X), transformada ˙ de la original mediante su propia funci´ on de distribuci´ on, tiene distribuci´ on U (0, 1). Rec´ıprocamente, si U ∈ U (0, 1) entonces la variable F −1 (U ) tiene funci´ on de distribuci´on F (la misma distribuci´ on que la de X). Demostraci´ on: Denotando por G la funci´on de distribuci´on de U y dado un valor u ∈ (0, 1), se tiene

G (u) = P (U ≤ u) = P (F (X) ≤ u) = P X ≤ F −1 (u) = F F −1 (u) = u Por otra parte es obvio que G (u) = 0 si u ≤ 0 y G (u) = 1 si u ≥ 1, con lo cual G es la funci´on de distribuci´on de una U (0, 1). Para la segunda parte, denotando por H la funci´on de distribuci´on de X = F −1 (U ), con U (0, 1),
H (x) = P (X ≤ x) = P F −1 (U ) ≤ x = P (U ≤ F (x)) = F (x)

23

´ 24 CAP´ITULO 3. METODOS UNIVERSALES PARA VARIABLES CONTINUAS El resultado anterior da pie al siguiente algoritmo gen´erico para simular cualquier variable continua con funci´on de distribuci´on F invertible: Algoritmo (m´ etodo de inversi´ on) 1. Generar U ∼ U (0, 1). 2. Devolver X = F −1 (U ). Ejemplo 3.1.1 Dar un algoritmo, basado en el m´etodo de inversi´ on, para simular la distribuci´ on exponencial de par´ametro λ > 0. La funci´ on de densidad de una exp (λ) es  λe−λx si x ≥ 0 f (x) = 0 si x < 0 y su funci´ on de distribuci´on:  F (x) =

1 − e−λx 0

si x ≥ 0 si x < 0

que es continua e invertible en el intervalo [0. ∞). Obs´ervese que x = F −1 (u) ⇔ F (x) = u ⇔ 1 − e−λx = u ln (1 − u) . ⇔ 1 − u = e−λx ⇔ x = − λ Como consecuencia, el algoritmo ser´ıa 1. Generar U ∼ U (0, 1). ) . 2. Devolver X = − ln(1−U λ que, en versi´ on simplificada, resulta: 0. Hacer L = −1/λ. 1. Generar U ∼ U (0, 1). 2. Devolver X = L · ln U 3. Repetir los pasos 1-2 tantas veces como se precise.

´ ´ 3.1. EL METODO DE INVERSION

3.1.1.

25

Ventajas e inconvenientes del m´ etodo de inversi´ on

La ventaja m´as importante del m´etodo de inversi´on es que, en general, es aplicable a cualquier distribuci´on continua. No obstante presenta algunos inconvenientes. Inconvenientes del m´ etodo de inversi´ on 1. En ocasiones la funci´on de distribuci´on no tiene una expresi´on expl´ıcita (por ejemplo para la distribuci´on normal). 2. A veces, a´ un teniendo una expresi´on expl´ıcita para F (x), es imposible despejar x en la ecuaci´on F (x) = u (es decir, encontrar una expresi´on expl´ıcita para F −1 (u)). 3. A´ un siendo posible encontrar x = F −1 (u), puede ocurrir que esta expresi´on sea complicada y conlleve una gran lentitud de c´alculo. El primero de los inconvenientes expuesto puede, a veces, subsanarse mediante el uso de aproximaciones de la distribuci´on en cuesti´on o mediante tabulaciones de la misma. El segundo suele abordarse mediante la utilizaci´on de m´etodos num´ericos para la resoluci´on aproximada de la ecuaci´on F (x) = u. El mayor problema pr´actico que esto conlleva es la necesidad de resolver num´ericamente una ecuaci´on cada vez que se desee generar un nuevo n´ umero aleatorio que siga esa distribuci´on (sin que los c´alculos hechos para el anterior valor simulado sean de ayuda).

3.1.2.

Algunas distribuciones que pueden simularse por el m´ etodo de inversi´ on

En lo que sigue u representa un n´ umero pseudoaleatorio en el intervalo (0,1).

Distribuci´ on exponencial (Exp(λ), λ > 0) Funci´ on de densidad: f (x) = λe−λx , x ≥ 0. Funci´ on de distribuci´ on: F (x) = 1 − e−λx , x ≥ 0. ln (1 − u) , Inversa de la distribuci´ on: F −1 (u) = − λ ln u Forma simplificada de la inversa: S(u) = − ,. λ Distribuci´ on de Cauchy 1 , x ∈ R. π (1 + x2 ) 1 arctan x Funci´ on de distribuci´ on: F (x) = + , x ∈ R. 2 π

Inversa de la distribuci´ on: F −1 (u) = tan π u − 12 . Forma simplificada de la inversa: S(u) = tan πu. Funci´ on de densidad: f (x) =

´ 26 CAP´ITULO 3. METODOS UNIVERSALES PARA VARIABLES CONTINUAS Distribuci´ on triangular en (0. a) 2 x 1− , 0 < x < a. a a   x2 2 x− , 0 < x < a. Funci´ on de distribuci´ on: F (x) = a 2a √
Inversa de la distribuci´ on: F −1 (u) = a 1 − 1 −√u . Forma simplificada de la inversa: S(u) = a (1 − u) . Funci´ on de densidad: f (x) =

Distribuci´ on de Pareto (a > 0 , b > 0) aba , x ≥ b. xa+1  a b , x ≥ b. Funci´ on de distribuci´ on: F (x) = 1 − x b . Inversa de la distribuci´ on: F −1 (u) = (1 − u)1/a b Forma simplificada de la inversa: S(u) = 1/a . u Funci´ on de densidad: f (x) =

Distribuci´ on de Weibull (λ > 0 , α > 0) α

Funci´ on de densidad: f (x) = αλα xα−1 e−(λx) , x ≥ 0 α Funci´ on de distribuci´ on: F (x) = 1 − e−(λx) , x ≥ 0. (− ln (1 − u))1/α −1 Inversa de la distribuci´ on: F (u) = . λ 1/α (− ln u) Forma simplificada de la inversa: S(u) = . λ La distribuci´ on de Laplace o doble exponencial (Lap(µ, λ) µ, λ > 0) Funci´ on de densidad: f (x) = λ2 e−λ|x−µ| , para todo x ∈ R. Funci´ on de distribuci´ on:  1 λ(x−µ) e si x < µ 2 F (x) = 1 −λ(x−µ) 1 − 2e si x ≥ µ Inversa de la distribuci´ on:  ln 2u   µ+ si u < 1/2 −1 λ F (u) =   µ − ln 2(1 − u) si u ≥ 1/2 λ pudi´endose generar por el m´etodo de inversi´on, usando como auxiliar T ∼ Exp(λ) : Algoritmo basado en el m´ etodo de inversi´ on 1. Generar U, V ∼ U (0. 1). ln U 2. Hacer T = . λ 3. Si V < 1/2, devolver X = µ + T. En caso contrario, hacer X = µ − T .

´ ´ 3.1. EL METODO DE INVERSION

3.1.3.

27

Inversi´ on aproximada

Como se coment´o anteriormente, en casos en los que no es posible determinar una expresi´on expl´ıcita para F (x) o en los que no se puede hallar la de su inversa, puede optarse por encontrar expresiones sencillas que aproximen razonablemente bien la funci´on F −1 (u). Ejemplo 3.1.2 A continuaci´on se detalla la aproximaci´ on encontrada por Odeh y Evans para la distribuci´on normal est´ andar. Estos autores consideran la funci´ on auxiliar √
√ A −2 ln v √
, −2 ln v g (v) = B −2 ln v siendo A (x) =

4 X

i

ai x y B (x) =

i=0

4 X

bi xi , con

i=0

a0 = −0. 322232431088 a2 = −0. 342242088547 a4 = −0. 0000453642210148 b1 = 0. 588581570495 b3 = 0. 103537752850

a1 = −1 a3 = −0. 0204231210245 b0 = 0. 0993484626060 b2 = 0. 531103462366 b4 = 0. 0038560700634

La aproximaci´on consiste en utilizar g (1 − u) en lugar de F −1 (u) para los valores de / [10−20 , 1 − 10−20 ] (que s´ olo ocurre u ∈ [10−20 , 12 ] y −g (u) si u ∈ [ 12 , 1 − 10−20 ]. Para u ∈ −20 con una probabilidad de 2 · 10 ) la aproximaci´ on no es recomendable. Algoritmo de Odeh y Evans 1. 2. 3. 4.

Generar U ∼ U (0, 1) . Si U < 10−20 o U > 1 − 10−20 entonces volver a 1. Si U < 0. 5 entonces hacer X = g (1 − U ) sino hacer X = −g (U ) . Devolver X.

´ 28 CAP´ITULO 3. METODOS UNIVERSALES PARA VARIABLES CONTINUAS

3.2.

El m´ etodo de aceptaci´ on/rechazo

Es un m´etodo universal alternativo al de inversi´on que est´a adaptado al caso en que, aunque se desconozca una f´ormula expl´ıcita para F (x) o sea dif´ıcil de resolver F (x) = u, s´ı se disponga de una expresi´on (preferiblemente sencilla) para la funci´on de densidad f (x). El m´etodo est´a basado en el siguiente resultado te´orico. Teorema 3.2.1 (de aceptacion/rechazo) Sea X una variable aleatoria con funci´ on de densidad f y sea U otra variable aleatoria, independiente de la anterior, con distribuci´ on U (0, 1). Entonces, para cada c > 0, la variable aleatoria bidimensional (X, c · U · f (X)) tiene distribuci´on uniforme en el recinto A = {(x, y) ∈ R2 /0 ≤ y ≤ cf (x)} Rec´ıprocamente, si dada una funci´ on de densidad f , un vector aleatorio (X, Y ) tiene distribuci´on uniforme sobre el conjunto A, entonces, su primera componente, X, es una variable aleatoria unidimensional con funci´ on de densidad f . El teorema anterior establece la equivalencia entre la simulaci´on de densidades unidimensionales y la simulaci´on de variables bidimensionales con distribuci´on uniforme sobre el hipografo de c · f (x) (el conjunto de puntos del plano que quedan por debajo de la gr´afica de c · f pero por encima del eje OX). La idea del algoritmo consistir´a en utilizar el rec´ıproco en el teorema para simular valores de ese tipo de distribuciones bidimensionales y luego tomar la primera componente. Para simular valores de esa distribuci´on bidimensional se usa tambi´en el teorema en sentido directo aplic´andolo a otra densidad auxiliar g, f´acil de simular. Sup´ongase que se desea simular una distribuci´on con densidad f y que no es factible hacerlo por el m´etodo de inversi´on. Consid´erese otra distribuci´on, con densidad g, f´acil de simular, de forma que exista cierta constante c > 0 tal que f (x) ≤ c · g (x) , para todo x ∈ R Teniendo en cuenta esta condici´on, Af ⊂ Acg , donde Af = {(x, y) /0 ≤ y ≤ f (x)} , Acg = {(x, y) /0 ≤ y ≤ c · g (x)} . son los hipografos de f y de c · g. Dado que la densidad g es f´acil de simular, puede aplicarse la primera parte del teorema de aceptaci´on/rechazo para encontrar una variable aleatoria bidimensional, (T, Y ), con distribuci´on uniforme sobre Acg . Aceptando tan s´olo los valores de (T, Y ) que pertenezcan a Af se tendr´a una variable bidimensional con distribuci´on uniforme sobre Af . T´ecnicamente hablando estamos afirmando que la distribuci´on condicionada (T, Y ) |(T,Y )∈Af es uniforme sobre Af . Finalmente la segunda parte del teorema permite obtener una variable con densidad f sin m´as que tomar la primera componente del par obtenido.

´ ´ 3.2. EL METODO DE ACEPTACION/RECHAZO

29

De forma m´as detallada, el m´etodo constar´ıa de los siguientes pasos: 1. Generar un valor T con densidad g. 2. Utilizar el teorema de aceptaci´on/rechazo para generar un par (T, Y ) con distribuci´on uniforme en Acg . 3. Comprobar si (T, Y ) ∈ Af . En caso afirmativo, hacer X = T , en caso contrario, volver al paso 1. El par de valores (T, Y ) se obtiene simplemente simulando U ∼ U (0, 1) y definiendo Y = c · U · g (T ). Adem´as, la condici´on (T, Y ) ∈ Af que hay que comprobar en el u ´ltimo paso equivale a Y ≤ f (T ). Teniendo todo esto en cuenta el algoritmo proceder´ıa como sigue: Algoritmo de aceptaci´ on/rechazo 1. Repetir 1.1. Generar U ∼ U (0, 1) y T con densidad g. 2. Hasta que c · U · g (T ) ≤ f (T ) . 3. Devolver X = T. Ejemplo 3.2.1 (densidades acotadas en un intervalo cerrado) Sea f una funci´ on de densidad con soporte en un intervalo cerrado [a, b] (i.e., {x/f (x) 6= 0} = [a, b]) de tal forma que ∃M > 0 tal que f (x) ≤ M ∀x (es decir, f es acotada superiormente). En este caso puede tomarse como densidad auxiliar g la de una U [a, b]. En efecto, tomando c = M (b − a) y teniendo en cuenta que  1 si x ∈ [a, b] b−a g (x) = 0 si x ∈ / [a, b] c se tiene que f (x) ≤ M = b−a = c · g (x), ∀x ∈ [a, b]. As´ı pues, el algoritmo quedar´ıa como sigue: 1. Repetir 1.1. Generar U, V ∼ U (0, 1). 1.2. Hacer T = a + (b − a) V . 2. Hasta que M · U ≤ f (T ). 3. Devolver X = T .

´ 30 CAP´ITULO 3. METODOS UNIVERSALES PARA VARIABLES CONTINUAS

3.2.1.

Eficiencia del algoritmo de aceptaci´ on/rechazo

Dado que el algoritmo de aceptaci´on/rechazo repite los pasos 1-2 un n´ umero aleatorio de veces, ser´a importante medir, de alguna forma, la eficiencia del mismo. En primer lugar, existen restricciones obvias para la constante c que ha de elegirse en el algoritmo. As´ı, debido al hecho de que f (x) ≤ c · g (x), se tiene Z Z 1 = f (x) dx ≤ c g (x) dx = c, luego c ≥ 1. Puede demostrarse adem´as que si c = 1 entonces f y g ser´ıan densidades correspondientes a la misma distribuci´on (iguales salvo en un conjunto de probabilidad cero) y, por tanto, si g es f´acil de simular igualmente f´acil lo ser´ıa f . As´ı pues, se tiene c > 1. La comprobaci´on que aparece en el paso 2 del algoritmo es c · U · g (T ) ≤ f (T ). La probabilidad de aceptaci´on de esta condici´on es R f (x) dx 1 area (Af ) =R = . p= area (Acg ) c c · g (x) dx De ´esta se obtiene la probabilidad de rechazo: q = c−1 . El flujo del algoritmo es c aleatorio y el n´ umero de repeticiones de los pasos 1-2 hasta poder generar un valor de f (paso 3) es una variable aleatoria, N , con distribuci´on geom´etrica (entendida ´esta como el n´ umero de pruebas necesarias hasta obtener el primer ´exito). En tales circunstancias el n´ umero medio de repeticiones de los pasos 1-2 es E (N ) =

1 =c p

luego c puede interpretarse como el n´ umero medio de comparaciones necesarias (o de repeticiones de los pasos 1-2, o de pares de variables (T, U ) que se necesitan generar) hasta obtener un valor simulado de la variable X. Es obvio, por tanto, que cuanto m´as cercano a 1 sea el valor de c m´as eficiente ser´a el algoritmo.

3.2.2.

Elecci´ on de c

Una vez fijada la densidad g es obvio que el mejor valor de c (que denotaremos por copt ) se obtiene al encontrar el m´as peque˜ no n´ umero real c que verifica f (x) ≤ c · g (x), es decir c≥

f (x) , para todo x del soporte de g (que ha de contener al de f ). g (x)

De esta forma, ha de cumplirse que f (x) 6= 0 ⇒ g (x) 6= 0 y adem´as f (x) . x/g(x)>0 g (x)

c ≥ m´ax

As´ı pues, el menor valor posible que cumple esta condici´on es f (x) . x/g(x)>0 g (x)

copt = m´ax

´ ´ 3.2. EL METODO DE ACEPTACION/RECHAZO

31

Ejemplo 3.2.2 (Simulacion de la normal mediante la doble exponencial) Se trata de simular la distribuci´on normal est´ andar, cuya funci´ on de densidad viene dada por x2 1 f (x) = √ e− 2 , para todo x ∈ R, 2π mediante aceptaci´on/rechazo, utilizando como densidad auxiliar la doble exponencial de par´ ametro 1 (o distribuci´on de Laplace, Lap(0,1) o Lap(1)), cuya funci´ on de densidad viene dada por 1 g (x) = e−|x| , para todo x ∈ R. 2 El valor ´optimo para c es 2

f (x) copt = m´ax = m´ax x∈R g (x) x∈R

x √1 e− 2 2π 1 −|x| e 2

r =

2 m´ax eϕ(x) = π x∈R

r

2 m´axx∈R ϕ(x) e , π

2

donde ϕ (x) = − x2 + |x|. Dado que esta funci´ on es sim´etrica, continua en toda la recta real y diferenciable tantas veces como se desee salvo en x = 0, bastar´ a encontrar su m´ aximo absoluto en el intervalo [0. ∞]: x > 0 ⇒ ϕ0 (x) = −x + 1, ϕ00 (x) = −1; {x > 0. ϕ0 (x) = 0} ⇔ x = 1 ϕ00 (1) < 0. De esta forma, ϕ alcanza un m´ aximo relativo en x = 1 y otro de id´entico valor en x = −1. Resulta f´acil demostrar que ambos son m´ aximos absolutos (por los intervalos de crecimiento y decrecimiento de la funci´ on). Consiguientemente, r r r 2 ϕ(1) 2 1/2 2e e = e = ' 1. 3155. copt = π π π Como consecuencia el algoritmo proceder´ıa del siguiente modo: 1. Repetir 1.1. Generar U, V, W ∼ U (0, 1). 1.2. Si W < 0. 5 hacer  2 T = ln V , sino hacer T = − ln V . 2. Hasta que U · exp T2 − |T | + 21 ≤ 1. 3. Devolver X = T . La condici´on que hay que comprobar para decidir si hay aceptaci´ on o rechazo surge de que r r  2   2  π T T 1 g (T ) 2e exp − |T | = U · exp − |T | + . c·U · = U 2 2 2 2 f (T ) π Dado que el n´ umero medio de repeticiones de los pasos 1-2 hasta que se obtiene un valor simulado on en el paso 2 es p π para X es c ' 1. 3155 y la probabilidad de aceptaci´ p = 1/c = 2e = 0. 76017, puede decirse que el algoritmo es bastante eficiente.

´ 32 CAP´ITULO 3. METODOS UNIVERSALES PARA VARIABLES CONTINUAS

3.2.3.

Elecci´ on de la densidad auxiliar g

Como se ha comentado anteriormente, un aspecto importante que influye en la eficiencia del m´etodo de aceptaci´on/rechazo es el valor de la constante c. Conocida la densidad auxiliar g sabemos c´omo elegir c de forma que el algoritmo sea lo m´as eficiente posible, sin embargo es obvio que algunas densidades auxiliares ser´ıan mejores candidatas que otras para conseguir un m´etodo eficiente. En general, cuanto m´as parecida sea la forma de g a la de f , m´as peque˜ no es el m´ınimo c necesario para conseguir que la gr´afica de c · g quede por encima de la de f . De todas formas, el problema de encontrar la densidad auxiliar g que ofrezca un c (´optimo) lo menor posible, no tiene soluci´on. Mejor dicho, tiene la soluci´on trivial g = f , que es absolutamente in´ util para la implementaci´on del algoritmo, pues si f era dif´ıcil de simular, no podemos tomar como g la propia f (ya que ser´ıa igual de dif´ıcil de simular). Una soluci´on intermedia al problema de elegir una funci´on de densidad auxiliar, g, adecuada consiste en tomar cierta familia param´etrica de densidades que presenten un abanico de formas entre las que haya alguna que se parece bastante a la de f : {gθ /θ ∈ Θ}, encontrar el valor de c ´optimo para cada densidad de esa familia: cθ = m´ax x

f (x) gθ (x)

y, finalmente, elegir el mejor valor del par´ametro, θ0 , en el sentido de ofrecer el menor posible cθ : f (x) . cθ0 = m´ın m´ax x gθ (x) θ∈Θ

´ ´ 3.2. EL METODO DE ACEPTACION/RECHAZO

33

Ejemplo 3.2.3 Sup´ongase que se desea utilizar como densidad auxiliar en el m´etodo de aceptaci´on/rechazo, para simular la normal est´ andar, la doble exponencial de par´ ametro λ > 0 (Lap(0,λ) o Lap(λ)) , gλ (x) =

λ −λ|x| e , para todo x ∈ R. 2

Si pretendemos encontrar el mejor valor de λ, en t´erminos de eficiencia del algoritmo, debemos calcular 2

f (x) = gλ (x)

cλ0 = m´ın m´ax λ>0 x∈R

x √1 e− 2 2π m´ın m´ax λ −λ|x| λ>0 x∈R e 2

.

De forma totalmente an´aloga a la vista para el caso λ = 1, se tiene 2

cλ =

x √1 e− 2 2π m´ax λ −λ|x| x∈R e 2

1 = λ

r

r

2 1 m´ax eϕλ (x) = π x∈R λ

2 m´axx∈R ϕλ (x) e , π

2

donde ϕλ (x) = − x2 + λ |x|. De forma totalmente similar a aquel caso puede probarse que ϕλ alcanza su m´aximo absoluto en los puntos x = ±λ, siendo dicho valor m´ aximo λ2 ϕλ (±λ) = 2 . Como consecuencia, 1 cλ = λ

r

λ2

2 ϕλ (±λ) e 2 e = π λ

r

2 . π

Ahora debemos encontrar λ0 tal que cλ0 = m´ınλ>0 cλ : dcλ = dλ

r

λ2

2 λe 2 λ − e π λ2

λ2 2

r =

2e π

λ2 2

(λ2 − 1) , λ2

r λ2 λ2 λ2 d2 cλ 2 [λe 2 (λ2 − 1) + e 2 2λ]λ2 − e 2 (λ2 − 1) 2λ = dλ2 π λ4 r λ2 r λ2 2 e 2 (λ5 + λ3 − 2λ3 + 2λ) 2 e 2 (λ5 − λ3 + 2λ) = , = π λ4 π λ4 dcλ = 0 ⇔ λ = 1, ya que λ > 0 dλ r d2 cλ 2e =2 > 0, luego λ = 1 es un punto de m´ınimo. 2 dλ λ=1 π De esto se deduce que la mejor doble exponencial, como densidad auxiliar en el algoritmo, es la correspondiente a λ = 1, la usada en el ejemplo anterior.

´ 34 CAP´ITULO 3. METODOS UNIVERSALES PARA VARIABLES CONTINUAS

3.2.4.

El m´ etodo de aceptaci´ on/rechazo “squeeze”

Esta variante del m´etodo de aceptaci´on/rechazo es de gran utilidad en aquellas situaciones donde, para llevar a cabo la comprobaci´on de la condici´on c·U ·g (T ) ≤ f (T ) , tiene un elavado coste computacional evaluar f (T ). La idea del m´etodo consiste en encontrar dos funciones h1 y h2 , f´aciles de evaluar, que “aprieten”a f (i.e. h1 (x) ≤ f (x) ≤ h2 (x), ∀x), de manera que reduzcamos considerablemente el n´ umero de evaluaciones de ´esta. As´ı se conseguir´ıa un algoritmo m´as eficiente, pues reemplazar´ıamos las evaluaciones de f por evaluaciones de h1 o de h2 (mucho menos costosas computacionalmente). Algoritmo de aceptaci´ on/rechazo “squeeze” 1. Generar U ∼ U (0, 1) y T con densidad g. 2. 2.1 Si c · U · g (T ) ≤ h1 (T ) , hacer X = T. (Aceptaci´ on r´ apida) 2.2 Si c · U · g (T ) > h2 (T ) , volver a 1. (Rechazo r´ apido) 3. Si no se verifica ninguna de las condiciones del paso 2, comprobar la condici´ on c · U · g (T ) ≤ f (T ) . Si c · U · g (T ) ≤ f (T ) , hacer X = T. Si no, volver a 1. Ejemplo 3.2.4 (Simulaci´ on “squeeze”de la distribuci´ on normal) Utilizando el teorema de Taylor, puede comprobarse que la funci´ onn de densidad de la normal est´ andar x2 1 f (x) = √ e− 2 , para todo x ∈ R, 2π 2

2

1 − x2 + 1− x √ puede “apretarse”por h1 (x) = √ 2 y h2 (x) = 2π 2π

x4 8

.

Cap´ıtulo 4 M´ etodos universales para la simulaci´ on de variables discretas En lo que sigue se expondr´an algunos m´etodos generales para simular distribuciones discretas. En concreto, se estudiar´a el m´etodo de la transformaci´on cuantil en su versi´on cl´asica y con etiquetados ´optimos, el m´etodo de las tablas gu´ıa y los m´etodos de truncamiento. El problema consiste en simular una variable aleatoria discreta, X, que toma los valores x1 , x2 , . . ., xn (. . .), con probabilidades pj = P (X = xj ), j = 1, 2, . . . , n (. . .). Un planteamiento est´andar, equivalente al anterior, consiste en resolver la cuesti´on de simular la variable aleatoria I que toma los valores 1, 2, . . . , n (. . .) con las mismas probabilidades pj , j = 1, 2, . . . , n (. . .).

4.1.

El m´ etodo de la transformaci´ on cuantil

Este m´etodo es una adaptaci´on del m´etodo de inversi´on (v´alido para el caso continuo) a distribuciones discretas. En primer lugar veamos porqu´e el m´etodo de inversi´on no es aplicable directamente en este caso. Dada una variable aletoria discreta, su funci´on de distribuci´on viene dada por X F (x) = pj , ∀x ∈ R. xj ≤x

Supondremos (por comodidad) que los valores que toma la variable ya est´an ordenados y nos ce˜ niremos al caso finito. De esta forma tendr´ıamos: x1 < x2 < · · · < xn . En este caso es obvio que el resultado dado por el teorema de inversi´on no es cierto ya que la variable aleatoria F (X) toma s´olo los valores p1 , p1 +p2 , . . ., p1 +p2 +· · ·+pn . Siendo, por tanto, discreta y no pudiendo tener distribuci´on U (0, 1).

35

´ CAP´ITULO 4. METODOS UNIVERSALES PARA VARIABLES DISCRETAS

36

De la misma forma, dada una variable U ∼ U (0, 1), tampoco puede ser cierto que F (U ) tenga la misma distribuci´on que X. De hecho F −1 no est´a definida de forma u ´nica pues las funciones de distribuci´on discretas no tienen inversa (para casi todo u ∈ [0, 1] no hay ning´ un x tal que F (x) = u y para un n´ umero finito (o infinito numerable) de u ∈ [0, 1] se tiene que existe todo un intervalo de valores para x cumpliendo F (x) = u). A pesar de ello puede definirse la llamada funci´on cuantil (o inversa generalidada) de una distribuci´on cualquiera F a partir de −1

Q (u) = ´ınf {x ∈ R/F (x) ≥ u} , ∀u ∈ (0, 1) . Es obvio que esta funci´on siempre est´a definida y que cuando F sea invertible, Q = F −1 . El siguiente teorema da un resultado que generaliza al teorema de inversi´on a situaciones en las que F no es invertible. Teorema 4.1.1 (de inversion generalizada) Sea X una variable aleatoria con funci´ on de distribuci´on F y con funci´on cuantil Q. Si U es una variable aleatoria con distribuci´ on U (0, 1), la variable Q (U ) tiene la misma distribuci´ on que X. Demostraci´ on: Sea G la funci´on de distribuci´on de Q (U ). Dado x ∈ R, se tiene G (x) = P (Q (U ) ≤ x) = P (´ınf {y ∈ R/F (y) ≥ U } ≤ x) Z F (x) = P (F (x) ≥ U ) = du = F (x) . 0

A partir del teorema de inversi´on generalizada puede obtenerse un algoritmo general para simular cualquier distribuci´on de probabilidad discreta. Es el llamado algoritmo de transformaci´on cuantil o de inversi´on generalizada. Algoritmo de transformaci´ on cuantil 1. Generar U ∼ U (0, 1). 2. Devolver X = Q (U ). La mayor dificultad en la implementaci´on del algoritmo radica en el c´alculo de ( , j ) X Q (U ) = ´ınf {x ∈ R/F (x) ≥ U } = ´ınf xj pi ≥ U i=1

= xk , tal que

k X i=1

pi ≥ U >

k−1 X i=1

pi .

´ ´ CUANTIL 4.1. EL METODO DE LA TRANSFORMACION

37

Todo el problema radica, por tanto, en encontrar el valor, k, de la variable, I, que guarda las etiquetas, para el cual la funci´on de distribuci´on supera o iguala por primera vez al valor de U . Este valor puede hallarse mediante una b´ usqueda secuencial, utilizando el siguiente algoritmo: Algoritmo de transformaci´ on cuantil con b´ usqueda secuencial 1. Generar U ∼ U (0, 1). 2. Hacer I = 1 y S = p1 . 3. Mientras U > S hacer 3.1. I = I + 1 y S = S + pI . 4. Devolver X = xI . Si se desea generar un gran n´ umero de valores de la variable X (que es P lo m´as habitual) puede resultar m´as eficiente calcular previamente las cantidades Sj = ji=1 pj de forma recursiva: S1 = p1 , Sj = Sj−1 + pj para j = 2, 3, . . . , n y hacer la comparaci´on U > SI en el paso 3 del algoritmo anterior. De esta forma se evita lo que podr´ıan ser c´alculos repetitivos de las mismas sumas de probabilidades al simular distintos valores de X. Ejemplo 4.1.1 (Simulacion de la distribucion de Poisson) T´ omese una variable, X, con distribuci´on de Poisson de par´ ametro λ, que toma los valores x1 = 0, x2 = 1, . . . con probabilidades pj = P (X = xj ) = P (X = j − 1) =

e−λ λj−1 , j = 1, 2, . . . (j − 1)!

El algoritmo de inversi´on con b´ usqueda secuencial viene dado por 1. Generar U ∼ U (0, 1). 2. Hacer I = 1 y S = e−λ . 3. Mientras U > S hacer −λ λI−1 3.1. I = I + 1 y S = S + e (I−1)! . 4. Devolver X = I − 1. Debido a que esta forma de etiquetar los valores de la variable conlleva el desfase de una unidad en los ´ındices, es recomendable ajustar el algoritmo para evitar este efecto.Tambi´en, para simplificar los c´ alculos que aparecen en el paso 3.1, es conveniente calcular las probabilidades de forma recursiva P (X = j) =

e−λ λj λ e−λ λj−1 λ = = P (X = j − 1) , j! j (j − 1)! j

As´ı, el algoritmo optimizado es 1. Generar U ∼ U (0, 1). 2. Hacer I = 0, p = e−λ y S = p. 3. Mientras U > S hacer 3.1. I = I + 1, p = λI p y S = S + p. 4. Devolver X = I.

38

´ CAP´ITULO 4. METODOS UNIVERSALES PARA VARIABLES DISCRETAS

4.1.1.

Eficiencia del algoritmo

Dada la forma del algoritmo general para simular una distribuci´on discreta mediante el m´etodo de la transformaci´on cuantil utilizando b´ usqueda secuencial, es f´acil probar que el n´ umero de comprobaciones de la forma U > S es precisamente igual a I, el valor de la variable que contiene las etiquetas. Como el valor de I es aleatorio y variar´a con cada ejecuci´on del algoritmo, una medida de la eficiencia del mismo ser´a el n´ umero medio de comparaciones del paso 3, es decir,  Pn si X toma un n´ umero finito (n) de valores j=1 jpj P∞ E (I) = si X toma un infinitos valores j=1 jpj Resulta pues evidente que, como no existe una u ´nica forma de etiquetar los valores que toma la variable en cuesti´on, habr´a quiz´a alg´ un etiquetado que ofrezca un menor n´ umero medio de comparaciones en el paso 3 del algoritmo que el etiquetado original (que obedece a la idea de ordenar de forma creciente los valores que toma la variable). Ejemplo 4.1.2 Consid´erese la variable aleatoria discreta X con distribuci´ on dada por P (X = 3) = 0. 1, P (X = 5) = 0. 3, P (X = 7) = 0. 6 Tomando x1 = 3, x2 = 5, x3 = 7, se tiene un etiquetado I con distribuci´ on P (I = 1) = 0. 1, P (I = 2) = 0. 3, P (I = 3) = 0. 6 y, por tanto, con media E (I) = 1 · 0. 1 + 2 · 0. 3 + 3 · 0. 6 = 2. 5. Si, por el contrario, consideramos el etiquetado x01 = 7, x02 = 5, x03 = 3, se tiene que P (I 0 = 1) = 0. 6, P (I 0 = 2) = 0. 3, P (I 0 = 3) = 0. 1 y as´ı E (I 0 ) = 1 · 0. 6 + 2 · 0. 3 + 3 · 0. 1 = 1. 5. Se observa que E (I 0 ) es sensiblemente inferior a E (I) y, por tanto, el segundo etiquetado proporciona un algoritmo m´ as eficiente que el dado por el etiquetado l. Como parece deducirse del ejemplo anterior, un etiquetado ser´a tanto mejor cuanto menores sean las etiquetas que se asignen a los valores que tienen mayor probabilidad. Dicho de otra forma, el etiquetado que se obtiene al ordenar los valores en orden decreciente de probabilidad. Cuando la variable a simular tiene un n´ umero finito de valores: x1 , x2 , . . ., xn , al implementar el m´etodo de la transformaci´ o n cuantil con b´ usqueda secuencial Pn−1 Pn directa, una vez comprobado que U > j=1 pj , no es necesario comprobar U > j=1 pj = 1 (que siempre es falso), sin´o que generamos xn sin necesidad de efectuar esa comparaci´on. Por ese motivo el n´ umero medio de comparaciones ser´ıa realmente: n−1 X j=1

jpj + (n − 1) pn .

´ ´ CUANTIL 4.1. EL METODO DE LA TRANSFORMACION

39

Ejemplo 4.1.3 Consideremos la variable aleatoria discreta con distibuci´ on P (X = 1) = 0. 11, P (X = 3) = 0. 3, P (X = 5) = 0. 25, P (X = 7) = 0. 21, P (X = 9) = 0. 13. Tomando el etiquetado x1 = 1, x2 = 3, x3 = 5, x4 = 7 y x5 = 9, el n´ umero medio de comparaciones del algoritmo es E (I) = 0. 11 · 1 + 0. 3 · 2 + 0. 25 · 3 + (0. 21 + 0. 13) · 4 = 2. 82 Mientras que, utilizando el etiquetado ´ optimo x1 = 3, x2 = 5, x3 = 7, x4 = 9 y x5 = 1, el n´ umero medio de comparaciones se reduce a E (I) = 0. 3 · 1 + 0. 25 · 2 + 0. 21 · 3 + (0. 13 + 0. 11) · 4 = 2. 39

4.1.2.

C´ alculo directo de la funci´ on cuantil

En ocasiones el m´etodo de la transformaci´on cuantil puede acelerarse computacionalmente porque, mediante c´alculos directos, es posible encontrar el valor de la funci´on cuantil en cualquier U , en un tiempo de computaci´on m´ınimo (evitando el bucle de b´ usqueda en el que se van acumulando las probabilidades). Ejemplo 4.1.4 (la distribuci´ on uniforme discreta en {1, 2, . . . , n}) En este caso la masa de probabilidad viene dada por pj =

1 , para j = 1, 2, . . . n. n

De esta forma se tiene k X i=1

pi ≥ U >

k−1 X i=1

pi ⇔

k k−1 ≥U > ⇔ k ≥ nU > k − 1. n n

Esta u ´ltima condici´on equivale a k = dnU e + 1, siendo dxe la parte entera de x. El algoritmo resulta: 1. Generar U ∼ U (0, 1). 2. Devolver X = dnU e + 1.

40

´ CAP´ITULO 4. METODOS UNIVERSALES PARA VARIABLES DISCRETAS

Ejemplo 4.1.5 (la distribuci´ on geom´ etrica) La distribuci´ on geom´etrica representa el n´ umero de fracasos antes del primer ´exito y tiene la siguiente masa de probabilidad P (X = j) = p (1 − p)j , j = 0. 1, . . . Para un valor j entero no negativo su funci´ on de distribuci´ on viene dada por F (j) =

j X

p (1 − p)i =

i=0

p (1 − p)j+1 − p = 1 − (1 − p)j+1 . 1−p−1

Como consecuencia se tiene F (k) ≥ U > F (k − 1) ⇔ 1 − (1 − p)k+1 ≥ U > 1 − (1 − p)k ⇔ (1 − p)k > 1 − U ≥ (1 − p)k+1 ⇔ k ln (1 − p) > ln (1 − U ) ≥ (k + 1) ln (1 − p) ln (1 − U ) ⇔k< ≤k+1 ln (1 − p) condici´ on que equivale a  ln (1 − U ) . k= ln (1 − p) 

El algoritmo proceder´ıa de la siguiente forma: 0. Hacer a = ln (1 − p). 1. Generar U ∼ U(0, 1).  2. Devolver X = lnaU .

´ 4.2. ALGORITMOS BASADOS EN ARBOLES BINARIOS

4.2.

41

Algoritmos basados en ´ arboles binarios

El uso de ´arboles binarios permite, en muchos casos, obtener algoritmos m´as eficientes que los basados en la b´ usqueda secuencial. Conceptos usados en este contexto: ´ Arbol: Grafo orientado, formado por un sistema de nodos conectados entre s´ı mediante una serie de arcos. Nodo ra´ız: Nodo del cual parten arcos pero al cual no llegan arcos. Nodo terminal: Nodo al cual llegan arcos pero del cual no parten arcos. Profundidad de un nodo: N´ umero de nodos que le preceden. ´ ´ Arbol binario: Arbol en el que todo nodo, a excepci´on de los nodos terminales, tiene dos nodos hijos.

Descripci´on de un ´arbol binario

ARCO 

 







aa  aa a 

NODO RA´ IZ

    

@ @ 













    

NODO TERMINAL



  c c c  





   

 l l l l 











42

´ CAP´ITULO 4. METODOS UNIVERSALES PARA VARIABLES DISCRETAS

Para la generaci´on de una variable aleatoria discreta, X, con funci´on de masa de probabilidad P (X = xi ) = pi , i = 1, 2, . . . , n se tratar´a de encontrar un ´arbol binario con n nodos terminales (uno para cada valor que se necesite generar), con profundidades di , i = 1, 2, . . . , n, de manera que n X

pi di

i=1

sea m´ınima. Es decir, se tratar´a de asignar mayor profundidad a los nodos correspondientes a valores de X con menor probabilidad.

4.2.1.

´ Arboles de Huffman

Un ´arbol de Huffman es un ´arbol binario en el que los nodos se generan siguiendo los siguientes pasos: 1. Agrupar los nodos con menor probabilidad en un solo nodo con probabilidad igual a la suma de ambos. 2. En el ´arbol resultante (con un nodo menos) proceder como en el paso anterior, repitiendo este proceso hasta finalizar con un ´arbol con solo dos nodos.

´ 4.3. EL METODO DE LA TABLA GU´IA

4.3.

43

El m´ etodo de la tabla gu´ıa

El mayor problema computacionalP del m´etodo de la on cuantil consiste Ptransformaci´ k k−1 en encontrar el ´ındice k que cumple i=1 pi ≥ U > i=1 pi . Como ya se ha visto en los dos u ´ltimos ejemplos existen distribuciones para las cuales este valor k se puede calcular directamente. El m´etodo de la tabla gu´ıa consiste en hacer uso de la rapidez de c´alculo de la funci´on cuantil para alguna de esas distribuciones (f´acilmente simulable mediante el m´etodo de inversi´on generalizada) para abreviar al aximo el n´ uP mero de Pm´ k k−1 comparaciones necesarias a la hora de comprobar la condici´on i=1 pi ≥ U > i=1 pi . Consid´erese una variable aleatoria discreta con masa de probabilidad dada por pj , j = 1, 2, . . . , n y def´ınanse las sumas acumulativas de estas probabilidades Pj (que no son otra cosa que los valores que toma la funci´on de distribuci´on), qj = i=1 pi , que, para evitar c´alculos innecesarios, deben calcularse de forma recursiva: q0 = 0, qj = qj−1 + pj , j = 1, 2, . . . , n. Dada la variable aleatoria I, asociada al etiquetado original (o a otro) la idea del m´etodo consiste en construir n subintervalos equiespaciados contenidos en [0, 1] de la , i ) para i = 1, 2, . . . , n y luego definir los valores de la tabla gu´ıa forma Ji = [ i−1 n n    i gi = m´ax j qj < , para i = 1, 2, . . . , n n es decir, para cada intervalo se considera el valor m´as alto del ´ındice entero tal que la suma acumulada de probabilidades hasta ´el es menor que el extremo superior de dicho intervalo. Ejemplo 4.3.1 Tomemos como ejemplo la distribuci´ on discreta dada por p1 = 0. 13, p2 = 0. 25, p3 = 0. 17, p4 = 0. 1, p5 = 0. 24 y p6 = 0. 11. Se tiene que q1 = 0. 13, q2 = 0. 38, q3 = 0. 55, q4 = 0. 65, q5 = 0. 89 y q6 = 1. Los valores de la tabla gu´ıa son g1 = 1, g2 = 1, g3 = 2, g4 = 4, g5 = 4, g6 = 5. A la hora de aplicar el m´etodo de la transformaci´on cuantil, dado el valor de U , es inmediato detectar en cu´al de los intervalos Ji ha ca´ıdo, basta con hacer i = dnU e + 1. Lo u ´nico que resta por hacer, una vez encontrado este ´ındice, es obtener el valor del ´ındice I a simular. Dicho valor ser´a gi + 1 si ya ocurre que U > qgi . En caso contrario deberemos encontrar el primer ´ındice j = gi − 1, gi − 2, . . . , 0, para el cual se cumple U > qj y luego hacer I = j + 1.

44

´ CAP´ITULO 4. METODOS UNIVERSALES PARA VARIABLES DISCRETAS

Algoritmo de simulaci´ on mediante una tabla gu´ıa 1. Generar U ∼ U (0, 1). 2. Hacer i = dnU e + 1. 3. Hacer j = gi . 4. Mientras U ≤ qj hacer j = j − 1. 5. Devolver I = j + 1. Por otra parte, los valores de la tabla gu´ıa pueden calcularse f´acilmente de forma r´apida seg´ un el siguiente algoritmo: Algoritmo de c´ alculo de la tabla gu´ıa 1. Desde i = 1 hasta n − 1 hacer gi = 0. 2. Hacer S = 0. 3. Desde i = 1 hasta n − 1 hacer 3.1. S = S + pi 3.2. j = dnSe + 1 3.3. gj = i 4. Desde i = 2 hasta n hacer gi = m´ax (gi−1 , gi ).

4.3.1.

Eficiencia del algoritmo

Cuando el valor U cae en el intervalo Ji , es obvio que el n´ umero medio de comparaciones en el paso 4 del algoritmo es menor o igual que 1 m´as el n´ umero de valores qj pertenecientes al intervalo Ji . Utilizando este hecho, la esperanza del n´ umero de comparaciones (N ) puede acotarse mediante n

n

1X 1X (1 + # {j /qj ∈ Ji }) = 1 + # {j /qj ∈ Ji } E (N ) ≤ n i=1 n i=1 = 1+

1 n−1 # {j /qj ∈ [0, 1)} = 1 + < 2. n n

En general, el m´etodo es aplicable para tablas gu´ıa de m elementos (donde m no tiene porqu´e ser necesariamente igual a n). En tal caso el intervalo [0, 1) se divide en m subintervalos, pudiendo acotar el n´ umero medio de comparaciones mediante E (N ) ≤ n umero exhorbitante de posibles 1+ m . Gracias a este argumento, para variables con un n´ valores, pueden utilizarse tablas gu´ıa de un n´ umero m´as moderado de elementos de forma que la tabla no ocupe demasiada memoria y que, a la vez, el n´ umero medio de comparaciones est´e acotado por un valor moderado. As´ı, por ejemplo, para una variable discreta con 1.000.000 de posibles valores podr´ıamos utilizar una tabla gu´ıa de s´olo 10.000 elementos (para que no ocupe demasiado en memoria) obteniendo que el n´ umero medio de comparaciones estar´ıa acotado por 101.

´ 4.4. METODO DE TRUNCAMIENTO

4.4.

45

M´ etodo de truncamiento

La idea general de este m´etodo consiste en hacer uso de una distribuci´on continua auxiliar cuya funci´on de distribuci´on se parezca (en cierto sentido que se precisar´a m´as adelante) a la funci´on de distribuci´on de la variable discreta que se desea simular. Sup´ongase, sin p´erdida de generalidad, que se desea simular la variable I, que toma los valores 1, 2, . . ., n, con probabilidades p1 , p2 , . . ., pn . En este caso, la funci´on de distribuci´on de I viene dada por X F (x) = pi . i≤x

Sup´ongase, adem´as, que tenemos otra variable aleatoria continua, con funci´on de distribuci´on G (x) y ciertos valores a0 = −∞ < a1 < a2 < · · · < an−1 < an = ∞, tales que F (i) − F (i− ) = pi = G (ai ) − G (ai−1 ), i = 1, 2, . . . , n. Esta u ´ltima condici´on viene a garantizar que la probabilidad de que la variable continua caiga en el intervalo [ai−1 , ai ) coincide con la probabilidad con la que la variable discreta original toma el valor i. Si la distribuci´on continua es f´acil de simular, simplemente deberemos generar valores de la misma y luego transformarlos en valores de la variable I. Algoritmo de simulaci´ on por truncamiento 1. Generar T con distribuci´ on G. 2. Encontrar el valor i tal que ai−1 ≤ T < ai . 3. Devolver I = i. El m´etodo se hace especialmente r´apido cuando el valor de i puede obtenerse de forma inmediata a partir del valor de T . Uno de los casos en los que esto es as´ı se da cuando G (0) = 0 y los valores ai = i, i = 0, 1, . . . , n (o, incluso, infinitos valores ai de esta forma). En este caso el algoritmo resulta: Algoritmo de simulaci´ on por truncamiento a la parte entera 1. Generar T con distribuci´ on G. 2. Hacer I = dT e + 1.

46

´ CAP´ITULO 4. METODOS UNIVERSALES PARA VARIABLES DISCRETAS

Ejemplo 4.4.1 (simulacion de la geom´ etrica por truncamiento) La masa de probabilidad de la distribuci´on geom´etrica es P (X = j) = P (I = j + 1) = p (1 − p)j , j = 0, 1, . . . Consid´erese como variable aleatoria continua auxiliar la exponencial, que tiene funci´ on de distribuci´ on dada por  1 − e−λx si x ≥ 0 G (x) = 0 si x < 0 Ahora bien,
G (i) − G (i − 1) = 1 − e−λi − 1 − e−λ(i−1) = e−λ(i−1) − e−λi
i−1

= p (1 − p)i−1 = e−λ(i−1) 1 − e−λ = 1 − e−λ e−λ siempre que tomemos p = 1 − e−λ . De esta forma se tiene G (i) − G (i − 1) = P (X = i − 1) = P (I = i) = pi y el algoritmo resultar´ıa: 0. 1. 2. 3.

Hacer λ = − ln (1 − p). Generar U ∼ U (0, 1). Hacer T = − lnλU . Devolver X = dT e.

Este es el mismo algoritmo que se obten´ıa anteriormente cuando razon´ abamos c´ omo calcular directamente el valor de la funci´ on cuantil para la distribuci´ on geom´etrica.

Cap´ıtulo 5 M´ etodos espec´ıficos para la simulaci´ on de distribuciones notables En este cap´ıtulo se estudiar´an algoritmos espec´ıficos para simular algunas de las distribuciones de probabilidad m´as importantes. La mayor´ıa de ellos son aplicaciones de los m´etodos generales ya expuestos, quiz´a con alguna particularidad.

5.1. 5.1.1.

Distribuciones continuas La distribuci´ on uniforme

Una vez generada la variable U con distribuci´on U (0, 1), la variable X con distribuci´on U (a, b) se obtiene haciendo X = a + (b − a)U. Algoritmo para generar la U (a, b) 1. Generar U ∼ U (0, 1). 2. Hacer X = a + (b − a)U.

47

´ DE DISTRIBUCIONES NOTABLES CAP´ITULO 5. SIMULACION

48

5.1.2.

La distribuci´ on normal

Se trata de simular X con distribuci´on normal est´andar, ya que la variable Y ∈ N (µ, σ) , con par´ametros arbitrarios (µ ∈ R, σ > 0), puede simularse mediante Y = µ + σX.

M´ etodo de Box-M¨ uller Se basa en la siguiente propiedad: √ Dadas E ∈ exp √ (1) y U ∈ U (0, 1) , variables aleatorias independientes, las variables 2E cos 2πU y 2Esen2πU son variables independientes con distribuci´on N (0, 1) . Algoritmo de Box-M¨ uller 1. Generar U, V√∼ U (0, 1). 2. Hacer W1 = −2 ln U y W2 = 2πV . 3. Devolver X1 = W1 cos W2 , X2 = W1 senW2 .

M´ etodo polar Se basa en una propiedad que da la distribuci´on condicionada a cierto suceso de un par de variables transformadas de otras uniformes: Dadas dos variables independientes V1 y V2 , con distribuci´on U (−1, 1), la distribuci´on condicionada s s ! 2 2 −2 ln (V1 + V2 ) −2 ln (V12 + V22 ) , V2 V1 2 2 V12 + V22 V12 + V22 V1 +V2 ≤1

 es N2

0 0

   1 0 , . 0 1

Algoritmo polar 1. Generar U1 , U2 ∼ U (0, 1). 2. Hacer V1 = 2U1 − 1, V2 = 2U2 − 1 y W = V12 + V22 . 3. Si W > 1 entonces volver a 1. q −2 ln W . 4. Hacer Y = W 5. Devolver X1 = V1 Y , X2 = V2 Y .

5.1. DISTRIBUCIONES CONTINUAS

49

M´ etodo del Teorema Central del L´ımite Como su propio nombre indica, este m´etodo se deriva a partir del Teorema Central del L´ımite: Dadas variables aleatorias T1 , T2 , . . ., Tn , independientes e id´enticamente distribuidas, con media µT y varianza σT2 finitas, se tiene que T − µT √ n 'N (0, 1) , σT si n es suficientemente grande. Este teorema puede aplicarse para simular una N (0, 1) tomando variables con otra distribuci´on m´as f´acil de simular. El caso m´as habitual es elegir Ti = Ui ∈ U (0, 1) y n = 12 (por simplicidad de c´alculo). De esta forma, la variable 12

X U − µU √ n= Ui − 6 σU i=1 tiene distribuci´on aproximadamente N (0, 1) . Algoritmo basado en el TCL 1. Generar U1 , U2 , . . . , U12 ∼ U (0, 1). 2. Devolver X = U1 + U2 + · · · + U12 − 6.

5.1.3.

La distribuci´ on de Cauchy

Esta distribuci´on puede definirse, de forma general, dependiendo de dos par´ametros: µ el de localizaci´on y σ > 0 el de escala. Su funci´on de densidad viene dada por f (x) =

π

σ2

σ
, para todo x ∈ R. + (x − µ)2

Un sencillo c´alculo permite hallar su funci´on de distribuci´on:   1 x−µ 1 F (x) = arctan + , π σ 2 pudi´endose implementar el m´etodo de inversi´on. Algoritmo basado en el m´ etodo de inversi´ on 1. Generar U ∼ U (0, 1). 2. Devolver X = σ tan (πU ) + µ.

´ DE DISTRIBUCIONES NOTABLES CAP´ITULO 5. SIMULACION

50

5.1.4.

La distribuci´ on exponencial

Se simula utilizando el m´etodo de inversi´on. Algoritmo basado en el m´ etodo de inversi´ on 1 1. Hacer L = − λ . 2. Generar U ∼ U (0, 1). 3. Devolver X = L · ln U .

5.1.5.

La distribuci´ on de Laplace o doble exponencial

Esta distribuci´on puede definirse, de forma general, dependiendo de dos par´ametros: µ el de localizaci´on y λ > 0 el de escala. Su funci´on de densidad viene dada por f (x) =

λ −λ|x−µ| e , para todo x ∈ R. 2

Su funci´on de distribuci´on es:  F (x) =

1

1 λ(x−µ) e 2 − 21 e−λ(x−µ)

si x < µ si x ≥ µ

pudi´endose generar por el m´etodo de inversi´on. Algoritmo basado en el m´ etodo de inversi´ on 1. Generar U, V ∼ U (0, 1). ln U . 2. Hacer T = λ 3. Si V < 1/2, devolver X = µ + T. En caso contrario, hacer X = µ − T .

5.1. DISTRIBUCIONES CONTINUAS

5.1.6.

51

Las distribuciones gamma y de Erlang

La distribuci´on gamma, Γ (a, p), depende de dos par´ametros: a > 0, par´ametro de escala, y p > 0, par´ametro de forma. La distribuci´on de Erlang no es m´as que la particularizaci´on de la gamma al caso en que p ∈ N. La funci´on de densidad de una Γ (a, p) viene dada por  ap p−1 −ax x e si x ≥ 0 Γ(p) f (x) = 0 si x < 0 Z ∞ xp−1 e−x dx es la llamada funci´on gamma de Euler. donde Γ (p) = 0

Puede demostrarse una relaci´on recursiva para Γ (p) sin m´as que hacer una integraci´on por partes: Tomando u = xp−1 y dv = e−x dx, se tiene, Z ∞ Z ∞  p−1

p−1 −x −x ∞ p−2 −x Γ (p) = x e dx = x −e − (p − 1) x −e dx 0 0 0 Z ∞ = (p − 1) xp−2 e−x dx = (p − 1) Γ (p − 1) 0

Esto permite reducir el c´alculo de Γ (p) al caso en que p ∈ (0, 1], ya que Γ (p) = (p − 1) (p − 2) · · · (p − [p]) Γ (p − [p]) Z Dado que Γ (1) =

∞

e−x dx = 1, la expresi´on anterior se simplifica cuando p ∈ N,

0

dando lugar a Γ (p) = (p − 1)! Cuando p = 1 la densidad de la gamma es  −ax ae si x ≥ 0 f (x) = 0 si x < 0 d es decir, Γ (a, 1) = exp (a).

52

´ DE DISTRIBUCIONES NOTABLES CAP´ITULO 5. SIMULACION

Una propiedad muy importante de la distribuci´on gamma es la llamada propiedad de reproductividad, que afirma que si se dispone de dos variables aleatorias independientes, X ∈ Γ (a, p1 ) e Y ∈ Γ (a, p2 ), la suma de ambas tambi´en es una gamma: X + Y ∈ Γ (a, p1 + p2 ). Este resultado se puede generalizar, por inducci´on, a la suma de cualquier n´ umero finito de variables gamma independientes con primer par´ametro, a, coincidente. En virtud de ello, si p es entero, dadas X1 , X2 , · · ·, Xp variables independientes Pp con distrinuci´on exp (a) (o, lo que es lo mismo, Γ (a, 1)) se tiene que su suma, i=1 Xi , tiene distribuci´on Γ (a, p). Como consecuencia, la distribuci´on de Erlang se puede simular f´acilmente como suma de exponeciales: Algoritmo reproductivo de simulaci´ on de la Erlang 1. Desde i = 1 hasta p hacer 1.1. Generar Ui ∼ U (0, 1). 1.2. Hacer Xi = −PlnaUi . 2. Devolver X = pi=1 Xi . Este algoritmo puede agilizarse computacionalmente definiendo previamente el valor L = − a1 y calculando un u ´nico logaritmo (en lugar de p) teniendo en cuenta que Pp ln Ui Pp Q 1 X = − = − ln ( pi=1 Ui ). As´ı se tiene: i i=1 a i=1 a Algoritmo reproductivo de simulaci´ on de la Erlang optimizado 1. Hacer L = − a1 . 2. Hacer S = 1. 3. Desde i = 1 hasta p hacer 3.1. Generar U ∼ U (0, 1). 3.2. Hacer S = S · U . 4. Devolver X = L · ln S. Los algoritmos anteriores s´olo son v´alidos para p entero, siendo adem´as muy lentos si p es grande. Por contra son muy simples y de f´acil implementaci´on. Como alternativa existen otros algoritmos m´as complicados que cubren tambi´en el caso en que p no sea entero.

5.1. DISTRIBUCIONES CONTINUAS

53

Los algoritmos que se describen a continuaci´on permiten generar la distribuci´on Γ(1, p). Si a 6= 1, la distribuci´on Γ(a, p) podr´a generarse a partir de la distribuci´on anterior, utilizando para ello la propiedad que afirma que si X ∈ Γ(1, p) entonces X/a tiene distribuci´on Γ(a, p). El algoritmo de Tadikamalla (1978), que s´olo es v´alido si p > 34 (a = 1), es un algoritmo de aceptaci´on/rechazo que usa como densidad auxiliar una doble exponencial centrada en p − 1 y con par´ametro de escala dado por λ=

1 2 √ = . θ 1 + 4p − 3

Para la implementaci´on del algoritmo debe definirse la funci´on   (θ − 1) x p−1 |x − (p − 1)| + (p − 1) (θ + 1) exp −x + T (x) = . θ (p − 1) θ

Algoritmo de Tadikamalla 1. Generar X, doble exponencial con media p−1 y par´ ametro de escala λ = 1θ = 1+√24p−3 . 2. Si X < 0 entonces volver a 1. 3. Generar U ∼ U (0, 1). 4. Si U ≤ T (X) entonces devolver X, sin´ o volver a 1.

54

´ DE DISTRIBUCIONES NOTABLES CAP´ITULO 5. SIMULACION

Como el anterior, el algoritmo de Cheng-Feast (1979) es un algoritmo de aceptaci´on/rechazo que es v´alido si p > 1 (a = 1). Algoritmo de Cheng-Feast 1. Generar U1 , U2 , independientemente, con distribuci´ on  1 p − 6p U1 U (0, 1) y hacer V = . (p − 1)U2 2(U2 − 1) 1 2. Si + V + ≤ 2 hacer X = (p − 1)V . p−1 V 2 log U2 En otro caso, si − log V + V ≤ 1 hacer X = (p − 1)V . p−1 3. Volver a 1. El algoritmo de Best-Ahrens-Dieter (1983) es tambi´en un algoritmo de aceptaci´on/rechazo que es v´alido si p < 1 (a = 1). Algoritmo de Best-Ahrens-Dieter √ pe−t . 0. Hacer t = 0,07 + 0,75 1 − p y b = 1 + t 1. Generar U1 , U2 , independientemente, con distribuci´ on U (0, 1) y hacer V = bU1 . 1 2. Si V ≤ 1 hacer W = tV p . En otro caso, ir a 3. 2−W 2.1. Si U2 ≤ , ir a 5. 2+W 2.2. Si U2 ≤ e−W , ir  a 5.  t(b − V ) 3. Hacer W = − log e Y = Wt . p 3.1. Si U2 (p + (1 − p)Y ) ≤ 1, ir a 5. 3.2. Si U2 ≤ Y p−1 , ir a 5. 4. Volver a 1. 5. Hacer X = W .

5.1. DISTRIBUCIONES CONTINUAS

5.1.7.

55

La distribuci´ on beta

Dadas dos variables aleatorias Y ∈ Γ (1, p) y Z ∈ Γ (1, q), independientes, se dice que la variable Y X= Y +Z tiene distribuci´on β (p, q), beta de par´ametros p y q. La funci´on de densidad de una β (p, q) viene dada por  p−1 q−1  x (1 − x) si x ∈ [0, 1] f (x) = β (p, q)  0 en otro caso Z siendo β (p, q) =

1

xp−1 (1 − x)q−1 dx.

0

Aunque existen multitud de algoritmos para simular la distribuci´on β (p, q) , probablemente, el m´as sencillo de todos es el que se obtiene, a partir de la distribuci´on gamma, como consecuencia de la propia definici´on. El algoritmo de Fox (1963) es adecuado para simular la distribuci´on beta cuando p, q ∈ N y son valores peque˜ nos. Algoritmo de Fox 1. Generar U1 , U2 , . . . , Up+q−1 ∼ U (0, 1). 2. Ordenarlos: U(1) ≤ U(2) ≤ · · · ≤ U(p+q−1) . 3. Devolver X = U(p) .

56

´ DE DISTRIBUCIONES NOTABLES CAP´ITULO 5. SIMULACION

Un m´etodo v´alido aunque p ´o q no sean enteros es el dado por el algoritmo de J¨ohnk (1964). Algoritmo de J¨ ohnk 1. Repetir. 1.1. Generar U, V ∼ U (0, 1). 1 1 1.2. Hacer Y = U p , Z = V q , S = Y + Z. 2. Hasta que S ≤ 1. 3. Hacer X = YS . Este m´etodo resulta extremadamente ineficiente para p ´o q mayores que 1. Esto es debido a que la condici´on S ≤ 1 del paso 2 puede tardar much´ısimo en verificarse. Por este motivo, el algoritmo de J¨ohnk s´olo es recomendable para p < 1 y q < 1. Como remedio a esto puede usarse el algoritmo de Cheng (1978) que es bastante m´as complicado de implementar pero tambi´en mucho m´as eficiente. Algoritmo de Cheng 1.1. Hacer α = p + q. 1.2. Si m´ın (p, q) ≤ 1 entonces hacer β = q α−2 . 2pq−α

1 , m´ın(p,q)

en otro caso hacer β =

1.3. Hacer γ = p + β1 . 2. Generar U1 , U2 ∼U (0,1). U1 3. Hacer V = β · ln 1−U y W = p · eV . 1   α 4. Si α · ln q+W + γV − ln 4 < ln (U12 U2 ) entonces volver a 1. 5. Devolver X =

W . q+W

5.1. DISTRIBUCIONES CONTINUAS

5.1.8.

57

La distribuci´ on chi-cuadrado de Pearson

Dadas variables aleatoriasP Z1 , Z2 , . . . , Zn independientes y con distribuci´on N (0, 1), diremos que la variable X = ni=1 Zi2 tiene distribuci´on chi-cuadrado con n grados de libertad (χ2n ). Su funci´on de densidad viene dada por n

x

x 2 −1 e− 2 f (x) = n n
, para todo x ≥ 0 22 Γ 2
d Es decir, χ2n = Γ 21 , n2 y los algoritmos vistos para la distribuci´on gamma son aplicables a este caso (para n grande es recomendable usar el algoritmo de Tadikamalla).
1 n d Adem´ a s, debido a la reproductividad de la distribuci´ o n gamma, se tiene que Γ , = 2
2
 
Γ 12 , n2 + Γ 12 , 12 , cuando n no sea par, siendo esta u ´ltima distribuci´on, Γ 21 , 12 , la del cuadrado de una normal est´andar. De esta forma, para n peque˜ no, se tiene el siguiente algoritmo para la simulaci´on de la chi-cuadrado: Algoritmo reproductivo para simular la chi-cuadrado n 1. Hacer m = 2 . 2. Hacer S = 1. 3. Desde i = 1 hasta m hacer 3.1. Generar U ∼ U (0, 1). 3.2. Hacer S = S · U . 4. Hacer X = −2 ln S. 5. Si n es impar hacer 6.1. Generar Z ∼ N (0, 1). 6.2. Hacer X = X + Z 2 . 7. Devolver X.

´ DE DISTRIBUCIONES NOTABLES CAP´ITULO 5. SIMULACION

58

5.1.9.

La distribuci´ on F de Fisher-Snedecor

Dadas dos variables aleatorias Y1 ∈ χ2m e Y2 ∈ χ2n independientes, la variable aleatoria Y1 /m X= Y2 /n tiene distribuci´on F de Fisher-Snedecor con m y n grados de libertad (Fm,n ). Su funci´on de densidad es
n n 2 f (x) =

β

m m
x 2 −1 m n , 2 2



n − m+n 2 , para todo x ≥ 0 x+ m

Adem´as de poder simularse a trav´es de algoritmos de simulaci´on de la chi-cuadrado, puede simularse mediante el uso de una distribuci´on beta. Algoritmo de simulaci´ on
de la F a trav´ es de la beta m n 1. Generar Z ∼ β 2 , 2 . nZ . 2. Hacer X = m(1−Z)

5.1.10.

La distribuci´ on t de Student

Dadas dos variables independientes Y1 ∈ N (0, 1) e Y2 ∈ χ2n , la variable aleatoria Y1

X=p

Y2 /n

tiene distribuci´on t de Student con n grados de libertad (tn ). Su funci´on de densidad es f (x) =

Γ Γ

n 2

n+1
2

√


 nπ

x2 1+ n

− n+1 2 , para todo x ∈ R.

La t de Student puede simularse f´acilmente teniendo en cuenta la relaci´on entre d esta distribuci´on y la F de Fisher-Snedecor: t2n = F1,n . Algoritmo de simulaci´ on de la t de Student a partir de la F 1. Generar U ∼ U (0, 1) y Z ∼ F1,n . √ √ 2. Si U < 0,5 entonces devolver X = Z, si no devolver X = − Z.

5.1. DISTRIBUCIONES CONTINUAS

5.1.11.

59

La distribuci´ on de Weibull

La distribuci´on de Weibull, W (λ, α), es una generalizaci´on de la distribuci´on exp (α). Su funci´on de densidad de probabilidad es α

f (x) = αλα xα−1 e−(λx) , para todo x ≥ 0 d En particular, W (λ, 1) = exp (λ). Puede simularse f´acilmente mediante el m´etodo de inversi´on (ligeramente optimizado). Algoritmo de inversi´ on para simular la distribuci´ on de Weibull 1. Generar U ∼ U (0, 1). 1 (− ln U ) α 2. Devolver X = . λ

5.1.12.

La distribuci´ on log´ıstica

Es la que tiene por funci´on de distribuci´on: F (x) =

1 1 + e−

x−a b

, ∀x ∈ R,

siendo a ∈ R y b > 0. Puede simularse f´acilmente mediante el m´etodo de inversi´on. Algoritmo para simular la distribuci´ on log´ıstica mediante inversi´ on 1. Generar U ∼ U (0, 1).
2. Devolver X = a − b ln U1 − 1 .

´ DE DISTRIBUCIONES NOTABLES CAP´ITULO 5. SIMULACION

60

5.1.13.

La distribuci´ on de Pareto

Tiene utilidad en ciencias como la Econom´ıa, donde en ocasiones sirve para modelizar distribuciones de rentas. Su densidad viene dada por ( f (x) =

aba xa+1 0

si x ≥ b si x < b

Como consecuencia, su funci´on de distribuci´on resulta  0
si x < b F (x) = a si x ≥ b 1 − xb y, por tanto, es simulable mediante inversi´on. Una versi´on optimizada del algoritmo es: Algoritmo de inversi´ on para simular la distribuci´ on de Pareto 1. Generar U ∼ U (0, 1). b 2. Devolver X = 1 . Ua

5.2. 5.2.1.

Distribuciones discretas La distribuci´ on uniforme discreta

Dado un conjunto finito de N elementos (que, sin p´erdida de generalidad, supondremos el conjunto {1, 2, . . . , N }) la distribuci´on uniforme discreta en dicho conjunto (o equiprobable sobre dicho conjunto) es la definida por P (X = i) = N1 , para i = 1, 2, . . . , N . Tanto el m´etodo de inversi´on (calculando expl´ıcitamente la funci´on cuantil) como el de truncamiento dan lugar al siguiente algoritmo. Algoritmo para simular la distribuci´ on uniforme discreta en {1, 2, . . . , N } 1. Generar U ∼ U (0, 1). 2. Devolver X = [N · U ] + 1.

5.2. DISTRIBUCIONES DISCRETAS

5.2.2.

61

La distribuci´ on binomial

La distribuci´on binomial de par´ametros n y p, B (n, p), se define como el n´ umero de ´exitos en n pruebas independientes, en las que la probabilidad de ´exito es p. Su masa de probabilidad es   n i P (X = i) = p (1 − p)n−i , para i = 0, 1, . . . , n. i

Puede simularse a partir de su definici´on: Algoritmo para la generaci´ on de la distribuci´ on B(n, p) 1. Hacer S = 0. 2. Repetir n veces 2.1. Generar U ∼ U (0, 1). 2.2. Si U ≤ p entonces hacer S = S + 1. 3. Devolver X = S. Este m´etodo es extremadamente lento cuando n es grande. Por eso, en ese caso, resulta m´as ventajoso utilizar el m´etodo de la tabla gu´ıa.

5.2.3.

La distribuci´ on de Poisson

Una variable aleatoria discreta, X, tiene distribuci´on de Poisson de par´ametro λ > 0 si su masa de probabilidad viene dada por e−λ λi , para i = 0, 1, . . . P (X = i) = i! La distribuci´on de Poisson puede simularse mediante el m´etodo de la transformaci´on cuantil con b´ usqueda secuencial.

´ DE DISTRIBUCIONES NOTABLES CAP´ITULO 5. SIMULACION

62

Tambi´en puede simularse haciendo uso de la relaci´on que guarda con la distribuci´on exponencial. As´ı, dadas variables aleatorias T1 , T2 , . . ., Tn , . . . independientes y con distribuci´on exp (λ), la variable aleatoria entera, X, que verifica X X

Ti ≤ 1 <

X+1 X

i=1

Ti

i=1

(definiendo X = 0 si T1 > 1) tiene distribuci´on Pois(λ). Las variables aleatorias Ti pueden simularse, utilizando valores Ui de una uniforme, mediante Ti = − lnλUi . En virtud de ello, se tiene X X

Ti ≤ 1 <

i=1

ln −

X Q

X+1 X i=1

 ln

Ui

i=1

λ X Y i=1

Ti ⇔ −

X X ln Ui i=1

X+1 Q i=1

≤ 1

λ X+1 Y

λ

≤1 ln

X+1 Y

! Ui

⇔

i=1

Ui .

i=1

As´ı, puede utilizarse el siguiente algoritmo: Algoritmo de simulaci´ on de la Poisson a traves de la exponencial 1. Hacer p = 1 y S = −1. 2. Repetir 2.1. Generar U ∼ U (0, 1). 2.2. Hacer p = p · U y S = S + 1. 3. Hasta que p < e−λ . 4. Hacer X = S. Tanto este algoritmo como el de la transformaci´on cuantil tienen el inconveniente de ser muy ineficientes cuando λ es grande. En ese caso, aunque la distribuci´on de Poisson tiene un n´ umero infinito de posibles resultados, es perfectamente aplicable el m´etodo de la tabla gu´ıa desembocando en una b´ usqueda secuencial cuando el intervalo elegido sea el u ´ltimo de la tabla. Esto mejora muy considerablemente la eficiencia del m´etodo.

5.2. DISTRIBUCIONES DISCRETAS

5.2.4.

63

La distribuci´ on geom´ etrica

Su masa de probabilidad es P (X = i) = p · (1 − p)i , para i = 0, 1, . . . Adem´as de poder simularse a partir de su definici´on (n´ umero de fracasos antes del primer ´exito), tambi´en puede hacerse por truncamiento. El algoritmo que resulta por este m´etodo es equivalente al basado en la expresi´on expl´ıcita de la funci´on cuantil. Algoritmo de truncamiento para la distribuci´ on geom´ etrica 1 . 0. Hacer L = − ln (1 − p) 1. Generar U ∼ U (0, 1). 2. Hacer T = L · ln U . 3. Devolver X = [T ].

5.2.5.

La distribuci´ on binomial negativa

La distribuci´on binomial negativa, BN (r, p) , generaliza a la geom´etrica, pudiendo interpretarse como el n´ umero de fracasos antes del r-´esimo ´exito. Su funci´on de masa de probabilidad es   i+r−1 r P (X = i) = p (1 − p)i , para i = 0, 1, . . . i

Debido a su reproductividad en el par´ametro r, puede simularse como suma de r variables geom´etricas, aunque este algoritmo puede ser muy costoso en tiempo de computaci´on si r es elevado. Existe tambi´en un m´etodo espec´ıfico basado en la propiedad   p , r ⇒ X ∈ BN (r, p) . X|Y ∈ Pois (Y ) , Y ∈ Γ 1−p Algoritmo condicional la binomial negativa  para simular  p 1. Simular L ∼ Γ ,r . 1−p 2. Simular X ∼Pois(L). 3. Devolver X.

64

´ DE DISTRIBUCIONES NOTABLES CAP´ITULO 5. SIMULACION

Cap´ıtulo 6 Simulaci´ on de distribuciones multidimensionales La simulaci´on de vectores aleatorios X = (X1 , X2 , . . . , Xd )0 que sigan cierto modelo de distribuci´on dado no es tarea siempre sencilla. En general, no resulta una extensi´on inmediata del caso unidimensional, aunque, si las variables que componen el vector son independientes, entonces bastar´a simular cada Xi con la distribuci´on marginal deseada (Fi ) y luego agrupar los valores simulados para cada componente en un vector. En la mayor parte de los casos de inter´es, las componentes del vector aleatorio son dependientes y el m´etodo anterior no es v´alido. A continuaci´on se ver´an algunos m´etodos generales para la simulaci´on de distribuciones multidimensionales.

6.1.

M´ etodo de las distribuciones condicionadas

Sup´ongase un vector aleatorio d-dimensional, con distribuci´on continua. Den´otese por f (x1 , x2 , . . . , xn ) su funci´on de densidad conjunta y consid´erese la primera densidad marginal, f1 (x1 ), y las sucesivas densidades condicionales f2 (x2 |x1 ), f3 (x3 |x1 , x2 ), . . ., fd (xd |x1 , x2 , . . . , xd−1 ). Gracias a la regla del producto, generalizada a funciones de densidad, se tiene f (x1 , x2 , . . . , xn ) = f1 (x1 ) · f2 (x2 |x1 ) · f3 (x3 |x1 , x2 ) · · · fd (xd |x1 , x2 , . . . , xd−1 ) y, como consecuencia, puede darse el siguiente algoritmo general: 1. Generar X1 con densidad f1 . 2. Desde i = 2 hasta d generar Xi con densidad fi (•|X1 , X2 , . . . , Xi−1 ). 3. Devolver X = (X1 , X2 , . . . , Xd )0 . Es inmediato comprobar que el m´etodo anteriormente expuesto es igualmente v´alido si las variables Xi son discretas o, incluso, si algunas son discretas y otras continuas. En tal caso se sustituir´ıa la densidad por la masa de probabilidad. As´ı pues, lo realmente importante para poder aplicar el m´etodo de las distribuciones condicionadas es conocer y saber simular la distribuci´on marginal de X1 y las distribuciones condicionadas del tipo Xi |X1 , X2 , . . . , Xi−1 para i = 2, 3, . . . , d. 65

´ DE DISTRIBUCIONES MULTIDIMENSIONALES 66CAP´ITULO 6. SIMULACION Ejemplo 6.1.1 (Algoritmo para simular la distribuci´ on normal bidimensional por el m´ etodo de las distribuciones condicionadas) Consideremos una    2  µ1 σ1 σ12 N2 , , µ2 σ12 σ22 por las propiedades de la distribuci´on normal, bastar´ a saber simular la distribuci´ on     2 0 σ1 σ12 N2 , σ12 σ22 0 

µ1 µ2

y luego sumarle el vector

 .

Dado que X1 ∈ N (0, σ1 ), se tiene que f1 (x1 ) =

σ1

1 √

x2 exp − 12 2σ1 2π 

Adem´ as





1

1 f (x1 , x2 ) = f (x) = p exp − x0 Σ−1 x 2 2π det (Σ) Como −1

Σ

1 = det (Σ)



σ22 −σ12 −σ12 σ12



1 = 2 2 2 σ1 σ2 − σ12





σ22 −σ12 −σ12 σ12



se tiene que 1 0 −1 σ 2 x2 − 2σ12 x1 x2 + σ12 x22 xΣ x= 2 1 2 2 2 (σ12 σ22 − σ12 ) y, por tanto, √   2 2  σ1 2π x21 σ2 x1 − 2σ12 x1 x2 + σ12 x22 f (x1 , x2 ) p = − 2 exp − f2 (x2 |x1 ) = 2 2 f1 (x1 ) 2 (σ12 σ22 − σ12 ) 2σ1 2π σ12 σ22 − σ12  2 2 2  2 4 2 2 ) x21 1 σ1 σ2 x1 − 2σ1 σ12 x1 x2 + σ1 x2 − (σ12 σ22 − σ12 = √ q 2 2 2 exp − 2 σ1 σ2 −σ12 2σ12 (σ12 σ22 − σ12 ) 2π 2 σ1   2 2 −2σ12 σ12 x1 x2 + σ14 x22 + σ12 x1 1 = √ q 2 2 2 exp − 2 σ1 σ2 −σ12 2σ12 (σ12 σ22 − σ12 ) 2π σ12   2 x2 σ12 2σ12 x1 x2 2 1 − + x + 2 1 σ12 σ14  = √ q 2 2 2 exp − 2 2 2 σ1 σ2 −σ12 σ1 σ2 −σ12 2 2π σ12 σ12   2  σ12 x1 x − 2 σ12 1   = √ q 2 2 2 exp −  2 2 2 σ1 σ2 −σ12 σ1 σ2 −σ12 2 σ2 2π 2 σ 1 1

que es la densidad de una N



σ 2 σ 2 −σ 2 σ12 x , 1 2σ2 12 σ12 1 1



.

´ 6.1. METODO DE LAS DISTRIBUCIONES CONDICIONADAS

67

En resumen, se tiene que si      2  X1 0 σ1 σ12 ∈ N2 , X2 0 σ12 σ22 q 2 2 2   σ1 σ2 −σ12 σ12 2 entonces X1 ∈ N (0, σ1 ) y X2 |X1 ∈ N σ2 X1 , . As´ı, el algoritmo de simuσ12 1 laci´ on consistir´ıa en los siguientes pasos: 1. 2. 3. 4. 5.

Simular Z1 , Z2 ∼ N (0, 1) independientes. Hacer Y1 = σ1 Z1 . q 2 2 2 σ1 σ2 −σ12 . Hacer Y2 = σσ122 Y1 + Z2 σ12 1 Hacer X1 = Y1 + µ1 , X2 = Y2 + µ2 . → − Devolver X = (X1 , X2 )t .

Ejemplo 6.1.2 (La distribuci´ on uniforme en el c´ırculo unitario). Se trata de la distribuci´on bidimensional continua cuya densidad es constante en dicho c´ırculo C = {(x1 , x2 )0 ∈ R2 /x21 + x22 ≤ 1}. Su funci´on de densidad viene dada por  f (x1 , x2 ) =

1 π

0

si si

(x1 , x2 )0 ∈ C (x1 , x2 ) 0 ∈ /C

La densidad marginal de la primera variable resulta p Z +√1−x21 2 1 − x21 1 dx2 = si x1 ∈ [−1, 1] f1 (x1 ) = √ π − 1−x21 π es decir,  f1 (x1 ) =

2 π

p

0

1 − x21

si si

x1 ∈ [−1, 1] x1 ∈ / [−1, 1]

Adem´as  q  q 1 f (x1 , x2 ) 1 π f2 (x2 |x1 ) = = √ 2 = p , si x2 ∈ − 1 − x21 , 1 − x21 2 1−x1 f1 (x1 ) 2 1 − x21 π

valiendo cero en otro caso. Se tiene entonces que  q  q 2 2 X2 |X1 ∈ U − 1 − X1 , 1 − X1 , siempre que X1 ∈ [−1, 1]. Finalmente, el algoritmo resulta: 1. Simular X1 con densidad f1 .h i p p 2. Simular X2 con densidad U − 1 − X12 , 1 − X12 . 3. Devolver X = (X1 , X2 )0 . Para el paso 1 puede utilizarse, por ejemplo, el m´etodo de aceptaci´ on/rechazo, pues se trata de una densidad acotada definida en un intervalo acotado.

´ DE DISTRIBUCIONES MULTIDIMENSIONALES 68CAP´ITULO 6. SIMULACION

6.2.

El m´ etodo de aceptaci´ on/rechazo

La idea general del m´etodo de aceptaci´on/rechazo es aplicable para simular variables aleatorias definidas en cualquier espacio (no s´olo en R). En particular puede usarse para simular vectores aleatorios de Rd . Sin embargo, en este contexto, resulta mucho m´as dif´ıcil encontrar una densidad auxiliar adecuada y, especialmente, conseguir que el n´ umero medio de comparaciones del m´etodo se mantenga dentro de unos l´ımites de eficiencia razonables cuando la dimensi´on es elevada. Ejemplo 6.2.1 (Simulaci´ on de puntos uniformemente distribu´ıdos sobre la “esfera” unitaria d-dimensional Cd )  Cd = (x1 , x2 , . . . , xd )0 ∈ Rd /x21 + x22 + · · · + x2d ≤ 1 Denotando por Vd (1), el “volumen” (la medida) de la esfera d-dimensional de radio 1 (en general, la de radio r verifica Vd (r) = rd Vd (1)), se tiene:  1 si (x1 , x2 , . . . , xd )0 ∈ Cd Vd (1) f (x1 , x2 , . . . , xd ) = 0 si (x1 , x2 , . . . , xd )0 ∈ / Cd Para simularvalores en Rd , con densidad f ,podemos utilizar como densidad auxiliar   d d la de una U [−1, 1] × [−1, 1] × · · · × [−1, 1] = U [−1, 1] , dada por  g (x1 , x2 , . . . , xd ) =

1 2d

0

si

xi ∈ [−1, 1], para todo i = 1, 2, . . . , d en otro caso

La constante c ´optima para la utilizaci´on del m´etodo de aceptaci´ on/rechazo es f (x) = c = m´ax → − x /g(x)>0 g (x)

1 Vd (1) 1 2d

2d = Vd (1)

y la condici´ on de aceptaci´on cU g (T) ≤ f (T) se convierte en 2d 1 1 U d 1[−1,1]d (T) ≤ 1C (T) − Vd (1) 2 Vd (1) d o, lo que es lo mismo, U 1[−1,1]d (T) ≤ 1Cd (T). Esta condici´ on equivale a que T ∈ Cd , es decir, que se verifique T12 + T22 + · · · + Td2 ≤ 1   Por otra parte, la simulaci´on de T ∼ U [−1, 1]d puede hacerse trivialmente mediante Ti ∼ U ([−1, 1]) para cada i = 1, 2, . . . , d, ya que las componentes son independientes. Como el valor de U es superfluo en este caso, el algoritmo queda: 1. 2. 3. 4.

Simular V1 , V2 , . . . , Vd ∼ U (0, 1) independientes. Para i = 1, 2, . . . , d hacer Ti = 2Vi − 1. Si T12 + T22 + · · · + Td2 > 1 entonces volver al paso 1. Devolver X = (T1 , T2 , . . . , Td )0 .

´ ´ O ETIQUETADO 6.3. METODOS DE CODIFICACION

69

Usando las f´ormulas del “volumen” de una “esfera” d-dimensional:  d/2 d π r   si d es par  (d/2)! Vd (r) = d d2 e+1 π d d2 e rd    2 si d es impar 1 · 3 · 5···d puede verse que el n´ umero medio de repeticiones de los pasos 1-3 del algoritmo, que d viene dado por la constante c = Vd2(1) , puede hacerse enormemente grande. As´ı, si d = 2 se tiene c =1.27, si d = 3 se tiene c =1.91, si d = 4 entonces c =3.24 y para d = 10 resulta c =401.5 que es un valor que hace que el algoritmo sea tremendamente lento en dimensi´on 10.

6.3.

M´ etodos de codificaci´ on o etiquetado

En el caso de que la funci´on de distribuci´on d-dimensional sea discreta existen m´etodos que permiten reducir la simulaci´on de dicha variable al contexto de simular una variable aleatoria discreta unidimensional. Estos m´etodos son conocidos como m´etodos de etiquetado o codificaci´on y la idea b´asica consiste en construir una funci´on h que codifique las posibles d-tuplas del conjunto donde toma valores la variable discreta, haciendo corresponder a cada uno un n´ umero entero no negativo diferente. Ejemplo 6.3.1 (Algoritmo para simular una variable bidimensional discreta (X1 , X2 )0 cada una de cuyas componentes toma valores enteros no negativos). El subconjunto de R2 en el que toma valores el vector aleatorio es Z+ × Z+ = Z+


2

= {(i, j) /i, j ∈ {0, 1, 2, . . .}}

Se tratar´a de definir una funci´ on biyectiva, h : Z+ × Z+ −→ Z+ , que permita etiquetar los pares de enteros. De esta forma, h induce sobre la variable transformada, C = h (X1 , X2 ), una masa de probabilidad (C)

pk


(X ,X ) := P (C = k) = P (h (X1 , X2 ) = k) = P (X1 , X2 ) = h−1 (k) =: ph−11 (k)2

Resulta inmediato, por tanto, obtener la masa de probabilidad de la variable discreta unidimensional C, a partir de la masa de probabilidad de la variable original (X1 , X2 ). De todas formas, debemos tener en cuenta que para que esto sea calculable en la pr´ actica en un tiempo razonable, la funci´ on h debe poder invertirse de forma r´ apida. As´ı pues para simular la variable (X1 , X2 ) podemos proceder mediante uno de los algoritmos posibles para simular C calculando en tantos pasos como sean necesarios los valores de la forma h−1 (k). Una posibilidad sencilla consiste en utilizar h (i, j) =

(i + j) (i + j + 1) +i 2

´ DE DISTRIBUCIONES MULTIDIMENSIONALES 70CAP´ITULO 6. SIMULACION Consideremos k ∈ Z+ , el valor (i, j) = h−1 (k) debe verificar (i + j) (i + j + 1) +i=k 2 Denotando ahora n = i + j, para encontrar (i, j) = h−1 (k) basta con hallar n e i, enteros positivos, con n ≥ i tales que h (i, j) = k ⇔

n (n + 1) +i=k 2 Debemos entonces encontrar el u ´nico n que cumple n (n + 1) n (n + 1) n (n + 1) + 2 (n + 1) (n + 1) (n + 2) ≤k≤ +n< = 2 2 2 2 2 2 Como adem´ as n < n (n + 1) y (n + 1) (n + 2) < (n + 2) , se tiene que ese valor n ha de verificar n2 < 2k < (n + 2)2 , es decir

m l√ m 2k − 2 < n ≤ 2k . l√ m l√ m Dicho de otro modo, se tiene que n ha de ser igual a 2k − 1 ´ o 2k . l√

n(n+1) Basta entonces con calcular la expresi´ l√ monl√2 m para esos posibles valoresl√de nm y comparar el resultado con 2k. As´ı, si 2k 2k + 1 > 2k entonces n = 2k − l√ m 2k . Finalmente se calcula 1 y, en caso contrario, n =

n (n + 1) y j =n−i 2 El c´ alculo de h−1 (k) es muy r´apido y el resto del algoritmo se reduce a la simulaci´ on de la variable unidimensional C. i=k−

−1 Ejemplo 6.3.2 Calculemos anterior √ por el procedimiento √  √  el valor h (16). Calculamos primeramente n = 2 · 16 = 2 · 16 = 32 = d5. 656 854 2e = 5. Luego calculamos 5 (5 + 1) = 30 ≤ 32 = 2 · 16, con lo cual n = 5. Adem´ as i = 16 − 5·6 =1y 2 −1 j = 5 − 1 = 4. As´ı pues se obtiene h (16) = (1, 4).

Aunque no entraremos con detalle en ello, conviene resaltar que es posible generad lizar este tipo de funciones de codificaci´on a (Z+ ) . Tambi´en es factible encontrar la inversa de tal funci´on generalizada (llamada funci´on de decodificaci´on) que se puede calcular eficientemente. Cuando la variable aleatoria X2 toma un n´ umero finito de valores (supongamos comprendidos entre 0 y M ), otra posible funci´on de codificaci´on, m´as sencilla es h (i, j) = (M + 1) i + j cuya inversa viene dada por −1

h

 (k) =

  k , kmod (M + 1) . M +1 d

Estas funciones de codificaci´on y decodificaci´on son generalizables a (Z+ ) y aplicables al caso en que el vector aleatorio X tome un n´ umero finito de valores.

´ ´ NORMAL MULTIVARIANTE71 6.4. METODOS PARA SIMULAR LA DISTRIBUCION

6.4.

M´ etodos para simular la distribuci´ on normal multivariante

Dado un vector µ = (µ1 , µ1 , . . . , µd )0 ∈ Rd y una matriz definida positiva   σ11 σ12 · · · σ1d  σ21 σ22 · · · σ1d    Σ=  .. .. . . ..   . . . .  σd1 σd2 · · · σdd la distribuci´on normal d-dimensional de par´ametros (µ , Σ) (que corresponden con su vector de medias y su matriz de varianzas-covarianzas), abreviadamente Nd (µ , Σ), es la que tiene densidad dada por   1 0 −d/2 −1/2 −1 f (x ) = (2π) det (Σ) exp − (x − µ ) Σ (x − µ) 2 Cuando Σ es diagonal    Σ= 

σ12 0 0 σ22 .. .. . . 0 0

··· ··· .. .

0 0 .. .

· · · σd2

   , 

se obtiene f´acilmente −d/2

f (x ) = (2π)

d Y

!−1/2 σi2

i=1





1 σ12

0

··· ··· ...

0 0 .. .



    1  0 σ122  0   (x − µ × exp  .. − 2 (x − µ )  ..  .   .  1 0 0 · · · σ2 d ! ! d d Y X (xi − µi )2 1 −d/2 = (2π) σi−1 exp − 2 i= σi2 i=1 !! d d Y Y 1 (xi − µi )2 √ exp − = = φµ1 ,σi (xi ) 2 2σ σ 2π i i i=1 i=1

   )  

siendo φµ1 ,σi la funci´on de densidad de una N (µi , σi ). De esta forma, cuando Σ es diagonal, las componentes son independientes y resulta trivial simular la Nd (µ , Σ) mediante el siguiente algoritmo: 1. Simular Z1 , Z2 , . . . , Zd ∼ N (0, 1) independientes. 2. Para i = 1, 2, . . . , d hacer Xi = µi + σi Zi . 3. Devolver X = (X1 , X2 , . . . , Xd )0 .

´ DE DISTRIBUCIONES MULTIDIMENSIONALES 72CAP´ITULO 6. SIMULACION Una propiedad que resulta muy u ´til para simular la distribuci´on Nd (µ , Σ) con Σ arbitraria es la siguiente. Proposici´ on 6.4.1 Si X ∈ Nd (µ , Σ) y A es una matriz de dimensi´ on p×d, de rango m´ aximo, con p ≤ d, entonces Y = AX ∈ Np (Aµ , AΣA0 ). Dada una variable aleatoria X ∈ Nd (µ , Σ), Y = X − µ ∈ Nd (0 , Σ). Si Σ es una matriz definida positiva, existe una matriz ortogonal H (es decir, tal que H−1 = H0 ) de forma que la matriz Λ = H0 ΣH es diagonal. De hecho, H es la matriz de cambio de base para que la matriz asociada a la correspondiente aplicaci´on lineal sea la matriz diagonal Λ (en lugar de la matriz de partida Σ). Las columnas de la matriz H son precisamente los autovectores linealmente independientes (y de m´odulo unitario) de la matriz Σ, es decir, d vectores linealmente independientes, x1 , x2 , . . . , xd , tales que xi 0 xi = 1 para todo i = 1, 2, . . . , n y con xi 0 xj = 0 si i 6= j, verificando adem´as que ∃λi ∈ R tal que Σxi = λi xi (condici´on de ser un autovector). Adem´as, los autovalores λ1 , λ2 , . . ., λd (que son todos positivos) son precisamente los elementos de la diagonal de la matriz Λ. Partiendo de una variable Z ∈ Nd (0, I) (f´acilmente simulable a partir de Z1 , Z2 , . . ., Zd ∼ N (0, 1) independientes), se tiene que Λ1/2 Z ∈ Nd (0, Λ), siendo  Λ

1/2

  = 

1/2

λ1 0 .. .

0 λ2 .. .

0

0

1/2

··· ··· ...

0 0 .. . 1/2

   . 

· · · λd

Multiplicando por la izquierda por la matriz H, se tiene HΛ1/2 Z ∈ Nd (0, HΛH0 ) ∈ Nd (0, Σ) Finalmente, basta sumar el vector µ para obtener X = µ + HΛ1/2 Z ∈ Nd (µ , Σ) Una vez obtenidos los autovalores, λ1 , λ2 , . . ., λd , y los autovectores asociados de la matriz Σ, que determinan las columnas de la matriz H, el algoritmo proceder´ıa como sigue: 1. Simular Z1 , Z2 , . . . , Zd ∼ N (0, 1)√ independientes. 2. Para i = 1, 2, . . . , d hacer Yi = λi Zi . 3. Devolver X = µ + HY.

´ ´ NORMAL MULTIVARIANTE73 6.4. METODOS PARA SIMULAR LA DISTRIBUCION Ejemplo 6.4.1 Dar un algoritmo para simular la distribuci´ on     1 2,36 −0,48 N2 , 3 −0,48 2,64 Para encontrar los autovalores y autovectores de Σ resolvemos det (Σ − λI) = 0, es decir, 2,36 − λ −0,48 2 −0,48 2,64 − λ = 0 ⇔ (2,36 − λ) (2,64 − λ) − (−0,48) = 0 √ 5 ± 52 − 6 · 4 2 ⇔ λ − 5λ + 6 = 0 ⇔ λ = 2 que ofrece como soluciones λ1 = 3 y λ2 = 2. Para encontrar autovectores de m´ odulo 1 correspondientes a esos autovalores no tenemos m´as que resolver los sistemas (Σ − λi I) µ = 0 para i = 1, 2 imponiendo la condici´on de m´odulo igual a 1, es decir x21 + x22 = 1. As´ı, resulta     -0.64 -0.48 16 12 Σ − λ1 I = = -0.04 , luego -0.48 -0.36 12 9 4 (Σ − λ1 I) x = 0 ⇔ x2 = − x1 , pero como x21 + x22 = 1, se tiene 3 25 2 3 4 x1 = 1, luego x1 = y x2 = − 9 5 5 (tambi´en es soluci´ on si cambiamos ambos de signo); 

   0.36 -0.48 9 −12 Σ − λ2 I = = 0.04 , luego -0.48 0.64 −12 16 3 (Σ − λ2 I) x = 0 ⇔ x2 = x1 , pero como x21 + x22 = 1, se tiene 4 25 2 4 3 x1 = 1, luego x1 = y x2 = 16 5 5 De esta forma, la matriz H resulta, entre otras posibilidades,     3/5 4/5 0.6 0.8 H= = . −4/5 3/5 -0.8 0.6 Ahora Y=Λ

1/2

 √ Z=

3 √0 2 0



Z1 Z2



 √  3Z 1 = √ 2Z2

y finalmente,  X = µ + HY =

1 + 0.6Y1 + 0.8Y2 3 − 0.8Y1 + 0.6Y2



As´ı, el algoritmo resultar´ıa 1. Simular Z1 ,√ Z2 ∼ N (0, 1) √ independientes. 2. Hacer Y1 = 3Z1 e Y2 = 2Z2 . 3. Obtener X1 = 1 + 0.6Y1 + 0.8Y2 y X2 = 3 − 0.8Y1 + 0.6Y2 . 4. Devolver X = (X1 , X2 )0 .

´ DE DISTRIBUCIONES MULTIDIMENSIONALES 74CAP´ITULO 6. SIMULACION

Cap´ıtulo 7 Dise˜ no de experimentos de simulaci´ on En el presente cap´ıtulo se abordar´an algunas de las cuestiones m´as importantes a la hora de dise˜ nar un estudio de simulaci´on: Similitudes y diferencias entre la simulaci´on y la experimentaci´on sobre el sistema real. Simulaci´on est´atica y din´amica. Simulaci´on por eventos y por cuantos. T´ecnicas de reducci´on de la varianza. Problemas de estabilizaci´on y dependencia.

7.1.

Diferencias y similitudes con la experimentaci´ on real

Teniendo en cuenta que la simulaci´on es la t´ecnica consistente en la realizaci´on de experimentos de muestreo sobre un modelo construido a partir de un sistema real, es obvio que la simulaci´on necesitar´a de gran cantidad de t´ecnicas estad´ısticas para obtener las muestras (muestreo) y para analizar los resultados obtenidos por la experimentaci´on artificial (estimaci´on, intervalos de confianza, contrastes de hip´otesis, etc.). Por todo ello, puede afirmarse que, en general, en cuanto a la utilizaci´on de t´ecnicas estad´ısticas es muy similar a la propia experimentaci´on sobre el sistema real. Entre las diferencias caben destacar las siguientes: 1. La utilizaci´on de t´ecnicas de estimaci´on puntual, construcci´on de intervalos de confianza y contrastes de hip´otesis es algo menos frecuente en la simulaci´on que en la experimentaci´on real. La raz´on es que algunos de los par´ametros (los de control) ya son conocidos en la simulaci´on y, por tanto, no es necesario hacer inferencia sobre ellos, aunque s´ı sobre los de salida (que miden, de alguna forma, el comportamiento del sistema).

75

˜ DE EXPERIMENTOS DE SIMULACION ´ CAP´ITULO 7. DISENO

76

2. La simulaci´on suele hacer un uso mucho m´as intensivo de t´ecnicas de ordenaci´on y optimizaci´on. Esto es debido a que, en el contexto de la simulaci´on, es factible comparar un gran n´ umero de escenarios (entre los que se desea optimizar, por ejemplo) en muy poco tiempo, cosa que se da muy raramente en la experimentaci´on real. 3. Una peculiaridad de la simulaci´on es que casi siempre es posible comparar distintas estrategias sobre las mismas muestras simuladas (simplemente utilizando la misma semilla en la simulaci´on, convenientemente planificada).

7.2.

Simulaci´ on est´ atica y din´ amica

La simulaci´on se dice est´atica si en el modelo no juega ning´ un papel el transcurso del tiempo mientras que es din´amica si el tiempo es una de las variables importantes del modelo. La simulaci´on est´atica se usa muy frecuentemente por los estad´ısticos para comprobar el comportamiento comparativo de diversos m´etodos estad´ısticos alternativos para tama˜ nos muestrales finitos (complementando los estudios te´oricos, casi siempre asint´oticos). En la simulaci´on din´amica, normalmente, se trata de ir analizando los distintos estados por los que va pasando un sistema que evoluciona en el tiempo. Esto provoca, en general, un mayor coste computacional y problemas de estabilizaci´on y dependencia. Existen dos grandes tipos de simulaci´on din´amica: la simulaci´ on continua, en la que se supone que el sistema cambia de estado constantemente, y la simulaci´ on discreta, para la cual los cambios se producen en ciertos instantes de tiempo singulares.

7.3.

Simulaci´ on por eventos y por cuantos

Con el nombre de simulaci´on por eventos, o as´ıncrona, designamos el tipo de simulaci´on din´amica discreta en la cual se controla la variable tiempo movi´endola hasta la ocurrencia del siguiente suceso (o evento). Esto implica la necesidad de controlar minuciosamente cu´al es dicho pr´oximo suceso: saber cu´ales son los posibles sucesos en un futuro inmediato y cu´al de ellos es el m´as inmediato. La simulaci´on por cuantos, responde a una filisof´ıa totalmente diferente. Se trata de examinar el sistema (que evoluciona en el tiempo) dejando pasar peque˜ nos intervalos de tiempo de longitud δ, fija, (llamados cuantos) en los cuales se supone que, a lo sumo, un solo suceso puede producirse. En general, la simulaci´on por eventos es exacta y de m´as dif´ıcil implementaci´on, pero de mucha m´as r´apida ejecuci´on que la simulaci´on por cuantos. Sin embargo esta u ´ltima es muchas veces la u ´nica posibilidad factible en la simulaci´on din´amica continua.

´ ´ DE LA VARIANZA 7.4. TECNICAS DE REDUCCION

7.4.

77

T´ ecnicas de reducci´ on de la varianza

Existen multitud de t´ecnicas encaminadas a reducir la varianza en un estudio de simulaci´on o bien a tratar de estimarla. Algunas de ellas son el uso de n´ umeros aleatorios comunes, la utilizaci´on de variables antit´eticas, la estratificaci´on, el uso de variables de control, el m´etodo Jackknife, los m´etodos de remuestreo (destacando entre ellos el m´etodo bootstrap), etc. En general conviene tener en cuenta que si uno de los objetivos de la simulaci´on es precisamente estimar la variabilidad, no conviene utilizar estas t´ecnicas de reducci´on ´ de la varianza. Estas son aplicables normalmente cuando la simulaci´on pretende ofrecer respuestas, lo m´as precisas posibles, s´olo sobre cantidades medias.

7.4.1.

N´ umeros aleatorios comunes

Sup´ongase que se desea comparar dos estrategias distintas, X e Y, mediante N repeticiones (o “trials”) de un experimento de simulaci´on, de las cuales se han obtenido los valores num´ericos de las variables de salida X1 , X2 , . . . , XN , para la primera, e Y1 , Y2 , . . ., YN , para la segunda. Si la comparaci´on se realiza estimando la diferencia de las medias de las variables de salida para ambas estrategias, E (X) − E (Y ) = E (X − Y ) , puede usarse X − Y = N 1 X (Xi − Yi ), cuya varianza viene dada por N i=1 V ar X − Y




N 1 X 1 = V ar (Xi − Yi ) = V ar (X1 − Y1 ) 2 N i=1 N

=

1 (V ar (X1 ) + V ar (Y1 ) − 2Cov (X1 , Y1 )) N

Usando los mismos n´ umeros aleatorios (es decir, repitiendo los c´alculos con la misma semilla) en las variables de entrada de la simulaci´on, se tiene que Cov (Xi , Yi ) > 0 y, por tanto,
1 (V ar (X1 ) + V ar (Y1 )) V ar X − Y ≤ N que es la varianza que tendr´ıa X − Y en caso de haber usado muestras independientes para cada estrategia.

78

7.4.2.

˜ DE EXPERIMENTOS DE SIMULACION ´ CAP´ITULO 7. DISENO

Variables antit´ eticas

Sup´ongase ahora que se desea evaluar el resultado de una u ´nica estrategia (sin compararla con ninguna otra alternativa). Despu´es de N repeticiones de la simulaci´on, tendremos N valores num´ericos de las variables X1 , X2 , . . ., XN , procediendo a estiN 1 X mar la media E (X) te´orica mediante X = Xi . Dado que ´este es un estimador N i=1
insesgado, su precisi´on puede medirse mediante V ar X . N
1 X Si las variables son independientes, la V ar X = 2 V ar (Xi ), mientras que, N i=1 en general, se tiene ! N X
1 X V ar X = 2 V ar (Xi ) + 2 Cov (Xi , Xj ) N i=1 i