Simulación de Monte Carlo

Introducción Algunos Resultados Estimar una Medida Integración de Monte Carlo Muestreador de Gibbs Referencia Simulación Estocástica Simulación de Monte Carlo Carlos V. Ramı́rez Ibáñez Carlos V. Ramı́rez Ibáñez Universidad La Salle 1 / 44 Introducción Algunos Resultados Estimar una Medida Integración de Monte Carlo Muestreador de Gibbs Referencia Contenido 1 Introducción 2 Algunos Resultados 3 Estimar una Medida Estimando un área Algoritmo Estimando π 4 Integración de Monte Carlo Construcción Estimador de Monte Carlo Algoritmo Ejemplos Integración Múltiple 5 Muestreador de Gibbs Caso Bivariado Algoritmo Ejemplos 6 Referencias Carlos V. Ramı́rez Ibáñez Universidad La Salle 2 / 44 Introducción Algunos Resultados Estimar una Medida Integración de Monte Carlo Muestreador de Gibbs Referencia Introducción ∗ Es una técnica que permite dar solución a problemas matemáticos que son complicados o incluso imposibles de resolver de forma analı́tica. ∗ Surge entre 1944 y 1946, su descubrimiento lo comparten Ulam y John Von Neumann. ∗ Usado por cientı́ficos para desarrollar la Bomba Atómica durante la Segunda Guerra Mundial. ∗ El término “Monte Carlo” hace referencia al Casino de Monte Carlo en Mónaco. Figura: Casino de Monte Carlo Carlos V. Ramı́rez Ibáñez Universidad La Salle 3 / 44 Introducción Algunos Resultados Estimar una Medida Integración de Monte Carlo Muestreador de Gibbs Referencia Algunos Resultados Algunos conceptos y resultados que utilizaremos a lo largo de la presentación: 1. Función de densidad condicional. 2. Ley Fuerte de los Grandes Números. 3. Teorema del Lı́mite Central. 4. Generador Congruencial Combinado Múltiple. 5. Generador MRG32k3a. 6. Método de la Transformada Inversa Carlos V. Ramı́rez Ibáñez Universidad La Salle 4 / 44 Introducción Algunos Resultados Estimar una Medida Integración Estimando de Monte un área Carlo Algoritmo Muestreador Estimando de Gibbs π Referencia Estimando un área Sean a1 , a2 , b1 , b2 ∈ R+ tales que a1 < b1 y a2 < b2 . Supongamos que estamos interesados en calcular el área de una región D en R2 tal que D ⊂ G := [a1 , b1 ] × [a2 , b2 ] Gráficamente b2 D G a2 a1 Carlos V. Ramı́rez Ibáñez b1 Universidad La Salle 5 / 44 Introducción Algunos Resultados Estimar una Medida Integración Estimando de Monte un área Carlo Algoritmo Muestreador Estimando de Gibbs π Referencia Estimando un área Consideremos dos variables aleatorias continuas e independientes X y Y tales que X ∼ U(a1 , b1 ) y Y ∼ U(a2 , b2 ), por lo que V = (X, Y ) es un vector aleatorio con V ∼ U((a1 , b1 ) × (a2 , b2 )) ‘ Entonces la función de densidad conjunta de V está dada por fV (x, y) = fX (x)fY (y)    1 1 = b1 − a 1 b2 − a 2 1 = (b1 − a1 )(b2 − a2 ) 1 = A(G) donde A(G) denota el área de G. Carlos V. Ramı́rez Ibáñez Universidad La Salle 6 / 44 Introducción Algunos Resultados Estimar una Medida Integración Estimando de Monte un área Carlo Algoritmo Muestreador Estimando de Gibbs π Referencia Estimando un área Como D ⊂ G, entonces podemos calcular probabilidades de V como ZZ fV (x, y)dA P[V ∈ D] = Z ZD 1 = dA A(G) D ZZ 1 dA = A(G) D A(D) = A(G) esto es A(D) = A(G) P[V ∈ D] A(G) es fácil de calcular pues G es un rectángulo, entonces su área está dada por A(G) = (b1 − a1 )(b2 − a2 ) El problema es hallar P[V ∈ D]. Carlos V. Ramı́rez Ibáñez Universidad La Salle 7 / 44 Introducción Algunos Resultados Estimar una Medida Integración Estimando de Monte un área Carlo Algoritmo Muestreador Estimando de Gibbs π Referencia Estimando P[V ∈ D] Supongamos que una región D en R2 está contenida en el rectángulo G como se muestra a continuación: b2 D G a2 a1 b1 Como V ∼ U((a1 , b1 ) × (a2 , b2 )), entonces podemos generar números pseudoaleatorios de V en el rectángulo G, a la cantidad de números pseudoaleatorios que generemos le llamaremos ensayos. Carlos V. Ramı́rez Ibáñez Universidad La Salle 8 / 44 Introducción Algunos Resultados Estimar una Medida Integración Estimando de Monte un área Carlo Algoritmo Muestreador Estimando de Gibbs π Referencia Estimando P[V ∈ D] Con n ensayos, veremos cuántos caen dentro de la región D, pensaremos esta cantidad como el número de éxitos. Entonces podemos estimar P[V ∈ D] como número de éxito P[V ∈ D] ∼ = número de ensayos Si n → ∞, entonces obtendremos una aproximación buena de la probabilidad real. Por lo tanto A(D) = A(G) P[V ∈ D] ∼ = (b1 − a1 )(b2 − a2 ) Carlos V. Ramı́rez Ibáñez  número de éxito número de ensayos Universidad La Salle  9 / 44 Introducción Algunos Resultados Estimar una Medida Integración Estimando de Monte un área Carlo Algoritmo Muestreador Estimando de Gibbs π Referencia Algoritmo Si uno desea hallar el área de una región D ⊂ G con G = [a1 , b1 ]×[a2 , b2 ], entonces puede seguir el siguiente algoritmo: 1. Generar una muestra pseudoaleatoria de V ∼ U((a1 , b1 ) × (a2 , b2 )) de tamaño n, esta muestra la denotaremos por {vi }ni=1 . 2. Hallar el número de éxitos, es decir, obtener el número de vi′ s ∈ D. 3. Estimar el área de D con (b1 − a1 )(b2 − a2 ) Carlos V. Ramı́rez Ibáñez  número de éxito número de ensayos Universidad La Salle  10 / 44 Introducción Algunos Resultados Estimar una Medida Integración Estimando de Monte un área Carlo Algoritmo Muestreador Estimando de Gibbs π Referencia Estimando π Sean D y G dos regiones en el plano descritas por D = {(x, y) ∈ R2 : (x − 1)2 + (y − 1)2 ≤ 9, 1 ≤ x, y ≤ 4} G = {(x, y) ∈ R2 : 1 ≤ x, y ≤ 4} es decir Figura: Regiones G y D Carlos V. Ramı́rez Ibáñez Universidad La Salle 11 / 44 Introducción Algunos Resultados Estimar una Medida Integración Estimando de Monte un área Carlo Algoritmo Muestreador Estimando de Gibbs π Referencia Estimando π Para estimar el valor de π tenemos que estimar el área de D, para hacerlo consideremos que   número de éxito ∼ (b1 − a1 )(b2 − a2 ) A(D) = número de ensayos Como A(D) = Entonces 4A(D) πr2 ⇐⇒ π = 4 r2 4(b1 − a1 )(b2 − a2 ) π∼ = r2  número de éxito número de ensayos  Obtengamos una muestra de tamaño n del vector V = (X, Y ) ∼ U((1, 4)2 ) Carlos V. Ramı́rez Ibáñez Universidad La Salle 12 / 44 Introducción Algunos Resultados Estimar una Medida Integración Estimando de Monte un área Carlo Algoritmo Muestreador Estimando de Gibbs π Referencia Estimando π Teniendo {vi }ni=1 , determinemos el número de vi′ s ∈ D, es decir, el número de éxitos Éxitos = #{vi : (xi − 1)2 + (yi − 1)2 ≤ 9, 1 ≤ xi , yi ≤ 4} Finalmente calculamos 4(b1 − a1 )(b2 − a2 ) π∼ = r2  número de éxito número de ensayos  Por ejemplo, para n = 500 se obtuvo n Éxitos A(D) Pi Monte Carlo Error 500 389 7.002 3.112 0.02959265 Carlos V. Ramı́rez Ibáñez Universidad La Salle 13 / 44 Introducción Algunos Resultados Estimar una Medida Integración Estimando de Monte un área Carlo Algoritmo Muestreador Estimando de Gibbs π Referencia Estimando π El resultado anterior puede verse gráficamente como sigue Figura: Estimación de π con n = 500 Carlos V. Ramı́rez Ibáñez Universidad La Salle 14 / 44 Introducción Algunos Resultados Estimar una Medida Integración Estimando de Monte un área Carlo Algoritmo Muestreador Estimando de Gibbs π Referencia Estimando π Para distintos valores de n se obtuvieron los siguientes resultados: n Éxitos A(D) Pi Monte Carlo Error 1, 000 2, 000 3, 000 5, 000 10, 000 70, 000 100, 000 150, 000 200, 000 786 1, 574 2, 359 3, 944 7, 814 55, 028 78, 668 117, 866 157, 063 7.074 7.083 7.077 7.0992 7.0326 7.075029 7.08012 7.07196 7.067835 3.144 3.148 3.145333 3.1552 3.1256 3.144457 3.14672 3.143093 3.14126 0.002407346 0.006407346 0.00374068 0.01360735 0.01599265 0.002864489 0.005127346 0.00150068 0.0003326536 Carlos V. Ramı́rez Ibáñez Universidad La Salle 15 / 44 Introducción Algunos Resultados Estimar una Medida Integración Estimando de Monte un área Carlo Algoritmo Muestreador Estimando de Gibbs π Referencia Estimando de π Gráficamente (a) n = 3, 000 (b) n = 5, 000 Figura: Estimación de π Carlos V. Ramı́rez Ibáñez Universidad La Salle 16 / 44 Introducción Algunos Resultados Estimar una Medida Integración Estimando de Monte un área Carlo Algoritmo Muestreador Estimando de Gibbs π Referencia Estimando de π Figura: Estimación de π con n = 100, 000 Carlos V. Ramı́rez Ibáñez Universidad La Salle 17 / 44 Introducción Algunos Resultados Estimar una Medida Integración Construcción de Monte Estimador Carlo Muestreador de Monte Carlo de Gibbs Algoritmo Referencia Ejem Integración de Monte Carlo Podemos utilizar simulación de Monte Carlo para estimar integrales que no pueden ser evaluadas de forma analı́tica. La integración de Monte Carlo es muy utilizada para estimar integrales de la forma Z Z ··· f (x1 , . . . , xm )dx1 · · · dxm Ω donde Ω es un subconjunto en Rm . Carlos V. Ramı́rez Ibáñez Universidad La Salle 18 / 44 Introducción Algunos Resultados Estimar una Medida Integración Construcción de Monte Estimador Carlo Muestreador de Monte Carlo de Gibbs Algoritmo Referencia Ejem Integración de Monte Carlo Supongamos que g : [a, b] → R es una función integrable en (a, b) y estamos interesados en evaluar la integral I= Z b g(x)dx a donde g no es analı́ticamente integrable. Construiremos un expresión para poder estimar esta integral. Carlos V. Ramı́rez Ibáñez Universidad La Salle 19 / 44 Introducción Algunos Resultados Estimar una Medida Integración Construcción de Monte Estimador Carlo Muestreador de Monte Carlo de Gibbs Algoritmo Referencia Ejem Construcción Sean a, b ∈ R con a < b y X una variable aleatoria tal que X ∼ U(a, b), entonces la función de densidad de X, denotada por fX (x), está dada por   1 si x ∈ (a, b) fX (x) = b − a 0 si x ∈ / (a, b) Y sea g : [a, b] → R una función integrable en (a, b) tal que g(X) es una variable aleatoria con esperanza finita, entonces Z g(x)fX (x)dx E[g(X)] = = Z R b g(x) a 1 = b−a esto es I= Z Z 1 dx b−a b g(x)dx a b a g(x)dx = (b − a) E[g(X)] Carlos V. Ramı́rez Ibáñez Universidad La Salle 20 / 44 Introducción Algunos Resultados Estimar una Medida Integración Construcción de Monte Estimador Carlo Muestreador de Monte Carlo de Gibbs Algoritmo Referencia Ejem Construcción Si consideramos X1 , X2 , . . . , Xn una muestra de variables aleatorias independientes e idénticamente distribuidas, entonces podemos transformar cada una de estas con la función g para ası́ obtener la muestra g(X1 ), g(X2 ), . . . , g(Xn ) y por la Ley Fuerte de los Grandes Números n g(X)n := 1X c.s. g(Xi ) −−→ E[g(X)] n i=1 cuando n → ∞, es decir " # n 1X P lı́m g(Xi ) = E[g(X)] = 1 n→∞ n i=1 Por lo tanto Z b a n g(x)dx ∼ = (b − a)g(X)n = Carlos V. Ramı́rez Ibáñez b−aX g(Xi ) n i=1 Universidad La Salle 21 / 44 Introducción Algunos Resultados Estimar una Medida Integración Construcción de Monte Estimador Carlo Muestreador de Monte Carlo de Gibbs Algoritmo Referencia Ejem Estimador de Monte Carlo Al estadı́stico n (b − a)g(X)n = b−aX g(Xi ) n i=1 se le denota por Ib y se le conoce como el estimador de Monte Carlo de I. Y no es difı́cil ver que h i E Ib = I y   σ2 Var Ib = n donde σ 2 denota la varianza de la variable aleatoria g(X). Carlos V. Ramı́rez Ibáñez Universidad La Salle 22 / 44 Introducción Algunos Resultados Estimar una Medida Integración Construcción de Monte Estimador Carlo Muestreador de Monte Carlo de Gibbs Algoritmo Referencia Ejem Estimador de Monte Carlo Dado que σ 2 es desconocida, entonces esta debe estimarse, se considera 2 1 X g(Xi ) − Ib n − 1 i=1 n S2 = Por el Teorema del Lı́mite Central Ib − I d √ −−−−−→ N (0, 1) n→+∞ S/ n Entonces un intervalo de (1 − α)100 % de confianza para I es S S Ib − z1−α/2 √ ≤ I ≤ Ib + z1−α/2 √ n n Carlos V. Ramı́rez Ibáñez Universidad La Salle 23 / 44 Introducción Algunos Resultados Estimar una Medida Integración Construcción de Monte Estimador Carlo Muestreador de Monte Carlo de Gibbs Algoritmo Referencia Ejem Algoritmo Si se desea evaluar la integral I= Z b g(x)dx a uno puede seguir el siguiente algoritmo: 1. Generar una muestra pseudoaleatoria de tamaño n de X ∼ U(a, b). 2. Transformar cada observación de la muestra como g(xi ). 3. Calcular el estimador de Montecarlo n b−aX g(xi ) Ib = n i=1 Carlos V. Ramı́rez Ibáñez Universidad La Salle 24 / 44 Introducción Algunos Resultados Estimar una Medida Integración Construcción de Monte Estimador Carlo Muestreador de Monte Carlo de Gibbs Algoritmo Referencia Ejem Ejemplo 1 Consideremos la integral Z 3π 2 cos(x)dx 0 Esta integral puede evaluarse analı́ticamente utilizando el Teorema Fundamental del Cálculo 3π   Z 3π 2 2 3π − sen(0) = −1 = sen cos(x)dx = sen(x) 2 0 0 Utilizando el generador MRG32k3a con una semilla x0 = 4 se obtuvo: n Integración Montecarlo Error 10 100 1,000 10,000 100,000 -1.7118879 -1.5958955 -0.990388585 -1.06288631 -1.006874681 0.7118879 0.5958955 0.009611415 0.06288631 0.006874681 Carlos V. Ramı́rez Ibáñez Universidad La Salle 25 / 44 Introducción Algunos Resultados Estimar una Medida Integración Construcción de Monte Estimador Carlo Muestreador de Monte Carlo de Gibbs Algoritmo Referencia Ejem Estimando la integral gaussiana Consideremos la integral gaussiana Z ∞ 2 e−x dx −∞ existen distintas técnicas “no elementales” para evaluar esta integral, en cualquiera de ellas puede demostrarse que Z ∞ √ 2 e−x dx = π −∞ Para utilizar Integración de Monte Carlo puede notarse que el problema son los lı́mites de integración. Una posible opción (aunque de momento inútil) es notando que el integrando es una función par, por lo que Z ∞ Z ∞ 2 2 e−x dx e−x dx = 2 −∞ Carlos V. Ramı́rez Ibáñez 0 Universidad La Salle 26 / 44 Introducción Algunos Resultados Estimar una Medida Integración Construcción de Monte Estimador Carlo Muestreador de Monte Carlo de Gibbs Algoritmo Referencia Ejem Estimando la integral gaussiana El problema con los lı́mites de integración sigue existiendo por lo que es natural pensar en hacer un cambio de variable, si consideramos u= entonces x= Por lo que Z ∞ e −∞ −x2 1−u u dx = 2 y Z 1 0 1 1+x dx = − 1 du u2   (1 − u)2 1 exp − du u2 u2 Ahora ya podemos utilizar Integración de Monte Carlo. Carlos V. Ramı́rez Ibáñez Universidad La Salle 27 / 44 Introducción Algunos Resultados Estimar una Medida Integración Construcción de Monte Estimador Carlo Muestreador de Monte Carlo de Gibbs Algoritmo Referencia Ejem Estimando la integral gaussiana Para distintos valores √ de n, utilizando el generador MRG32k3a con semilla x0 = 9 y considerando que π ∼ = 1.772454 obtenemos n Integración Montecarlo Error 800 1,000 2,000 4,000 5,000 10,000 70,000 100,000 1.822809 1.833958 1.833301 1.77081 1.755689 1.744135 1.763236 1.767595 0.05035508 0.06150422 0.06084676 0.001643709 0.01676474 0.02831926 0.009217844 0.004859035 Carlos V. Ramı́rez Ibáñez Universidad La Salle 28 / 44 Introducción Algunos Resultados Estimar una Medida Integración Construcción de Monte Estimador Carlo Muestreador de Monte Carlo de Gibbs Algoritmo Referencia Ejem Estimando la integral gaussiana Figura: Estimación de la integral para distintos valores de n Carlos V. Ramı́rez Ibáñez Universidad La Salle 29 / 44 Introducción Algunos Resultados Estimar una Medida Integración Construcción de Monte Estimador Carlo Muestreador de Monte Carlo de Gibbs Algoritmo Referencia Ejem Estimando la integral gaussiana Figura: Estimación de la integral para distintos valores de n Carlos V. Ramı́rez Ibáñez Universidad La Salle 30 / 44 Introducción Algunos Resultados Estimar una Medida Integración Construcción de Monte Estimador Carlo Muestreador de Monte Carlo de Gibbs Algoritmo Referencia Ejem Integración Múltiple Supóngase que ahora se está interesado en evaluar la integral Z β2 Z α2 Z ϕ2 g(x1 , x2 , . . . , xm )dx1 dx2 · · · dxm ··· I= ϕ1 β1 α1 entonces uno puede seguir el siguiente algoritmo: 1. Generar una muestra pseudoaleatoria de tamaño n del vector (X1 , X2 , . . . , Xm ) ∼ U ([α1 , α2 ] × [β1 , β2 ] × · · · × [ϕ1 , ϕ2 ]) A la muestra pseudoaleatoria de tamaño n proveniente de la variable aleatoria Xj (j = 1, 2, . . . , m) se le denotará por {xij }ni=1 . 2. Trasformar cada observación de la muestra como g(xi1 , xi2 , . . . , xim ) i = 1, 2, . . . , n 3. Calcular el estimador de Montecarlo n (α2 − α1 )(β2 − β1 ) · · · (ϕ2 − ϕ1 ) X g(xi1 , xi2 , . . . , xim ) Ib = n i=1 Carlos V. Ramı́rez Ibáñez Universidad La Salle 31 / 44 Introducción Algunos Resultados Estimar una Medida Integración Caso Bivariado de MonteAlgoritmo Carlo Muestreador Ejemplos de Gibbs Referencia Muestreador de Gibbs Es uno de los algoritmos que pertenecen a los conocidos Métodos de Monte Carlo vı́a Cadenas de Markov. Sea X = (X1 , X2 , . . . , Xn ) un vector aleatorio con función de densidad conjunta fX1 ,X2 ,...,Xn (x1 , x2 , . . . , xn ) y supongamos que estamos interesados en obtener una muestra de X. Este algoritmo es útil cuando no se puede obtener una muestra directamente de la distribución pero sı́ sabemos cuáles son todas las distribuciones condicionales y podemos obtener muestras de ellas. Carlos V. Ramı́rez Ibáñez Universidad La Salle 32 / 44 Introducción Algunos Resultados Estimar una Medida Integración Caso Bivariado de MonteAlgoritmo Carlo Muestreador Ejemplos de Gibbs Referencia Caso Bivariado Supongamos que estamos interesados en generar una muestra del vector aleatorio continuo (X, Y ) donde fX,Y (x, y) es la función de densidad conjunta asociada. Y supongamos que conocemos la distribución de las variables aleatorias X|Y = y y Y |X = x y las funciones de densidad condicionales fX|Y (x|y) y fY |X (x|y). Carlos V. Ramı́rez Ibáñez Universidad La Salle 33 / 44 Introducción Algunos Resultados Estimar una Medida Integración Caso Bivariado de MonteAlgoritmo Carlo Muestreador Ejemplos de Gibbs Referencia Algoritmo Para generar una muestra pseudoaleatoria de un vector aleatorio (X, Y ) con función de densidad conjunta fX,Y (x, y) se sigue el siguiente algoritmo: 1. Se define x0 . 2. Generar yi |xi−1 ∼ fY |X (y|xi−1 ). 3. Generar xi |yi ∼ fX|Y (x|yi ). 4. Se repiten los pasos 2 y 3 hasta generar la muestra del tamaño deseado. Carlos V. Ramı́rez Ibáñez Universidad La Salle 34 / 44 Introducción Algunos Resultados Estimar una Medida Integración Caso Bivariado de MonteAlgoritmo Carlo Muestreador Ejemplos de Gibbs Referencia Ejemplo 1 Consideremos el vector aleatorio continuo (X, Y ) con distribución normal bivariada, entonces su función de densidad conjunta es   1 (x − µ)T Σ−1 (x − µ) √ fX,Y (x, y) = exp − 2 2π det Σ para x, y ∈ R donde µ=  µX µY  y Σ=  2 σX ρσX σY ρσX σY σY2  siendo ρ el coeficiente de correlación. Consideremos cuando σX = σY = 1 y µX = µY = 0, en tal caso   2 1 x − 2ρxy + y 2 p fX,Y (x, y) = exp − 2(1 − ρ2 ) 2π 1 − ρ2 Carlos V. Ramı́rez Ibáñez Universidad La Salle 35 / 44 Introducción Algunos Resultados Estimar una Medida Integración Caso Bivariado de MonteAlgoritmo Carlo Muestreador Ejemplos de Gibbs Referencia Ejemplo 1 Notemos que  x2 − 2ρxy + y 2 p fX,Y (x, y) = exp − 2(1 − ρ2 ) 2π 1 − ρ2  2  1 x − 2ρxy + y 2 (1 − ρ2 + ρ2 ) p = exp − 2(1 − ρ2 ) 2π 1 − ρ2   2 1 x − 2ρxy + ρ2 y 2 + y 2 (1 − ρ2 ) p = exp − 2(1 − ρ2 ) 2π 1 − ρ2   1 (x − ρy)2 y2 p = exp − − 2(1 − ρ2 ) 2 2π 1 − ρ2   2  1 1 y (x − ρy)2 p = √ exp − exp − 2 2(1 − ρ2 ) 2π 2π(1 − ρ2 ) {z }| | {z }  1 ∗ Carlos V. Ramı́rez Ibáñez ∗∗ Universidad La Salle 36 / 44 Introducción Algunos Resultados Estimar una Medida Integración Caso Bivariado de MonteAlgoritmo Carlo Muestreador Ejemplos de Gibbs Referencia Ejemplo 1 Notemos que ∗ corresponde a la función de densidad de Y ∼ N (0, 1) y ∗∗ corresponde a la función de densidad de X|Y = y ∼ N (ρy, 1 − ρ2 ), es decir fX,Y (x, y) = fY (y)fX|Y (x|y) por lo que ya sabemos quién es fX|Y (x|y) y similarmente sabremos quién es fY |X (y|x). Por lo tanto X|Y = y ∼ N (ρy, 1 − ρ2 ) Y |X = x ∼ N (ρx, 1 − ρ2 ) Carlos V. Ramı́rez Ibáñez Universidad La Salle 37 / 44 Introducción Algunos Resultados Estimar una Medida Integración Caso Bivariado de MonteAlgoritmo Carlo Muestreador Ejemplos de Gibbs Referencia Ejemplo 1 Por lo que para generar muestras de esta distribución, el algoritmo se reduce a: 1. Se define y0 . 2. Generar xi |yi−1 ∼ N (ρyi−1 , 1 − ρ2 ). 3. Generar yi |xi ∼ N (ρxi , 1 − ρ2 ). 4. Se repiten los pasos 2 y 3 hasta generar la muestra del tamaño deseado. Considerando y0 = 1/5, el generador de R con semilla 1998, entonces los primeros cinco pares generados del vector (X, Y ) con distribución normal bivariada con σX = σY = 1, µX = µY = 0 y ρ = 1/4 son: 1 2 3 4 5 xi yi -1.6870310 1.0040635 0.6992634 -0.3308307 0.6162344 -1.94417194 -2.40980887 -0.42746603 -0.02150579 -0.99492508 Carlos V. Ramı́rez Ibáñez Universidad La Salle 38 / 44 Introducción Algunos Resultados Estimar una Medida Integración Caso Bivariado de MonteAlgoritmo Carlo Muestreador Ejemplos de Gibbs Referencia Ejemplo 2 Consideremos el vector aleatorio (X, Y ) con función de densidad bivariada r   1 y λy(x − µ)2 λ α−1 fX,Y (x, y) = y exp − − Γ(α)β α 2π β 2 para x ∈ R y y > 0. Puede demostrarse que (x − µ)2 exp − fX|Y (x|y) = q 1 1 ) 2( λy ) 2π( λy 1 fY |X (y|x) = ! 1 y α−1 e−y/β Γ(α)β α son las funciones de densidad de las variables X|Y = y y Y |X = x respectivamente. Carlos V. Ramı́rez Ibáñez Universidad La Salle 39 / 44 Introducción Algunos Resultados Estimar una Medida Integración Caso Bivariado de MonteAlgoritmo Carlo Muestreador Ejemplos de Gibbs Referencia Ejemplo 2 Por lo tanto X|Y = y ∼ N  1 µ, λy Y |X = x ∼ Γ(α, β)  En este caso el algoritmo se reduce a: 1. Se define x0 . 2. Generar yi |xi−1 ∼ Γ(α, β).   3. Generar xi |yi ∼ N µ, λy1 i . 4. Se repiten los pasos 2 y 3 hasta generar la muestra del tamaño deseado. Carlos V. Ramı́rez Ibáñez Universidad La Salle 40 / 44 Introducción Algunos Resultados Estimar una Medida Integración Caso Bivariado de MonteAlgoritmo Carlo Muestreador Ejemplos de Gibbs Referencia Ejemplo 2 Considerando α = 2, β = 3, µ = 0, λ = 4, x0 = 0.45 y el generador de R con semilla 1998, entonces los primeros cinco pares generados del vector (X, Y ) con la función de densidad mencionada son: 1 2 3 4 5 xi yi -0.39730305 -0.23157286 -0.26751160 -0.72914950 -0.02682224 0.7809389 0.5624251 0.7040122 0.4134725 0.7074945 Carlos V. Ramı́rez Ibáñez Universidad La Salle 41 / 44 Introducción Algunos Resultados Estimar una Medida Integración Caso Bivariado de MonteAlgoritmo Carlo Muestreador Ejemplos de Gibbs Referencia Ejemplo 2 (Observación) Considerando un tamaño de n = 10, 000 podemos construir el histograma de xi para ver si su distribución es similar a la de una normal. Figura: Histograma de las x′i s Del gráfico anterior puede verse que el histograma es simétrico respecto a la media 0 por lo que podrı́a sospecharse que las x′i s provienen de una distribución normal. Carlos V. Ramı́rez Ibáñez Universidad La Salle 42 / 44 Introducción Algunos Resultados Estimar una Medida Integración de Monte Carlo Muestreador de Gibbs Referencia Referencias Básicas Ross, Sheldon M., Simulation, Quinta Edición, Academic Press, EE.UU., 2013. Law, Averill, Simulation Modeling and Analysis, Quinta Edición, Mc Graw Hill, EE.UU., 2014. Fishman, G. A, First Course in Monte Carlo, Thomson Brooks, EE.UU., 2006. Anderson, D. R., Metodos Cuantitativos Para Los Negocios, Cengage Learning, México, 2016. Carlos V. Ramı́rez Ibáñez Universidad La Salle 43 / 44 Introducción Algunos Resultados Estimar una Medida Integración de Monte Carlo Muestreador de Gibbs Referencia Referencias Complementarias 1. Método Monte Carlo y Reducción de Varianza. Consultado en abril de 2021 en http://sistemas.fciencias.unam.mx/~silo/Cursos/simulacion/ tareas/cap6.pdf 2. Monte Carlo vı́a Cadenas de Markov. Consultado en febrero de 2021 en http://sistemas.fciencias.unam.mx/~silo/Cursos/simulacion/ tareas/cap7.pdf 3. El método Monte-Carlo y su aplicación a finanzas. Consultado en febrero de 2021 en http://mat.izt.uam.mx/mat/documentos/notas%20de% 20clase/cfenaoe3.pdf Carlos V. Ramı́rez Ibáñez Universidad La Salle 44 / 44