Intervalos de confianza bootstrap

Traballo Fin de Grao Intervalos de conanza bootstrap Alberto Portela Peleteiro 2020/2021 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Intervalos de conanza bootstrap Alberto Portela Peleteiro Julio 2021 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Trabajo propuesto Área de Coñecemento: Estatística e Investigación Operativa Título: Intervalos de conanza bootstrap Breve descrición do contido El método bootstrap permite aproximar la distribución de un estadístico, lo cual es de gran utilidad para llevar a cabo procedimientos de inferencia, como intervalos de conanza o contrastes de hipótesis. En este trabajo nos centraremos en la aplicación del bootstrap para la construcción de intervalos de conanza. Se llevará a cabo una revisión del método bootstrap en términos generales, así como de los métodos principales propuestos en la literatura para obtener intervalos de con- anza. Se estudiarán las propiedades teóricas de los distintos métodos y se realizarán estudios de simulación para comprobar su funcionamiento de datos simulados. Además, se ilustrarán con aplicaciones en datos reales. iii Índice general Resumen viii Introducción xi 1. Intervalos bootstrap percentil 1 1.1. Motivación y explicación del método . . . . . . . . . . . . . . . . . . . . . . 1 1.2. Comparación entre los intervalos bootstrap percentil y los intervalos normalesestándar.................................... 2 1.2.1. Cobertura de los intervalos . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.2. Restricciones de rango . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.3. Buen comportamiento frente a transformaciones . . . . . . . . . . . . 4 1.3. Ejemplos...................................... 4 1.4. SoftwaredeRasociado.............................. 10 1.4.1. Comandoboot .............................. 10 1.4.2. Comandoboot.ci............................. 11 2. Intervalos bootstrap-t 13 2.1. Motivación y explicación del método . . . . . . . . . . . . . . . . . . . . . . 13 2.2. Ventajas e inconvenientes del método bootstrap-t . . . . . . . . . . . . . . . 14 2.2.1. Malos resultados para muestras pequeñas . . . . . . . . . . . . . . . 15 2.2.2. Problemas computacionales derivados del bootstrap . . . . . . . . . 15 2.2.3. Cálculo de ˆσ cuando no tiene una expresión conocida . . . . . . . . . 15 2.2.4. Sensibilidad a cambios de escala . . . . . . . . . . . . . . . . . . . . 15 2.3. Relación entre el método bootstrap-t y el método bootstrap percentil . . . . 16 2.4. Estudiodesimulación .............................. 17 2.5. SoftwaredeRasociado.............................. 20 v vi ÍNDICE GENERAL 3. Intervalos BCa 21 3.1. Motivación y explicación del método . . . . . . . . . . . . . . . . . . . . . . 21 3.1.1. Parámetro z0 ............................... 22 3.1.2. Parámetro a ................................ 23 3.2. Propiedades de los intervalos BCa ....................... 24 3.2.1. Comportamiento frente a transformaciones . . . . . . . . . . . . . . . 24 3.2.2. Buenas propiedades de cobertura . . . . . . . . . . . . . . . . . . . . 25 3.3. Ejemplos...................................... 26 3.4. SoftwaredeRasociado.............................. 30 4. Intervalos ABC 31 4.1. Motivación y explicación del método . . . . . . . . . . . . . . . . . . . . . . 31 4.1.1. Coeciente ˆ b ............................... 33 4.1.2. Coeciente ˆcq ............................... 33 4.1.3. ABCcuadrático ............................. 33 4.2. Propiedades de los intervalos . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.3. Estudiodesimulación .............................. 35 4.4. SoftwaredeRasociado.............................. 37 Anexos 38 A. Código de R de las simulaciones 41 A.1. Código de R del primer estudio de simulación . . . . . . . . . . . . . . . . . 41 A.2. Código de R del segundo estudio de simulación . . . . . . . . . . . . . . . . 45 A.3. Código de R del tercer estudio de simulación . . . . . . . . . . . . . . . . . 51 Bibliografía 57 xiv INTRODUCCIÓN Sin embargo, al estudiar cómo se comportaba su método, Efron descubrió que era más exible que el propio método jackknife, pues para muestras de gran tamaño era mucho más ecaz. En este trabajo nos centraremos principalmente en el contexto no paramétrico del método bootstrap, que es más general. Veamos matemáticamente en qué consiste el método bootstrap en este caso: Como ya anticipamos previamente, el método bootstrap crea, a partir de una muestra inicial X , un número B sucientemente grande de submuestras del mismo tamaño que la muestra original. Estas submuestras se crean mediante la extracción con reemplazamiento de datos de la muestra original. Frente a las estimaciones teóricas de un parámetro, el bootstrap no asume que los datos que conforman la muestra proceden de una distribución particular, y solo utiliza la distribución empírica de los propios datos (que da un peso de 1 n a cada uno de ellos), por lo que es bastante menos restrictivo que los procedimientos teóricos de los que hablamos al principio de este texto. Suponiendo que, dado un conjunto de datos X , queremos calcular una estimación para la desviación típica de un estadístico T(X) empleando bootstrap, el procedimiento que deberíamos llevar a cabo es el siguiente: 1. Construir B muestras mediante extracciones con reemplazamiento que denotaremos por X∗(i) con i= 1, . . . , B . 2. Aplicar el estadístico T a cada uno de las B submuestras creadas, obteniendo T(X∗(i)) con i= 1, . . . , B . 3. Calcular la desviación típica muestral de estos B valores, y esa será la estimación de la desviación típica del estadístico T . Una de las principales ventajas del bootstrap frente a las estimaciones que emplean métodos teóricos es que, aunque la distribución bootstrap es, a menudo, difícil de calcular, se puede aproximar mediante simulación, por medio de algoritmos como el que se detalla arriba. Esto supone una gran mejoría con respecto a los métodos teóricos, pues no requiere de un tratamiento especíco según el conjunto de datos y el parámetro de interés. En este trabajo trataremos la construcción de intervalos de conanza empleando el método bootstrap, analizando 4 tipos de intervalo distintos, con su algoritmo de utilización, sus ventajas y desventajas, su aplicación computacional y ejemplos prácticos asociados, tanto con datos reales como simulados. INTRODUCCIÓN xv En el Capítulo 1 trataremos el primero de estos tipos de intervalos, conocidos como intervalos bootstrap percentil. En el Capítulo 2 trataremos los intervalos bootstrap-t o bootstrap estudentizados. En el Capítulo 3 trataremos los intervalos BCa y, nalmente, en el Capítulo 4 hablaremos de los intervalos ABC. xvi INTRODUCCIÓN Capítulo 1 Intervalos bootstrap percentil 1.1. Motivación y explicación del método Los intervalos bootstrap percentil son conceptualmente bastante sencillos, pero nos permitirán introducir otro tipo de intervalos más complejos como los BCa o los ABC , que veremos posteriormente. La idea principal en la que se basa el método es muy intuitiva a partir de la denición del concepto de bootstrap, y es obtener, mediante la muestra ordenada de las estimaciones bootstrap, los percentiles adecuados para calcular el intervalo de conanza que buscamos. Veámoslo matemáticamente: Sea X una muestra de tamaño n y θ nuestro parámetro de interés, cuya estimación muestral es ˆ θ . Construimos B submuestras bootstrap, todas de tamaño n , mediante extracciones con reemplazamiento y calculamos sus correspondientes estimaciones de θ que denotaremos por ˆ θ∗(i) , i= 1, . . . , B . Si denotamos por ˆ θ∗(α/2) al valor ordenado de los ˆ θ∗(i) que ocupa la posición B·α/2 , el intervalo percentil con un nivel de conanza (1 −α) será: (ˆ θ∗(α/2),ˆ θ∗(1−α/2)) Nótese que los valores ˆ θ∗(α/2) y ˆ θ∗(1−α/2) no son más que los cuantiles muestrales de orden α/2 y (1 −α/2) , respectivamente, de la muestra {ˆ θ∗(i)}B i=1 . Si, por ejemplo, α= 0.05 y B= 1200 , el extremo inferior del intervalo percentil asociado para un parámetro de interés θ será el valor ordenado de los ˆ θ∗(i) que ocupa la posición B·α/2 = 1200 ·0.025 = 30 , que se corresponde con el cuantil muestral de orden 0.025 , también llamado percentil del 2.5 % . 1 2 CAPÍTULO 1. INTERVALOS BOOTSTRAP PERCENTIL Observación 1.1 . Aquí surge un pequeño problema, pues si se emplea un valor de α demasiado preciso, o bien el número de muestras bootstrap B no es sucientemente redondo, puede ocurrir que B·α no sea un valor entero. Para solventar este inconveniente denimos el siguiente criterio: Dado k=b(B+ 1)α/2c , es decir, k es el mayor número natural menor o igual que (B+ 1)α/2 , tomamos ˆ θ∗(α/2) el valor ordenado de los valores ˆ θ∗(i) que ocupa la posición k y ˆ θ∗(1−α/2) el que ocupa la posición (B+ 1 −k) . 1.2. Comparación entre los intervalos bootstrap percentil y los intervalos normales estándar Introduzcamos en primer lugar el siguiente lema: Lema 1.2. Sea m:R−→ R una transformación continua y estrictamente monótona. Si ˆ φ=m(ˆ θ) normaliza perfectamente la distribución de ˆ θ , es decir, ˆ φ∼N(φ, σ2) , entonces el intervalo bootstrap percentil para θ es: m−1(ˆ φ−zα/2σ), m−1(ˆ φ+zα/2σ) La conclusión que podemos extraer de este resultado es que los intervalos bootstrap percentil tienen buen comportamiento frente a transformaciones. Para poder comparar las similitudes y las diferencias de los intervalos estándar y los intervalos bootstrap percentil analizaremos los dos casos siguientes: 1. ˆ θ∼N(θ, σ2) 2. m(ˆ θ)∼N(m(θ), σ2) En el primer caso en el que el estimador se aproxima a una distribución normal, los intervalos estándar y los bootstrap percentil producirán resultados similares, ambos bastante precisos. Sin embargo, en el segundo caso existe una notable mejoría de los intervalos bootstrap percentil frente a los intervalos estándar, pues, como hemos visto en el lema 1.2, el método percentil se comporta bien frente a transformaciones. Más concretamente, para que los intervalos estándar sean adecuados, se necesita conocer una transformación m:R−→ R tal que el estimador transformado m(ˆ θ) siga una distribución que se aproxima a una normal. Por el contrario, para que el método bootstrap percentil proporcione un intervalo correcto, basta con conocer que dicha transformación que normaliza las estimaciones existe, 1.2. INTERVALOS BOOTSTRAP PERCENTIL VS INTERVALOS ESTÁNDAR 3 pues, como hemos visto en el lema anterior, el intervalo bootstrap percentil será análogo al obtenido al realizar el intervalo estándar de las estimaciones normalizadas y aplicarle posteriormente la función inversa de la transformación m . El problema de estos dos tipos de intervalos es que, para casos más generales que los dos mencionados previamente, ninguno de los dos obtiene resultados óptimos. Sigamos ahora analizando ventajas y desventajas de los intervalos bootstrap percentil: Cobertura de los intervalos Restricciones de rango 1.2.1. Cobertura de los intervalos Denamos en primer lugar qué entendemos por cobertura de un método para el cálculo de intervalos de conanza. Denición 1.3. La cobertura de un método para el cálculo de intervalos de conanza es la probabilidad efectiva de que dicho intervalo contenga al parámetro en cuestión cuando se construye con un nivel de conanza (1−α) . La cobertura se suele aproximar en estudios de simulación mediante el porcentaje de intervalos que contienen al parámetro cuando se generan muchas muestras. Aunque por lo general los intervalos de conanza bootstrap percentil son más equilibrados (en el sentido de que la probabilidad de que el parámetro real sea inferior al extremo superior al intervalo es similar a la probabilidad de que sea mayor que el extremo superior) que los intervalos estándar, la cobertura de este tipo de métodos no suele ser tan precisa como nos gustaría. Estos problemas de cobertura serán solucionados posteriormente al introducir el método bootstrap-t. 1.2.2. Restricciones de rango Consideremos que el parámetro de interés, θ , tiene una restricción de rango, es decir, toma valores en un subintervalo K⊂R y no tienen sentido valores que no pertenezcan a dicho intervalo. Este tipo de restricciones de rango son muy comunes pues, por ejemplo, si el parámetro de interés, θ , es una probabilidad, tendrá que tomar valores pertenecientes al intervalo [0,1] , y si el parámetro θ es un coeciente de correlación, necesariamente tomará valores en el intervalo [−1,1] . 4 CAPÍTULO 1. INTERVALOS BOOTSTRAP PERCENTIL Las restricciones de rango suponen un problema para varios tipos de intervalos de conanza, ya que, como es evidente, lo que se busca es que los extremos de un intervalo de conanza también respeten las restricciones de rango que tenemos para el parámetro. Ésta es una de las importantes ventajas de los métodos bootstrap percentil, puesto que, por construcción, el intervalo de conanza para el parámetro buscado respetará la restricción de rango. Esto se tiene porque, para construir el intervalo, estamos empleando las estimaciones bootstrap, que respetan la restricción de rango, y considerando sus percentiles, que trivialmente también conservan la restricción. Sin embargo, los intervalos estándar no verican, en general, esta condición, por lo que puede ocurrir que alguno de los extremos del intervalo no pertenezca al rango de valores en los que el parámetro tiene sentido, y esto haga que el intervalo de conanza no proporcione información válida. Este tipo de intervalos que preservan el rango tienen, en general, buenas propiedades y son más precisos y ecaces que los que no lo preservan. 1.2.3. Buen comportamiento frente a transformaciones El método boostrap percentil se comporta bien frente a transformaciones, en el sentido de que si (ˆ θ1,ˆ θ2) es el intervalo bootstrap percentil para θ de nivel (1 −α) y m es una aplicación monótona creciente que transforma el propio parámetro, entonces el intervalo bootstrap percentil de nivel (1 −α) para m(θ) es (m(ˆ θ1), m(ˆ θ2)) . 1.3. Ejemplos Para ilustrar la comparación entre los intervalos bootstrap percentil y los intervalos estándar teóricos, que provienen de la suposición de que el estimador sigue una distribución normal, detallaremos los dos ejemplos siguientes: Ejemplo 1.4. En este ejemplo utilizaremos la base de datos rock del paquete datasets de R y analizaremos una muestra de 48 rocas de una reserva petrolera. El objetivo es construir un intervalo de conanza para la media de las áreas de los poros de cada roca, medidas en píxeles. Estas áreas se distribuyen según se puede ver en el histograma de la gura 1.1. A la vista del histograma, los datos se distribuyen de forma acampanada por lo que no parece descabellada la idea de que las áreas de los poros de las rocas sigan una distribución normal. 1.3. EJEMPLOS 5 Figura 1.1: Histograma de la muestra de datos rock de R. Para comprobar nuestra primera intuición estudiaremos brevemente la posible normalidad de la muestra empleando un test de normalidad de Shapiro-Wilk y una gráca cuantil-cuantil (también llamada qqplot). shapiro.test(x) ## ## Shapiro-Wilk normality test ## ## data: x ## W = 0.97944, p-value = 0.5555 El test de Shapiro-Wilk nos devuelve un nivel crítico alto, del 56 % , por lo podemos suponer que los datos de la muestra provienen de una distribución normal para niveles de signicación usuales ( 1 %,5 % y 10 % ). Veámoslo también grácamente con la gráca qq-plot de la gura 1.2, que nos dice que si los datos se distribuyen sobre la recta que dene una normal, podemos aceptar que los datos de la muestra provienen de una distribución de este tipo. Nuevamente, como los datos están en su mayoría sobre la recta, podemos aceptar que provienen de una normal. Estamos entonces en una situación similar a la descrita en la introducción del trabajo, donde se conocen intervalos exactos para parámetros especícos como la media o la 6 CAPÍTULO 1. INTERVALOS BOOTSTRAP PERCENTIL Figura 1.2: Gráca cuantil-cuantil sobre los datos rock . varianza. Veremos de esta forma cómo se comportan los intervalos bootstrap percentil frente a los intervalos de la T de Student, en un contexto que, a priori, es favorable para estos últimos, pues la suposición de normalidad parece cumplirse. Recordemos que se quiere estimar el valor medio del área de las rocas de la reserva petrolera, dando un intervalo de conanza del 95 % para él. Construyamos primero el intervalo de conanza de la T de Student. Para ello, suponemos que los datos provienen de una distribución normal de la que no conocemos su desviación típica. Como detallamos en la introducción, el intervalo de conanza de la T de Student para un nivel de conanza de (1 −α) viene dado por: ¯ X−tα/2 Sc √n,¯ X+tα/2 Sc √n Para los datos de nuestra muestra y con nivel 1−α= 0.95 obtenemos el intervalo (6408,7967) , mientras que la media muestral es 7188 . Hemos obtenido por simulación estimaciones bootstrap para la media, ¯ X∗(i), i = 1, . . . , B , con B muestras bootstrap simuladas. El histograma correspondiente a estas B réplicas de la media muestral se tiene en la gura 1.3. Observamos un comportamiento que se puede aproximar al de una normal y por lo tanto estaríamos en el primer caso de los detallados en la sección 1.2 de la página 2 donde también el intervalo bootstrap percentil tiene propiedades aceptables. 1.3. EJEMPLOS 7 Figura 1.3: Histograma de las estimaciones bootstrap de la media sobre los datos rock . Las líneas verticales que aparecen en el histograma son los cuantiles de nivel α/2=0.025 y 1−α/2=0.975 , que determinan el intervalo de conanza bootstrap percentil para la media. De esta manera el intervalo es (6438,7926) . Denamos ahora un concepto sencillo que nos permitirá comparar estos dos intervalos. Denición 1.5. Denimos la forma de un intervalo de conanza para θ , [θinf , θsup] , como: θsup −ˆ θ ˆ θ−θinf donde ˆ θ es la estimación muestral de θ . La forma de un intervalo de conanza mide la simetría del intervalo con respecto a ˆ θ . Los intervalos con forma igual a 1 tienen una simetría total con respecto a ˆ θ . Los intervalos con forma menor que 1 presentan más longitud a la izquierda de ˆ θ que a la derecha y, por el contrario, los intervalos con forma mayor que 1 presentan más longitud a la derecha de ˆ θ . Para poder comparar estos intervalos que hemos obtenido, elaboraremos una tabla que contenga los siguientes elementos: Extremos del intervalo Longitud del intervalo Forma del intervalo 14 CAPÍTULO 2. INTERVALOS BOOTSTRAP-T Cada conjunto de datos X∗ nos proporciona una estimación ˆ θ∗ de nuestro parámetro de interés y una estimación de la desviación típica ˆσ∗ . Así, construímos los siguientes estadísticos T∗ : T∗=ˆ θ∗−ˆ θ ˆσ∗ Ahora tenemos B valores para T∗ , uno para cada muestra bootstrap, y, tras haberlos ordenado de menor a mayor, denimos ˆ T(α) como el T∗ que ocupa la posición B·α . Observación 2.1 . Al igual que para los intervalos bootstrap percentil, puede ocurrir que una mala elección de B o de α propicie que B·α no sea entero. Para solucionar este inconveniente procederíamos igual que en la observación 1.1 Arreglado este problema técnico, para un nivel de conanza (1−α) tenemos el intervalo siguiente: (ˆ θ−ˆ T(1−α/2)ˆσ, ˆ θ−ˆ T(α/2)ˆσ) 2.2. Ventajas e inconvenientes del método bootstrap-t Al contrario que los intervalos normales estándar, los intervalos bootstrap-t no son, en general, simétricos. Este hecho permite que los intervalos bootstrap-t puedan ser más precisos que los intervalos estándar, sobre todo en muestras que presenten un alto grado de asimetría o muestras con un gran número de datos. El método bootstrap-t trabaja mejor si θ es una medida de localización, es decir, si θ es un parámetro que verica que, al aumentar en una constante k todos los datos de la muestra, θ también aumentará en k. Ejemplos de medidas de localización son la media, la mediana o la moda. Además, generalmente obtiene resultados ligeramente menos precisos que los métodos BCa y ABC que veremos posteriormente y que son algo más complejos. Sin embargo, pese a que el método bootstrap-t parece mejorar la precisión de los intervalos estándar en bastantes casos, también presenta unos notables inconvenientes en su ejecución: Malos resultados para muestras pequeñas. Problemas propios del bootstrap como el nivel de computación. Cálculo de la estimación ˆσ cuando no tiene una expresión conocida. Sensibilidad a cambios de escala. 2.2. VENTAJAS E INCONVENIENTES DEL MÉTODO BOOTSTRAP-T 15 Nos centraremos individualmente en cada uno de estos problemas y cómo poder solucionarlos si es posible. 2.2.1. Malos resultados para muestras pequeñas Como ya hemos comentado previamente, el método bootstrap-t trabaja mejor cuando tenemos una muestra sucientemente grande, para que tenga sentido crear réplicas bootstrap que proporcionen información relevante sobre la muestra inicial. 2.2.2. Problemas computacionales derivados del bootstrap Si el cálculo de ˆ θ o de ˆσ es computacionalmente costoso, como hay que realizarlo un número alto B de veces, pueden surgir problemas con el tiempo de cálculo computacional de los intervalos. 2.2.3. Cálculo de ˆσ cuando no tiene una expresión conocida En el cálculo de los valores T∗ aparece ˆσ en el denominador. Para parámetros θ comunes, como la media, existen expresiones conocidas y sencillas para la obtención de la desviación típica. Sin embargo, cuando θ es un parámetro más general, lo más habitual es que no exista una expresión conocida. Es común en estos casos emplear técnicas bootstrap para estimar estas desviaciones típicas. Aunque lo habitual es que un número Bσ= 50 de réplicas bootstrap proporcionen una estimación de ˆσ sucientemente buena, es necesario realizar estas réplicas para cada una de las B muestras que necesitábamos para calcular los percentiles. Así, si tomamos B= 1000 , que es una cantidad estándar de réplicas bootstrap, necesitaríamos B·Bσ= 1000 ·50 = 50000 muestras, que es un número muy considerable. Por este hecho decimos que el método bootstrap-t trabaja mejor cuando el parámetro de interés, θ , es una medida de traslación, pues para este tipo de parámetros es común conocer expresiones para la desviación típica. 2.2.4. Sensibilidad a cambios de escala Al contrario que los métodos BCa y ABC , el bootstrap-t es muy sensible a cambios de escala, es decir, los intervalos de conanza dieren bastante si la muestra está sometida a una transformación o no. Este hecho provoca que una mala elección de escala genere un intervalo de conanza impreciso. Para algunos parámetros, como el coeciente de correlación, es conocido el tipo de transformación adecuada para obtener un intervalo de conanza preciso, pero en general 16 CAPÍTULO 2. INTERVALOS BOOTSTRAP-T la transformación correcta es desconocida, lo que provoca una barrera a priori insalvable para este método, pues, como hemos dicho anteriormente, una mala elección de la escala puede tener graves consecuencias. Para poder solventar este inconveniente se puede emplear bootstrap y algunos conceptos algo más avanzados para conseguir determinar la trasformación adecuada para nuestro parámetro de interés θ y nuestro conjunto de datos. 2.3. Relación entre el método bootstrap-t y el método bootstrap percentil Como hemos visto, el pivote a partir del cual construimos el intervalo bootstrap-t es: T=ˆ θ−θ ˆσ Si suponemos que ˆσ= 1 , el pivote sería T=ˆ θ−θ . Construyendo ahora el intervalo bootstrap-t obtenemos: (ˆ θ−[ˆ θ∗−ˆ θ](1−α/2),ˆ θ−[ˆ θ∗−ˆ θ](α/2)) Y tras una simple modicación tenemos: (2ˆ θ−ˆ θ∗(1−α/2),2ˆ θ−ˆ θ∗(α/2)) (2.1) Podemos compararlo con el intervalo bootstrap percentil para el mismo parámetro de interés θ : (ˆ θ∗(α/2),ˆ θ∗(1−α/2)) (2.2) A la vista de las expresiones 2.1 y 2.2 observamos que si uno de los intervalos presenta una asimetría negativa (positiva), el otro presentará una asimetría positiva (negativa). Esto nos puede llevar a pensar que al menos uno de los dos tipos de intervalos debe ser erróneo, pues parece que ambos llegan a resultados contradictorios. Sin embargo, dependiendo de la situación, será más ecaz un tipo de intervalo que otro, es decir, ninguno de ellos es incorrecto. Veamos ahora estas y otras propiedades del método bootstrap-t con un estudio de simulación. 2.4. ESTUDIO DE SIMULACIÓN 17 2.4. Estudio de simulación En esta sección simularemos 1000 muestras de tamaño n= 50 de datos procedentes de distintas distribuciones y calcularemos intervalos de conanza para la media empleando tres métodos: T de Student, bootstrap percentil y bootstrap-t, en los tres casos con un nivel de conanza del 95 % . La idea será estudiar el comportamiento de los distintos tipos de intervalos según la procedencia del conjunto de muestras con el que trabajemos. El lenguaje de R asociado a este ejemplo se expondrá en el anexo para el caso de que las muestras procedan de una distribución normal estándar y será análogo para el resto de tipos de muestra. La tabla 2.1 contiene los valores teóricos de las medias de las distribuciones con las que hemos trabajado anteriormente vienen dadas por la siguiente tabla: N(0,1) N(0,42) Exp(1) Γ(2,2) B( 1 2,1 2 ) B(5,1) Media teórica 0 0 1 1 1 2 5 6 Tabla 2.1: Medias teóricas de cada tipo de distribución: Normal estándar (N(0,1)), Normal de media 0 y desviación típica 4 ( N(0,42) ), Exponencial de parámetro 1 (Exp(1)), Gamma con parámetros de forma y escala 2 ( Γ(2,2) ), Beta con parámetros 1/2 y 1/2 (B( 1 2,1 2 )) y Beta con parámetros 5 y 1 (B(5,1)). Como generamos 1000 muestras, obtenemos otros tantos intervalos de cada tipo. Para dar una idea de cómo será el intervalo promedio de cada tipo de método, trataremos en cada caso el extremo inferior medio y el extremo superior medio de cada uno de ellos.Recogemos los resultados obtenidos en la tabla 2.2. N(0,1) N(0,42) Exp(1) Γ(2,2) B( 1 2,1 2 ) B(5,1) T de Stud. (-0.29,0.27) (-1.17,1.09) (0.72,1.28) (0.80,1.20) (0.40,0.60) (0.79,0.87) Percentil (-0.28,0.26) (-1.12,1.05) (0.75,1.28) (0.82,1.20) (0.41,0.60) (0.79,0.87) Boot-t (-0.29,0.27) (-1.17,1.09) (0.76,1.35) (0.82,1.24) (0.40,0.60) (0.79,0.87) Tabla 2.2: Intervalos promedio que proporcionan los distintos métodos. Observamos que, como se podía prever, todos los intervalos contienen a la media teórica de cada distribución . A la vista de la tabla 2.2, el método bootstrap-t produce intervalos similares a los de la T de Student en poblaciones normales. En otras distribuciones como 18 CAPÍTULO 2. INTERVALOS BOOTSTRAP-T la exponencial o la gamma, los intervalos de la T de Student y los del bootstrap-t son ligeramente diferentes. Para valorar cuál de ellos es más preciso estudiaremos luego sus correspondientes coberturas. En vista del comportamiento frente a poblaciones normales, donde el método bootstrapt obtiene un intervalo similar al obtenido con el estadístico de la T de Student y ligeramente distinto al bootstrap percentil, parece que el método bootstrap-t reeja una cierta mejoría en el trato de poblaciones normales con respecto al bootstrap percentil, pues se asemeja más al intervalo exacto de la T de Student. Para poder analizar cómo se comportan los distintos métodos en los casos en que la muestra no procede de una distribución normal, emplearemos las tablas 2.3, 2.4 y 2.5, que contienen las longitudes medias, formas promedio y coberturas de los tres métodos. N(0,1) N(0,42) Exp(1) Γ(2,2) B( 1 2,1 2 ) B(5,1) T de Student 0.565 2.262 0.557 0.400 0.201 0.079 Bootstrap percentil 0.543 2.174 0.534 0.383 0.193 0.076 Bootstrap-t 0.564 2.257 0.592 0.412 0.201 0.081 Tabla 2.3: Longitudes de los intervalos promedio obtenidos por los distintos métodos. N(0,1) N(0,42) Exp(1) Γ(2,2) B( 1 2,1 2 ) B(5,1) T de Student 1 1 1 1 1 1 Bootstrap percentil 1.003 1.003 1.12 1.08 1.002 0.93 Bootstrap-t 1.005 1.005 1.47 1.32 1 0.79 Tabla 2.4: Forma promedio de los distintos tipos de intervalos. N(0,1) N(0,42) Exp(1) Γ(2,2) B( 1 2,1 2 ) B(5,1) T de Student 0.953 0.953 0.939 0.951 0.948 0.946 Bootstrap percentil 0.944 0.944 0.933 0.941 0.944 0.938 Bootstrap-t 0.954 0.954 0.952 0.957 0.957 0.948 Tabla 2.5: Coberturas de los intervalos obtenidos por los distintos métodos. De la tabla 2.3, lo que más llama la atención es que el método bootstrap percentil devuelve intervalos con una longitud considerablemente inferior a los otros dos, y esto 2.4. ESTUDIO DE SIMULACIÓN 19 se ve reejado luego en la tabla 2.5 en una cobertura claramente peor que los otros dos métodos. Por otro lado, para los datos procedentes de distribuciones no normales, el método bootstrap-t obtiene intervalos ligeramente más amplios que el intervalo de la T de Student. Sin embargo, para las poblaciones normales, los intervalos de la T de Student son ligeramente más amplios. Esto, que se ve reejado en la cobertura de los intervalos que, aunque en general un intervalo más amplio no tiene por qué vericar una mejor cobertura, en estos casos sí que se tiene. En los casos que hemos tratado, los métodos producen en general intervalos bastante simétricos (con forma próxima a 1), aunque hay algunas excepciones, como los intervalos dados por el método bootstrap-t para la distribución exponencial, la gamma, y la B(5,1) , lo cual se debe a la asimetría de dichas distribuciones en torno a su media teórica. Este buen comportamiento del método bootstrap frente a distribuciones asimétricas es una de las principales ventajas de los métodos bootstrap y del método bootstrap-t en particular. Nos centraremos nalmente en la tabla de las coberturas, que es la más importante de todas pues nos da una idea de si los intervalos que estamos construyendo verican verdaderamente el nivel de conanza nominal del 95 % para el que se buscaba construir el intervalo. Aquí se observa un comportamiento excelente del bootstrap-t, pues obtiene mejor cobertura de los intervalos que los otros dos métodos en las seis distribuciones que estamos considerando. Además en cinco de ellos obtiene un nivel de cobertura superior a 0.95 , lo que indica que los intervalos calculados son realmente intervalos para la media con un nivel de conanza del 95 % . Nótese además que para el caso de la distribución B(5,1) en el que no se obtiene el umbral del 0.95 , el valor obtenido ( 0.948 ) es muy cercano a dicho valor. Por otro lado se observa bien la pobre cobertura de los intervalos bootstrap percentil y los intervalos de la T de Student para algunas distribuciones, en concreto para la distribución exponencial. Este hecho conrma lo que ya se comentó en el apartado 1.2.1 de que los intervalos bootstrap percentil presentan ciertos problemas de cobertura. En resumen, en base a esta simulación y a la explicación teórica, determinamos que el método bootstrap-t perfecciona la idea del método bootstrap percentil y obtiene, en líneas generales, mejores resultados. Sin embargo, el método bootstrap-t para la construcción de los intervalos de conanza también tiene defectos. Para corregir estos pequeños defectos de este método, surgieron 20 CAPÍTULO 2. INTERVALOS BOOTSTRAP-T nuevos tipos de intervalos que consiguen mejorar alguna propiedad del método bootstrap-t como su mal comportamiento frente a transformaciones. 2.5. Software de R asociado Los comandos de R que se pueden emplear para la construcción de intervalos de con- anza por el método bootstrap-t son completamente análogos a los empleados para el caso de los bootstrap percentil, con la salvedad de que en el argumento type deberemos poner stud. Capítulo 3 Intervalos BCa En este capítulo nos centraremos en el primero de los métodos que consigue mejorar considerablemente las propiedades de los procedimientos estudiados hasta ahora: el método BCa . 3.1. Motivación y explicación del método El método BCa surge con la idea de conseguir reunir las ventajas de los métodos bootstrap-t y bootstrap percentil en un solo método. Su nombre es una abreviatura de bias-corrected and accelerated y es básicamente un perfeccionamiento del método bootstrap percentil. El método BCa también construye intervalos de conanza basados en percentiles del histograma de las réplicas bootstrap. Sin embargo, la gran diferencia con los intervalos percentil es que no toma los valores α/2 y 1−α/2 , sino que realiza una ligera modicación que le permite ajustarse más al conjunto de datos. Como veremos posteriormente, este método supone una ligera mejoría frente a los dos métodos bootstrap que hemos estudiado hasta ahora, aunque también presenta ciertos inconvenientes como una mala precisión ocasional cuando la muestra es de pequeño tamaño y un alto coste computacional, pues precisa bastantes más operaciones que el propio método bootstrap percentil. Recordemos que el intervalo de conanza bootstrap percentil para un parámetro θ con un nivel de conanza de (1 −α) viene dado por: (ˆ θ∗(α/2),ˆ θ∗(1−α/2)) El intervalo BCa no se aleja de esta idea de intervalo. De hecho, el intervalo de este método viene dado por: 21 22 CAPÍTULO 3. INTERVALOS BCA (ˆ θ∗(α1),ˆ θ∗(α2)) La diferencia con respecto al bootstrap percentil radica únicamente en los valores α1 y α2 que, aunque no son excesivamente difíciles de calcular, sí que requieren un cierto cálculo computacional. Estos valores α1 y α2 de los que estamos hablando se calculan como sigue: α1= Φ z0+z0+z(α/2) 1−a(z0+z(α/2))! α2= Φ z0+z0+z(1−α/2) 1−a(z0+z(1−α/2))! donde z0 es un parámetro que corrige el sesgo y a es un parámetro llamado aceleración en el que profundizaremos posteriormente. Además, Φ es la función de distribución de una normal estándar y z(α/2) y z(1−α/2) son los percentiles 100 ·α/2 y 100 ·(1 −α/2) de esta distribución. Estos valores a partir de los cuales se calculan los percentiles parecen tener una expresión muy compleja. Veamos que no lo es tanto: En primer lugar, si tomamos la mayor simplicación posible, es decir, a=z0= 0 , entonces: α1= Φ(z(α/2)) = α/2 α2= Φ(z(1−α/2)) = 1 −α/2 Es decir, para el caso particular en que a= 0 y z0= 0 el intervalo BCa coincide con el intervalo bootstrap percentil visto en el primer capítulo. Esto nos permitirá, según qué casos, mejorar el intervalo dado por el método bootstrap percentil mediante la introducción de un parámetro de corrección de sesgo z0 y un parámetro de aceleración a que mida la velocidad de cambio del error típico. Observación 3.1 . Efron y Tibshirany(1986) arman que un número B= 1000 réplicas bootstrap es, en general, suciente para obtener un intervalo de conanza sucientemente bueno empleando el método BCa . 3.1.1. Parámetro z0 Como se ha dicho anteriormente, el parámetro z0 es una medida de correción del sesgo. Su estimación más conocida, que denotaremos por ˆz0 , viene dada por: 3.1. MOTIVACIÓN Y EXPLICACIÓN DEL MÉTODO 23 ˆz0= Φ−1 PB b=1,ˆ θ∗(b)<ˆ θ1 B  (3.1) donde B es el número de réplicas bootstrap, ˆ θ es el estimador del parámetro θ , ˆ θ∗(b), b = 1, . . . , B son las réplicas bootstrap del estimador y Φ−1 es la inversa de la función de distribución normal estándar. Este parámetro, que se extrae directamente de las réplicas bootstrap, mide la diferencia entre ˆ θ y ˆ θ∗ a través de la proporción de réplicas ˆ θ∗ menores que el estimador ˆ θ . Tiene una interpretación muy sencilla pues ˆz0= 0 cuando la mitad de los valores ˆ θ∗(b) son menores que ˆ θ . El parámetro z0 , como hemos dicho, busca corregir el sesgo, dando sentido al nombre de la primera parte del método ( bias-corrected ). Veamos ahora el segundo parámetro principal, que requiere un cálculo más complejo. 3.1.2. Parámetro a El parámetro a , también conocido como parámetro de aceleración, tiene como objetivo analizar lo rápido que cambia el error típico de ˆ θ con respecto al verdadero valor del parámetro de interés θ . Hay distintas aproximaciones del valor de a en la literatura. La que emplearemos en este trabajo se puede ver en Efron(1987) y evita aproximaciones complejas en el cálculo de derivadas. Aunque probablemente no sea la aproximación más precisa, presenta buenas propiedades en la práctica y reduce considerablemente el cálculo computacional. Además, no deja de lado el objetivo de los métodos bootstrap, que es realizar inferencia apoyándose en las submuestras creadas. La idea de esta aproximación se basa en el método jackknife, que introdujimos previamente. Como breve recordatorio, dada una muestra de tamaño n, el método jackknife hace n estimaciones del parámetro de interés θ , suprimiendo en cada una de ellas un dato diferente de la muestra inicial. Sea X una muestra de datos de tamaño n y ˆ θ el estimador de θ obtenido con la muestra X. Para poder introducir la función de inuencia jackknife, que nos ayudará en la estimación de a , denotaremos ˆ θ(i) al estimador de θ obtenido con la muestra X eliminando el dato i-ésimo. A su vez, consideraremos la media de las estimaciones ˆ θ(i) , que denotamos: ˆ θ(·)= n X i=1 ˆ θ(i) n Estamos ahora en condiciones de denir la función de inuencia jackknife, que nos permitirá escribir de forma sencilla la aproximación del parámetro a . La función de inuencia 30 CAPÍTULO 3. INTERVALOS BCA implementar y requieren de un alto coste computacional. Para el caso que nos ocupa, las aproximaciones ˆa y ˆz0 son sucientemente buenas, pues logran mejorar considerablemente la cobertura del método bootstrap percentil y el comportamiento frente a transformaciones del método bootstrap-t. Hemos observado así con estos dos ejemplos el motivo de la creación del método BCa , que consigue reunir en cierto modo las ventajas de los dos métodos anteriores. 3.4. Software de R asociado Para construir los intervalos de tipo BCa en R, emplearemos el paquete boot que hemos tratado en los dos capítulos anteriores, con la única modicación de cambiar el valor del argumento type del comando boot.ci por bca. Capítulo 4 Intervalos ABC 4.1. Motivación y explicación del método El método BCa que hemos analizado en la sección anterior tiene un notable problema a la hora de aplicarlo en la práctica, y es su alto coste computacional. Para sobrellevar este inconveniente surge el método ABC para la creación de intervalos de conanza. El nombre de este método es una abreviatura de approximate bootstrap condence intervals y busca, como se puede extraer de su nombre, aproximar los extremos del intervalo BCa de forma analítica. El hecho de recurrir al mundo analítico para construir intervalos de conanza hace que algunos autores no lo consideren un método bootstrap, pues sustituye en gran medida el cálculo computacional derivado de las réplicas bootstrap por aproximaciones empleando cálculos de derivadas. Sin embargo, en este trabajo sí que lo consideraremos un método bootstrap, pues no es más que una reformulación de uno de los métodos vistos, empleando herramientas propias del análisis, que veremos posteriormente. Este método solventa en gran medida el alto coste computacional que requería el método BCa , pues en algunos casos logra reducirlo en más de un 90 % . Veamos ahora como construir los extremos del intervalo. Sea X una muestra de datos de tamaño n. Si el estadístico de interés es una función T de la muestra, se puede probar (Efron y Tibshirani, 1993) que la aproximación de la desviación típica de ˆ θ=T(X) empleando el método delta es: ˆσ= n X i=1 ˙ T2 i n2!1 2 donde ˙ Ti es la conocida como función de inuencia de primer orden y que se dene como sigue: 31 32 CAPÍTULO 4. INTERVALOS ABC ˙ Ti= l´ım →0 T((1 −)P0+~ei)−T(P0)  siendo ~ei el vector i-ésimo de la base canónica. Si (ˆ θABC[α/2],ˆ θABC[1−α/2]) es el intervalo de conanza para θ con nivel (1−α) dado por el método ABC, sus extremos se calculan como sigue: w= ˆz0+z(α/2) λ=w (1 −ˆaw)2 ˆ δ=˙ T(P0) ˆ θABC[α/2] = T P0+λˆ δ ˆσ! A ˆ δ se le conoce como dirección menos favorable. La principal ventaja del método ABC frente al método BCa es que las estimaciones ˆz0 y ˆa se pueden expresar como derivadas de primer y segundo orden, evitando el elevado cálculo computacional que se tenía en el BCa propiciado por el remuestreo. Este es el hecho por el que, como decíamos anteriormente, algunos autores no consideran este método como un método bootstrap. La estimación ˆa es muy similar a la que realizábamos en el capítulo anterior, sustituyendo en este caso la función de inuencia jackknife, U , por la función de inuencia de primer orden, ˙ T , denida previamente. ˆa=Pn i=1 ˙ T3 i 6Pn i=1 ˙ T2 i3 2 (4.1) La estimación de z0 no es tan sencilla como la anterior y requiere de tres coecientes: La estimación ˆa . Una aproximación del sesgo, ˆ b . Un coeciente de no-linelidad, que denotamos por ˆcq . De la aproximación ˆa ya nos hemos ocupado anteriormente. Nos centraremos ahora en los parámetros ˆ b y ˆcq . 4.1. MOTIVACIÓN Y EXPLICACIÓN DEL MÉTODO 33 4.1.1. Coeciente ˆ b Este valor es una aproximación del sesgo (b=E(ˆ θ)−θ) y viene dado por el desarrollo en serie de Taylor de θ=T(P0) , ˆ b= n X i=1 ¨ Ti 2n2 (4.2) donde ¨ Ti= l´ım →0 T((1 −)P0+~ei)−2T(P0) + T((1 −)P0−~ei) 2 es la componente i-ésima de la función de inuencia de segundo orden. 4.1.2. Coeciente ˆcq El valor ˆcq , conocido como coeciente cuadrático, mide la no-linealidad de la función θ=T(P) , cuando nos desplazamos en la dirección de ˆ δ . Su expresión analítica es: ˆcq= l´ım →0 T(1 −)P0+˙ T n2ˆσ−2T(P0) + T(1 −)P0−˙ T n2 2 (4.3) Ahora, con las expresiones (4.2) y (4.3) denimos: ˆγ=ˆ b ˆσ−ˆcq (4.4) Estamos nalmente en condiciones de dar la estimación ˆz0 , que viene dada por: ˆz0=φ−1{2φ(ˆa)·φ(−γ)}= ˆa−ˆγ Esta forma del método ABC es la más sencilla que existe, aunque no soluciona del todo el problema computacional pues requiere de la evaluación de las funciones T y ˙ T para el cálculo de ˆ δ y del propio extremo del intervalo. 4.1.3. ABC cuadrático Para solventar este inconveniente nace el ABC cuadrático o ABCq que emplea las estimaciones ˆa , ˆz0 y ˆcq para denir: w= ˆz0+z(α/2) λ=w (1 −ˆaw)2 34 CAPÍTULO 4. INTERVALOS ABC ξ=λ+ ˆcqλ2 ˆ θABC[α/2] = ˆ θ+ ˆσξ Sin embargo, esta versión del método ABC tampoco es perfecta, pues el ABCq es un método local. Esto hace que en algunos contextos surjan problemas computacionales derivados de trabajar fuera del conjunto donde el método produce resultados coherentes. 4.2. Propiedades de los intervalos De manera análoga al BCa , el método ABC produce intervalos con precisión de segundo orden y que respetan transformaciones monótonas. De hecho cuando B crece, los intervalos ABC y BCa tienden a ser el mismo. Denamos ahora un par de conceptos que nos ayudarán a relacionar los tipos de intervalo estudiados hasta el momento. Denición 4.1 (Exactitud de primer orden) . Sea (θ1, θ2) un intervalo de conanza exacto para θ con nivel de conanza (1 −α) . Un intervalo de conanza (ˆ θ1,ˆ θ2) se dice que tiene exactitud de primer orden si: ˆ θ1=θ1+k1 n ˆ θ2=θ2+k2 n donde k1 y k2 son constantes que no dependen de n. Denición 4.2 (Exactitud de segundo orden) . Sea (θ1, θ2) un intervalo de conanza exacto para θ con nivel de conanza (1 −α) . Un intervalo de conanza (ˆ θ1,ˆ θ2) se dice que tiene exactitud de segundo orden si: ˆ θ1=θ1+k1 √n3 ˆ θ2=θ2+k2 √n3 donde k1 y k2 son constantes, al igual que en la denición anterior. A la vista de las deniciones anteriores, se puede llegar a la siguiente proposición: Proposición 4.3. Si un intervalo (θ1, θ2) tiene exactitud de orden m, entonces también tiene precisión de orden m. 4.3. ESTUDIO DE SIMULACIÓN 35 Esta proposición nos permite comparar los distintos tipos de intervalos. Como ya vimos en el capítulo anterior, salvo en casos muy particulares, los intervalos normales asintóticos, T de Student y bootstrap percentil tienen precisión de primer orden, por lo que, a lo sumo, tendrán exactitud de primer orden. Esto supone una notable diferencia tanto con los intervalos BCa como con su aproximación empleando el método ABC, ya que, en la gran mayoría de los casos, estos métodos obtienen intervalos con exactitud de segundo orden. Veamos ahora si el método ABC ofrece una buena aproximación de los intervalos BCa con un nuevo estudio de simulación. 4.3. Estudio de simulación En este estudio generaremos, al igual que en la sección 2.4, 1000 muestras de tamaño n= 50 de seis distribuciones distintas, que serán las mismas que empleamos en dicho estudio. El objetivo será comparar los intervalos promedio dados por el método BCa con los correspondientes del método ABC , analizando también sus longitudes, coberturas y formas promedio. De manera análoga a los estudios de simulación anteriores, el código en lenguaje R asociado se puede ver en el anexo. En primer lugar, recordemos que las medias teóricas de las seis distribuciones vienen dadas por la tabla 4.1. N(0,1) N(0,42) Exp(1) Γ(2,2) B( 1 2,1 2 ) B(5,1) Media teórica 0 0 1 1 1 2 5 6 Tabla 4.1: Medias teóricas de cada tipo de distribución: Normal estándar (N(0,1)), Normal de media 0 y desviación típica 4 ( N(0,42) ), Exponencial de parámetro 1 (Exp(1)), Gamma con parámetros de forma y escala 2 ( Γ(2,2) ), Beta con parámetros 1/2 y 1/2 (B( 1 2,1 2 )) y Beta con parámetros 5 y 1 (B(5,1)). Lo primero que estudiaremos serán los intervalos promedio que obtiene cada método y que viene recogido en la tabla 4.2. Observamos que los intervalos son prácticamente idénticos con precisión de dos decimales para las seis distribuciones. Esto nos da una idea de que, a priori, el método ABC devuelve intervalos muy similares a los obtenidos mediante el método BCa . Sin embargo, aunque los intervalos promedio coincidan, podría ocurrir que los intervalos ABC presentasen problemas de cobertura o longitudes y asimetrías muy distintas a las obtenidas 36 CAPÍTULO 4. INTERVALOS ABC N(0,1) N(0,42) Exp(1) Γ(2,2) B( 1 2,1 2 ) B(5,1) BCa (-0.28,0.27) (-1.11,1.07) (0.77,1.32) (0.83,1.22) (0.40,0.60) (0.79,0.87) ABC (-0.28,0.27) (-1.11,1.08) (0.77,1.32) (0.83,1.22) (0.40,0.60) (0.79,0.87) Tabla 4.2: Intervalos promedio que proporcionan ambos métodos. empleando el método BCa . Para poder analizar estos hechos, estudiaremos las tablas 4.3, 4.4 y 4.5 que contienen las longitudes medias, formas promedio y coberturas de los dos métodos, respectivamente. N(0,1) N(0,42) Exp(1) Γ(2,2) B( 1 2,1 2 ) B(5,1) BCa 0.544 2.177 0.546 0.390 0.194 0.077 ABC 0.547 2.187 0.549 0.392 0.194 0.077 Tabla 4.3: Longitudes de los intervalos promedio obtenidos por los dos métodos. N(0,1) N(0,42) Exp(1) Γ(2,2) B( 1 2,1 2 ) B(5,1) BCa 1.004 1.004 1.410 1.295 1.002 0.802 ABC 1.004 1.004 1.409 1.295 1 0.800 Tabla 4.4: Forma promedio de los dos tipos de intervalos. Observamos en las tablas 4.3 y 4.4 valores muy similares de simetría y de longitud de los intervalos para los dos métodos. Esta semejanza en estas dos propiedades de los intervalos continúa en la línea de lo visto tabla 4.2, por lo que hasta ahora podemos observar un comportamiento notable de los intervalos ABC . Finalmente nos centraremos en la tabla 4.5 que es la más importante pues reeja cómo de correctos son los intervalos que estamos construyendo. A la vista de la tabla, observamos que la cobertura del método ABC para las seis distribuciones es similar al obtenido empleando el método BCa e incluso en algún caso, como en el de la Beta B(5,1), es notablemente mejor. Por tanto y a la vista de estas tablas podemos concluir que en nuestro estudio de simulación el método ABC se comporta satisfactoriamente al conseguir resultados similares, y en algún caso mejores, que el método BCa , además de reducir considerablemente su coste computacional. 4.4. SOFTWARE DE R ASOCIADO 37 N(0,1) N(0,42) Exp(1) Γ(2,2) B( 1 2,1 2 ) B(5,1) BCa 0.948 0.948 0.943 0.938 0.964 0.932 ABC 0.949 0.949 0.941 0.939 0.962 0.944 Tabla 4.5: Coberturas de los intervalos obtenidos por los dos métodos. 4.4. Software de R asociado Esta sección diere un poco de las tres anteriores, pues el comando boot.ci no permite construir intervalos del tipo ABC . Sin embargo, existe otro comando que nos permite realizar esta función, y ese comando es el abc.ci . Al igual que para el comando boot.ci, tenemos muchos argumentos, de entre los que destacamos: data: La muestra de datos con la que trabajamos. statistic: Función que calcula el estadístico. conf: Nivel de conanza para el que queremos calcular el intervalo. También existen multitud de argumentos secundarios como: index, strata, eps, . . . 38 CAPÍTULO 4. INTERVALOS ABC Anexos 39 46 APÉNDICE A. CÓDIGO DE R DE LAS SIMULACIONES ns=1000 # Número de muestras simuladas interv_teor=matrix(0,nrow=ns,ncol=2) interv_boot=interv_teor interv_perc=interv_teor interv_bca=interv_teor form_teor=c() form_perc=c() form_boot=c() form_bca=c() for (is in 1:ns){ # Bucle de las simulaciones x=rnorm(n) # Datos procedentes de la normal mediat=0 # Media teórica #x=rexp(n) # Datos procedentes de la Exponencial #mediat=1 #x=rbeta(n,1/2,1/2) # Datos procedentes de una beta (1/2,1/2) #mediat=1/2 #x=rgamma(n,2,2) # Datos procedentes de una gamma (2,2) #mediat=1 m=mean(x) dt=sd(x) # Intervalo de la T de Student interv_teor[is,]=c(m-t*dt,m+t*dt) # Intervalo bootstrap estudentizado nb=1000 k=0 tb=c() A.2. CÓDIGO DE R DEL SEGUNDO ESTUDIO DE SIMULACIÓN 47 mb=c() for (ib in 1:nb){ # Bucle de las réplicas bootstrap ind=sample(n,n,replace=TRUE) xb=x[ind] # Muestra bootstrap mb[ib]=mean(xb) if (mb[ib]<m){ k=k+1 } dtb=sd(xb) tb[ib]=(mb[ib]-m)/dtb # Pivote bootstrap estudentizado } z0=qnorm(k/nb) #Cálculo de la estimación del parámetro z0 mj=c() for (j in 1:n){ mj[j]=mean(x[-j]) } mjt=sum(mj)/n U=c() for(l in 1:n){ #Cálculo de la función de influencia jackknife U[l]=mjt-mj[l] } a=sum(U^3)/(6*((sum(U^2))^(3/2))) #Aproximación del parámetro a alfa1=pnorm(z0+(z0+qnorm(alfa/2))/(1-a*(z0+qnorm(alfa/2)))) alfa2=pnorm(z0+(z0+qnorm(1-alfa/2))/(1-a*(z0+qnorm(1-alfa/2)))) qb1=quantile(tb,prob=alfa/2) qb2=quantile(tb,prob=1-alfa/2) 48 APÉNDICE A. CÓDIGO DE R DE LAS SIMULACIONES interv_perc[is,]=c(quantile(mb,alfa/2),quantile(mb,1-alfa/2)) interv_boot[is,]=c(m-qb2*dt,m-qb1*dt) interv_bca[is,]=c(quantile(mb,alfa1),quantile(mb,alfa2)) ########## FORMA DE LOS INTERVALOS ########## form_teor[is]=(interv_teor[is,2]-m)/(m-interv_teor[is,1]) form_perc[is]=(interv_perc[is,2]-m)/(m-interv_perc[is,1]) form_boot[is]=(interv_boot[is,2]-m)/(m-interv_boot[is,1]) form_bca[is]=(interv_bca[is,2]-m)/(m-interv_bca[is,1]) ############################################# if (is==100*round(is/100)){ cat( ' is= ' ,is,interv_teor[is,],interv_perc[is,], interv_boot[is,],interv_bca[is,], ' \n ' ) } } ## is= 100 -0.139 0.470 -0.132 0.471 -0.136 0.485 -0.133 0.469 ## is= 200 -0.370 0.296 -0.351 0.302 -0.372 0.297 -0.343 0.320 ## is= 300 -0.515 0.065 -0.505 0.056 -0.507 0.075 -0.505 0.057 ## is= 400 -0.263 0.268 -0.267 0.254 -0.270 0.277 -0.269 0.254 ## is= 500 -0.238 0.289 -0.224 0.278 -0.240 0.275 -0.217 0.284 ## is= 600 -0.050 0.608 -0.049 0.584 -0.064 0.602 -0.064 0.579 ## is= 700 -0.439 0.074 -0.446 0.065 -0.412 0.138 -0.437 0.075 ## is= 800 -0.546 0.129 -0.514 0.119 -0.541 0.112 -0.505 0.140 ## is= 900 -0.316 0.304 -0.301 0.302 -0.324 0.304 -0.292 0.307 ## is= 1000 -0.402 0.258 -0.383 0.247 -0.381 0.260 -0.399 0.226 ########## COBERTURA DE LOS INTERVALOS ########## A.2. CÓDIGO DE R DEL SEGUNDO ESTUDIO DE SIMULACIÓN 49 cobertura_teor=sum((interv_teor[,1]<mediat)&(interv_teor[,2]>mediat))/ns cobertura_teor ## [1] 0.956 cobertura_perc=sum((interv_perc[,1]<mediat)&(interv_perc[,2]>mediat))/ns cobertura_perc ## [1] 0.943 cobertura_boot=sum((interv_boot[,1]<mediat)&(interv_boot[,2]>mediat))/ns cobertura_boot ## [1] 0.957 cobertura_bca=sum((interv_bca[,1]<mediat)&(interv_bca[,2]>mediat))/ns cobertura_bca ## [1] 0.948 ########## PROMEDIO DE LOS INTERVALOS ########## prom_teor=c(sum(interv_teor[,1])/ns,sum(interv_teor[,2])/ns) prom_teor ## [1] -0.2873495 0.2783529 prom_perc=c(sum(interv_perc[,1])/ns,sum(interv_perc[,2])/ns) prom_perc ## [1] -0.2763331 0.2671190 prom_boot=c(sum(interv_boot[,1])/ns,sum(interv_boot[,2])/ns) prom_boot ## [1] -0.2865110 0.2783045 prom_bca=c(sum(interv_bca[,1])/ns,sum(interv_bca[,2])/ns) prom_bca ## [1] -0.2767618 0.2673660 50 APÉNDICE A. CÓDIGO DE R DE LAS SIMULACIONES ########## FORMA DE LOS INTERVALOS ########## form_teor_prom=sum(form_teor)/ns form_teor_prom ## [1] 1 form_perc_prom=sum(form_perc)/ns form_perc_prom ## [1] 1.001297 form_boot_prom=sum(form_boot)/ns form_boot_prom ## [1] 1.007747 form_bca_prom=sum(form_bca)/ns form_bca_prom ## [1] 1.004354 ########## LONGITUD DE LOS INTERVALOS ########## long_teor=prom_teor[2]-prom_teor[1] long_teor ## [1] 0.5657024 long_perc=prom_perc[2]-prom_perc[1] long_perc ## [1] 0.5434521 long_boot=prom_boot[2]-prom_boot[1] long_boot ## [1] 0.5648155 long_bca=prom_bca[2]-prom_bca[1] long_bca ## [1] 0.5441279 A.3. CÓDIGO DE R DEL TERCER ESTUDIO DE SIMULACIÓN 51 A.3. Código de R del tercer estudio de simulación En esta sección trataremos el código R del estudio de simulación de la sección 4.3. Por simplicidad trataremos solo el caso en el que la muestra provenga de una distribución normal N(0,1). Los demás casos serán análogos. set.seed(1707) #Creamos una semilla para aleatorizar los datos library(boot) nivel=0.95 #El nivel de confianza que queremos para nuestro intervalo alfa=1-nivel #El nivel de significación del mismo n=50 #Tamaño de la muestra t=qt(1-alfa/2,df=n-1)/sqrt(n) #Pivote teórico ns=1000 # Número de muestras simuladas interv_bca=matrix(0,nrow=ns,ncol=2) interv_abc=interv_bca form_bca=c() form_abc=c() for (is in 1:ns){ # Bucle de las simulaciones x=rnorm(n) # Datos procedentes de la normal mediat=0 # Media teórica #x=rnorm(n,mean=0,sd=4) # Datos procedentes de la normal N(0,4^2) #mediat=0 # Media teórica #x=rexp(n) # Datos procedentes de la Exponencial #mediat=1 #x=rgamma(n,2,2) # Datos procedentes de una gamma (2,2) #mediat=1 52 APÉNDICE A. CÓDIGO DE R DE LAS SIMULACIONES #x=rbeta(n,1/2,1/2) # Datos procedentes de una beta (1/2,1/2) #mediat=1/2 #x=rbeta(n,5,1) # Datos procedentes de una beta (5,1) #mediat=5/6 m=mean(x) dt=sd(x) i=rep(1/n,n) media<-function(x,i){ sum(x*i) } interv_abc[is,]=abc.ci(x,media,conf=0.95,eps=0.001/n)[c(2,3)] nb=1000 k=0 tb=c() mb=c() for (ib in 1:nb){ # Bucle de las réplicas bootstrap ind=sample(n,n,replace=TRUE) xb=x[ind] # Muestra bootstrap mb[ib]=mean(xb) if (mb[ib]<m){ k=k+1 } } z0=qnorm(k/nb) #Cálculo de la estimación del parámetro z0 mj=c() for (j in 1:n){ mj[j]=mean(x[-j]) } A.3. CÓDIGO DE R DEL TERCER ESTUDIO DE SIMULACIÓN 53 mjt=sum(mj)/n U=c() for(l in 1:n){ #Cálculo de la función de influencia jackknife U[l]=mjt-mj[l] } a=sum(U^3)/(6*((sum(U^2))^(3/2))) #Aproximación del parámetro a alfa1=pnorm(z0+(z0+qnorm(alfa/2))/(1-a*(z0+qnorm(alfa/2)))) alfa2=pnorm(z0+(z0+qnorm(1-alfa/2))/(1-a*(z0+qnorm(1-alfa/2)))) interv_bca[is,]=c(quantile(mb,alfa1),quantile(mb,alfa2)) ########## FORMA DE LOS INTERVALOS ########## form_bca[is]=(interv_bca[is,2]-m)/(m-interv_bca[is,1]) form_abc[is]=(interv_abc[is,2]-m)/(m-interv_abc[is,1]) ############################################# if (is==100*round(is/100)){ cat( ' is= ' ,is,interv_bca[is,],interv_abc[is,], ' \n ' ) } } ## is= 100 -0.1326414 0.4685368 -0.1175815 0.4710775 ## is= 200 -0.3431239 0.3204149 -0.3537932 0.2887397 ## is= 300 -0.5047584 0.05687106 -0.4984203 0.06173165 ## is= 400 -0.2688639 0.2542557 -0.2566898 0.2561227 ## is= 500 -0.2168557 0.2843857 -0.2332563 0.2754397 ## is= 600 -0.06392886 0.5792595 -0.05800721 0.5783905 ## is= 700 -0.4369793 0.07525489 -0.4097825 0.08829007 ## is= 800 -0.5053228 0.1402596 -0.5312045 0.1205397 54 APÉNDICE A. CÓDIGO DE R DE LAS SIMULACIONES ## is= 900 -0.2916933 0.3070372 -0.2993186 0.2997511 ## is= 1000 -0.3986938 0.2256725 -0.3829054 0.2549409 ########## COBERTURA DE LOS INTERVALOS ########## cobertura_bca=sum((interv_bca[,1]<mediat)&(interv_bca[,2]>mediat))/ns cobertura_bca ## [1] 0.948 cobertura_abc=sum((interv_abc[,1]<mediat)&(interv_abc[,2]>mediat))/ns cobertura_abc ## [1] 0.949 ########## PROMEDIO DE LOS INTERVALOS ########## prom_bca=c(sum(interv_bca[,1])/ns,sum(interv_bca[,2])/ns) prom_bca ## [1] -0.2767618 0.2673660 prom_abc=c(sum(interv_abc[,1])/ns,sum(interv_abc[,2])/ns) prom_abc ## [1] -0.2776082 0.2691058 ########## FORMA DE LOS INTERVALOS ########## form_bca_prom=sum(form_bca)/ns form_bca_prom ## [1] 1.004354 form_abc_prom=sum(form_abc)/ns form_abc_prom ## [1] 1.004251 A.3. CÓDIGO DE R DEL TERCER ESTUDIO DE SIMULACIÓN 55 ########## LONGITUD DE LOS INTERVALOS ########## long_bca=prom_bca[2]-prom_bca[1] long_bca ## [1] 0.5441279 long_abc=prom_abc[2]-prom_abc[1] long_abc ## [1] 0.546714