Contrastes de especificación para modelos de distribución

Traballo Fin de Grao Contrastes de especicación para modelos de distribución Carlos Rodríguez Ameal 2018/2019 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Contrastes de especicación para modelos de distribución Carlos Rodríguez Ameal Julio de 2019 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Trabajo propuesto Área de Coñecemento: Estadística e Investigación Operativa Título: Contrastes de especicación para modelos de distribución Breve descrición do contido Se trata de revisar los contrastes de especicación más notables para modelos de distribución de una variable aleatoria discreta o continua. 1. Contrastes de especicación basados en la estimación de una distribución multinomial 2. Contrastes de especicación basados en la estimación de la función de distribución 3. Contrastes de especicación basados en la estimación de la función de densidad 4. Ilustración en bases de datos simulados iii Índice general Resumen viii Introducción xi 1. Contrastes basados en la estimación de una distribución multinomial 1 1.1. Introducción: el contraste χ2 dePearson.................... 1 1.2. Marcoteórico................................... 2 1.2.1. La distribución multinomial . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.2. Test de razón de verosimilitudes . . . . . . . . . . . . . . . . . . . . 5 1.2.3. Test de razón e verosimilitudes aplicado a la multinomial . . . . . . 6 1.3. Contrastes de especicación basados en la multinomial . . . . . . . . . . . . 9 1.3.1. Hipótesissimple ............................. 9 1.3.2. Hipótesis compuesta . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.4. La familia de estadísticos de divergencia . . . . . . . . . . . . . . . . . . . . 12 1.5. Comentarios sobre el método para distribuciones continuas . . . . . . . . . . 14 2. Contrastes basados en la función de distribución empírica 15 2.1. La función de distribución empírica . . . . . . . . . . . . . . . . . . . . . . . 15 2.2. Elprocesoempírico................................ 17 2.3. El estadístico de Kolmogorov-Smirnov . . . . . . . . . . . . . . . . . . . . . 19 2.3.1. El estadístico KS para hipótesis compuestas . . . . . . . . . . . . . . 21 2.3.2. El estadístico KS aplicado al caso discreto . . . . . . . . . . . . . . . 22 2.4. Estadísticos de Anderson-Darling . . . . . . . . . . . . . . . . . . . . . . . . 24 2.4.1. Comentariosnales............................ 25 3. Contrastes basados en la estimación de la función de densidad 27 3.1. Estimación de la función de densidad . . . . . . . . . . . . . . . . . . . . . . 28 3.1.1. Elhistograma............................... 28 v vi ÍNDICE GENERAL 3.1.2. Estimación de densidad kernel . . . . . . . . . . . . . . . . . . . . . 29 3.1.3. Análisis del estimador . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.2. Contraste de hipótesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.2.1. Estadísticos de contraste . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.2.2. Comportamiento asintótico . . . . . . . . . . . . . . . . . . . . . . . 35 3.2.3. Hipótesis compuesta . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.3. Comentariosnales................................ 37 4. Ilustración en bases de datos simulados 39 4.1. Hipótesissimple.................................. 40 4.1.1. Contraste usando los estadísticos de divergencia . . . . . . . . . . . . 40 4.1.2. Contraste con los estadísticos basados en el proceso empírico . . . . 43 4.1.3. Contrastes con los estadísticos basados el la estimación de la función dedensidad................................ 46 4.2. Hipótesiscompuesta ............................... 49 Conclusión 53 Código de R 57 .1. Hipótesissimple.................................. 57 .1.1. Primercapítulo.............................. 57 .1.2. Segundocapítulo............................. 59 .1.3. Tercercapítulo .............................. 63 .2. Hipótesiscompuesta ............................... 65 .2.1. Primercapítulo.............................. 65 .2.2. Capítulo2................................. 68 Bibliografía 75 xiv INTRODUCCIÓN Capítulo 1 Contrastes basados en la estimación de una distribución multinomial El más antiguo y quizás mejor conocido de entre todos los contrastes de especicación existentes es el de la χ2 de Pearson, introducido en su famoso trabajo de 1900. Su relevancia histórica es indiscutible, ya que durante varias décadas de principios del siglo XX fue el único test de bondad de ajuste disponible, empleándose universalmente no solo para el contraste de multinomiales, sino que fue más tarde modicado para distribuciones continuas. Estas eran transformadas en multinomiales dividiendo su conjunto de llegada en distintos intervalos y hallando las probabilidades asociadas a cada uno. Es evidente que esta forma de proceder supone una pérdida importante de información, y, en general, es más aconsejable utilizar los otros tipos de contrastes que estudiaremos en los próximos capítulos, mejor adaptados a distribuciones continuas. Sin embargo, la importancia de la distribución multinomial, junto con la relevancia histórica y conceptual del método y el hecho de que varios de los tests más modernos aún se basen en él, justican la inclusión para el estudio de esta familia de contrastes en el presente trabajo. 1.1. Introducción: el contraste χ2 de Pearson Supongamos que tenemos una muestra que sospechamos proviene de una variable aleatoria discreta cuya distribución conocemos. Nos referiremos a cada posible valor de la muestra a través de un índice i , y por Oi denotamos el número de veces que hemos obtenido ese resultado. Para estudiar si se cumple nuestra hipótesis, lo natural sería comparar de alguna forma dichas observaciones con los resultados más probables que obtendríamos en el caso de que nuestra suposición sea cierta. A estas cantidades las denotamos Ei . Esta idea sencilla es la que tiene Pearson cuando propone su famoso estadístico, que 1 2 CAPÍTULO 1. CONTRASTES BASADOS EN LA MULTINOMIAL tiene la siguiente expresión: X2=X i (Oi−Ei)2 Ei (1.1) La intuición de Pearson es bastante clara. El estadístico acumula las diferencias al cuadrado entre las observaciones y los resultados esperados, pero multiplicadas por un factor de ponderación: el inverso de lo esperado. El motivo de la elección de este factor de ponderación es bastante fácil de comprender: no es lo mismo que la diferencia sea de 3 cuando esperamos 4 observaciones que cuando esperamos 40 . Si bien la fórmula del estadístico de Pearson parece bastante directa, lo novedoso de su trabajo fue que en su célebre artículo de 1900 llegó a describir correctamente el comportamiento límite que tomaba: una χ2 . Aunque el descubrimiento en sí de esta distribución no se le concede, su trabajo supuso un punto crucial en la estadística matemática, que en su época se encontraba estancada en la cuestión de hasta qué punto las distribuciones estándares, concretamente la normal, servían para modelizar universalmente los diferentes procesos aleatorios. Alejándose de esa postura, llega al test de la χ2 , dando lugar al primer contraste de especicación propiamente dicho. Pearson no llegó a hablar de distribución multinomial, pero esta está implícita en su comprensión de los datos esperables frente a los observados. Nosotros, desde un enfoque más moderno, utilizaremos las propiedades de la multinomial para demostrar la convergencia asintótica del estadístico de Pearson a una χ2 . Además, daremos un enfoque alternativo basándonos en el test de razón de verosimilitudes, popularizado décadas más tarde por Fisher. Justamente, este último debe su interés por la bioestadística (la genética mendeliana y la biología evolutiva fueron dos campos que motivaron numerosos avances en las técnicas estadísticas de principios del siglo XX) a la lectura de las obras de Pearson, y ambos protagonizaron una sonada polémica al discrepar acerca del comportamiento del X2 cuando se usa para contrastar hipótesis compuestas. Estas cuestiones teóricas las introduciremos a continuación. 1.2. Marco teórico 1.2.1. La distribución multinomial Supongamos que tenemos una variable aleatoria X que puede tomar k resultados diferentes, de forma que cada suceso i tiene una probabilidad de ocurrir pi . Cada uno de sus resultados obtenidos viene modelizado por una variable categórica, mientras que el vector que mide las frecuencias observadas de los resultados obtenidos en n intentos se denomina 1.2. MARCO TEÓRICO 3 multinomial. Así, una variable categórica supone una generalización de una Bernouilli, y una multinomial de una binomial, cuando pasamos de 2 a k posibles resultados. Al igual que podemos denir la distribución binomial como una suma de variables Bernouilli independientes, la multinomial puede expresarse como una suma de variables categóricas i.i.d. 1 Más formalmente: Denición 1.1. Denimos la variable categórica asociada al suceso Xr de una variable aleatoria X , r∈ {1,2, ..., n} , como el vector: ξr= (I{Xr=1}, I{Xr=2}, ..., I{Xr=k}), (1.2) donde I{Xr=i} es la función indicatriz asociada al resultado i de la variable aleatoria Xr . Por tanto, el vector toma el valor 1 en la posición i -ésima y 0 en todas las demás con probabilidad pi=P(Xr=i) . La suma de n variables categóricas resulta en una multinomial, que expresamos como sigue: Denición 1.2. Denimos la multinomial de tamaño k y con parámetros n y p = (p1, p2, ..., pn) tal que Pn i=1 pi= 1 , como el vector aleatorio N = (N1, N2, ..., Nk) con función de probabilidad: P{N1=n1, N2=n2, ..., Nk=nk}=n! n1!n2!...nk!pn1pn2...pnk (1.3) Una forma alternativa de denir ambas variables consiste en suprimir la última coordenada de los vectores (que pasarán así a ser k−1 dimensionales), de forma que si el resultado aleatorio r toma un valor en la k -ésima categoría, la variable categórica correspondiente tenga todas las coordenadas iguales a 0 . En cuanto a la variable multinomial, si queremos ver cuántos resultados han ocurrido en la categoría k tras n intentos, no tenemos más que hallar la resta n−Pk−1 i=1 ni . Ambas representaciones son equivalentes y las utilizaremos indistintamente según nos convenga. Ahora enunciaremos las propiedades de la distribución categórica cuya demostración es inmediata: Lema 1.3. Sea ξr= (I{Xr=1}, I{Xr=2}, ..., I{Xr=k}) una variable categórica con parámetro p = (p1, p2, ..., pn) . Se verica: 1. E(ξr) = p 2. Su matriz de covarianzas Σ=(σ)ij cumple: 1 Independientes e idénticamente distribuidas 4 CAPÍTULO 1. CONTRASTES BASADOS EN LA MULTINOMIAL σij =E[I{Xr=i}·I{Xr=j}]−pipj=(−pipjsi i 6=j pi(1 −pi)si i =j A partir de la expresión de la multinomial como suma de variables categóricas i.i.d. N =Pn i=1 ξi y de las propiedades de la esperanza y la covarianza podemos probar el siguiente lema: Lema 1.4. Sea N = (N1, N2, ..., Nk) una multinomial de parámetros n y p = (p1, p2, ..., pn) . Se verica: 1. Ni∈Binomial(n;pi) para i= 1,2, ..., k . 2. E(Ni) = npi para i= 1,2, ..., k 3. V ar(Ni) = npi(1 −pi) para i= 1,2, ..., k 4. Cov(Ni, Nj) = −npipj Para terminar, presentaremos un resultado sobre la convergencia de la multinomial que nos será útil más adelante para estudiar el comportamiento asintótico nulo del estadístico de Pearson. Dado que tenemos la restricción: Pk i=1 Ni=n , la matriz de covarianzas del vector tiene rango n−1 . Para solventar este contratiempo, emplearemos la denición alternativa de la multinomial bajo la forma de un vector de dimensión k−1 . Teorema 1.5. Sea ξr= (I{Xr=1}, I{Xr=2}, ..., I{Xr=k−1}) la variable categórica asociada al resultado Xr y sean p y Σ su esperanza y matriz de covarianzas respectivamente. Si N= (N1, N2, ..., Nk−1) = n X r=1 ξr es el vector multinomial que acumula n resultados, se cumple cuando n→ ∞ que: N−np √nΣ−1/2d −→ Nk−1(0, I) . 2 , Demostración. Denamos ηr= (ξr− p )Σ−1/2 , que tiene media igual a 0 y matriz de covarianzas Ik−1 . Entonces, por el Teorema Central del Límite: N−np √nΣ−1/2=1 √n n X r=1 ηrd −→ Nk−1(0, I) 2 Por Xn d →X denotamos la convergencia de la sucesión de variabes Xn a la variable X en distribución. 1.2. MARCO TEÓRICO 5 1.2.2. Test de razón de verosimilitudes Al principio de este capítulo hemos introducido el estadístico de Pearson, nacido de su intuición matemática aplicada al problema de cómo medir la discrepancia entre una muestra y una distribución esperada. Pearson fue capaz de llegar a la distribución límite bajo la nula de este estadístico particular usando sus conocimientos de probabilidad, y de crear uno de los primeros test de contraste: el de la χ2 . Sin embargo, a pesar de la enorme importancia que su trabajo tuvo en su época y de la inuencia que sigue teniendo hoy en día, si el mismo problema se plantease a un estadístico actual, probablemente la primera forma de abordar el contraste que se le ocurriese sería utilizando uno de los métodos más universales: el test de razón de verosimilitudes. Lo interesante es que cuando este método fue creado y estudiado, sirvió para conrmar el resultado al que Pearson había llegado tres décadas antes. El test se aplica siempre que tengamos una muestra que sigue una distribución Fθ con fución de densidad fθ , especicada excepto por un parámetro θ que varía en un espacio paramétrico Θ⊂Rk , y queremos contrastar H0:θ∈Θ0 frente a H1:θ∈Θ\Θ0 donde Θ0⊂Θ . El método se basa en el concepto de verosimilitud de la muestra: fθ(x1, x2, ..., xn) , que funciona como una medida de lo bien que explica θ los resultados obtenidos. Así, sup θ∈Θ0 fθ(x1, ..., xn) supone un índice de la mejor explicación de la muestra bajo H0 , y el cual vamos a querer comparar con sup θ∈Θ fθ(x1, ..., xn) , que se corresponder con la mejor explicación posible que existe entre todos los valores posibles del parámetro. Los estimadores máximo verosímil serán los parámetros asociados a estos valores, respectivamente: ˆ θ0 y ˆ θ tales que fˆ θ0(x1, ..., xn) = sup θ∈Θ0 fθ(x1, ..., xn) y fˆ θ(x1, ..., xn) = sup θ∈Θ fθ(x1, ..., xn) . Para comparar ambas cantidades calculamos el cociente entre ellas, que es lo que denominamos rezón de verosimilitudes: Λ(x1, ..., xn) = sup θ∈Θ0 fθ(x1, ..., xn) sup θ∈Θ fθ(x1, ..., xn)=fˆ θ0(x1, ..., xn) fˆ θ(x1, ..., xn) (1.4) Un test de razón de verosimilitudes rechazará la hipótesis nula cuando Λ(x1, ..., xn)< c para un c que jaremos en función del nivel de signicación que busquemos. Para ello podemos intentar hallar la distribución exacta de Λ , que en ocasiones puede ser expresado como un estadístico más simple. Sin embargo, esto es muy difícil en general, por lo que es necesario contar con algún resultado que nos dé información sobre su comportamiento asintótico. Más precisamente, usaremos el fenómeno de Wilks, que caracteriza la convergencia del logritmo de este cociente. 6 CAPÍTULO 1. CONTRASTES BASADOS EN LA MULTINOMIAL Teorema 1.6. Supongamos que la hipótesis nula dependa de q parámetros, es decir: Θ0= {θ∈Θ : θi=gi(w1, ..., wq)para i = 1, ..., k;con (w1, ..., wq)∈Ω} siendo Ω un abierto de Rq y gi funciones con derivadas parciales de orden 1 continuas. Bajo ciertas hipótesis de regularidad (se puede consultar el manual Principios de Inferencia Estadística de Ricardo Vélez para más información), y cuando n→ ∞ , si la hipótesis nula es cierta se cumple que : −2logΛ(x1, ..., xn)d →χ2 k−q Es decir, conforme la muestra se hace más grande, el estadístico −2logΛ(x1, ..., xn) converge a una χ2 cuyos grados de libertad se corresponden a la diferencia de las dimensiones entre Θ y Θ0 . No hemos explicitado las hipótesis de regularidad del manual que estamos siguiendo (el de Ricardo Vélez) puesto que el propio autor arma que estas son sucientes pero no necesarias, y que la validez del resultado se suele aceptar sin reparos en los contextos usuales. Gracias a este resultado, ya podemos denir nuestro test, que rechazará la hipótesis nula con un nivel de signicación α siempre que: −2logΛ> χ2 k−1;α(⇔Λ< e −χ2 k−1;α 2) 1.2.3. Test de razón e verosimilitudes aplicado a la multinomial Ahora nos interesa aplicar la teoría de la sección anterior al caso particular de una multinomial. Supongamos que nuestra muestra x1, ..., xn proviene de una variable discreta que toma k valores diferentes y tiene a N como vector multinomial con parámetro p , y que queremos contrastar H0: p = p 0 frente a H1: p 6= p 0 . Usando la notación que hemos introducido antes, tenemos que Θ y Θ0 tienen dimensión k−1 ( p ∈Rk con la restricción Pk i=1 pi= 1 ) y 0 respectivamente. En primer lugar, para calcular la razón de verosimilitudes, veamos que el estimador máximo verosímil del parámetro pi es ˆpi=Ni/n . De aquí en adelante denotaremos de la misma manera a cada una de las componentes del vector aleatorio Ni , como al número de observaciones de la muestra iguales a i (es decir, no distinguiremos entre Ni y ni ). Entonces, el problema de hallar el estimador máximo verosímil se corresponde con hallar los parámetros p de N que hacen más probable los Ni , es decir, que maximizan la función de probabilidad: ˆ p =ArgMax p ∈Rk,Pk i=1 pi=1,pi≥0 n! N1!N2!...Nk!pN1 1pN2 2...(1 −p1−p2−... −pk−1)Nk Aplicando logaritmo a la función de probabilidad (el logaritmo es creciente e inyectivo) e igualando las derivadas parciales respecto de los pi a cero para buscar los argumentos máximos, obtenemos las ecuaciones de verosimilitud: 1.2. MARCO TEÓRICO 7 (Nj pj−Nk 1−p1−p2−...−pk−1= 0 j= 1,2, ..., k −1 por lo que la solución verica: nk 1−ˆp1−ˆp2−...−ˆpk−1=n1 ˆp1=n2 ˆp2=... =nk−1 ˆpk−1 De un sencillo cálculo obtenemos la relación ˆpi=Ni/n . La matriz de derivadas segundas es denida positiva, con lo que es un máximo relativo, y como en la frontera alguna de las probabilidades se anula, así lo hace la función de densidad y concluimos que el máximo es global. Ahora podemos obtener la razón de verosimilitudes: Λ(N1, N2, ..., Nk) = (p0 1)N1(p0 2)N2...(p0 k)Nk ˆpN1 1ˆpN2 2...ˆpNk k = k Y i=1 p0 i ˆpiNi (1.5) Podemos denir entonces un región crítica: {Λ< c}={−2 log Λ > h} , que es la región de valores para los cuales rechazamos la hipótesis nula. Para poder elegir valor k en función del nivel de signicación que queramos nos basamos en el estadístico: G2=−2 log Λ = 2 k X i=1 Ni(log ˆpi−log p0 i), (1.6) que tiene una distribución asintótica χ2 k−1 suponiendo la hipótesis nula H0:p=p0 como consecuencia del teorema 1.6. Es decir, este estadístico se comporta de misma forma que el ideado por Pearson en el contexto de los contrastes de especicación. Recordemos que este último se denía como la suma de las diferencias al cuadrado entre los valores observados y los esperados divididas por los esperados. La formalización de la multinomial nos aporta una nueva notación para este estadístico más acorde con el tono del trabajo: X2= k X i=1 (Ni−np0 i)2 np0 i (1.7) Es sorprendente que, a pesar de nacer en dos contextos totalmente diferentes, ambos estadísticos son muy parecidos en la práctica. La clave está en la intuición de Pearson al elegir como factores de ponderación los 1/np0 i , pues admitiendo que (p0 i−ˆpi)/ˆpi sean pequeños: logp0 i ˆpi=log 1 + p0 i−ˆpi ˆpi≃p0 i−ˆpi ˆpi−1 2p0 i−ˆpi ˆpi2 8 CAPÍTULO 1. CONTRASTES BASADOS EN LA MULTINOMIAL y por tanto: logΛ≃ k X i=1 Ni p0 i−ˆpi ˆpi−1 2 k X i=1 (p0 i−ˆpi)2 ˆp2 i =n k X i=1 (P0 i−ˆpi)−1 2 k X i=1 n2 Ni (ˆpi−p0 i)2= −1 2 k X i=1 (Ni−np0 i)2 Ni de forma que: −2logΛ≃ k X i=1 (Ni−np0 i)2 Ni≃ k X i=1 (Ni−np0 i)2 np0 i (1.8) puesto que Ni≃np0 i para i= 1, ..., k . Este similitud provoca, como ya hemos dicho, un mismo comportamiento asintótico para ambos estadísticos. Para G2 este resultado se prueba echando mano de la teoría del test de razón de verosimilitudes. Pero es fácil de ver para el estadístico de Pearson con unos pocos cálculos. Ello da lugar al siguiente teorema: Teorema 1.7. Cuando n→ ∞ y si la hipótesis nula es cierta: X2= k X i=1 (Ni−np0 i)2 np0 i d →χ2 k−1 Demostración. Por el teorema 1.5 y simplemente usando la denición de la χ2 tenemos: N−np √nΣ−1N−np √nd −→ χ2 k−1 Pero calculando la matriz Σ−1 y desarrollando esta producto (consúltese el manual de Ricardo Vélez para más detalles) se llega precisamente a la igualdad: N−np √nΣ−1N−np √n=X2 Contamos entonces con dos test que rechazarán la hipótesis nula con un nivel de signicación α cuando su estadístico correspondiente tome un valor mayor que χ2 1−k;α . La precisión de estos test depende en buena medida de cómo de buena sea la aproximación de la χ2 para las distribuciones exactas de los estadísticos. El criterio más extendido que se suele pedir es que los valores esperados cumplan np0 i>5 , que se suele combinar con restricciones sobre el número de obervaciones (normalmente se pedirá n > 30 ) y el de intervalos (cuando no imposibilite la primera pauta se cogerá k > 5 ). 1.3. CONTRASTES DE ESPECIFICACIÓN BASADOS EN LA MULTINOMIAL 9 1.3. Contrastes de especicación basados en la multinomial 1.3.1. Hipótesis simple Supongamos que tenemos n datos de una muestra aleatoria simple que sigue una distribución desconocida F , y que deseamos contrastar la hipótesis nula H0:F=F0 , donde F0 es una función completamente especicada (no depende de ningún parámetro de valor desconocido). De esta forma la hipótesis nula es simple, y su alternativa H1:F6=F0 está compuesta por todas las distribuciones distintas a F0 . Recordamos que este es el caso más sencillo, en el que en vez de contrastar si la muestra ha sido obtenida de una familia de distribuciones se hace para una F0 de la cuál conocemos su función de distribución, y por tanto la probabilidad p0 i asociada a cualquier intervalo Ai en el que tome valores. Dividiendo entonces el recorrido en A1 , A2 ,..., Ak intervalos disjuntos, inmediatamente podemos considerar el vector aleatorio N = (N1 , N2 ,..., Nk ), con Ni igual al número de observaciones de la muestra en cada subconjunto Ai , y que sigue una distribución multinomial con cada parámetro pi igual a la probabilidad de que F dé un valor en Ai . La idea consiste en sustituir la hipótesis inicial por H0:p=p0 frente a H1:p6=p0 , siendo p0= (p0 1, p0 2, ..., p0 k) el vector con las probabilidades de que F0 caiga en cada intervalo, y dando lugar al contraste paramétrico ya visto, y que podemos realizar usando los estadísticos X2 o G2 . Más tarde presentaremos una familia general de estadísticos de contraste y haremos ciertos comentarios sobre sus propiedades. Está claro que esta forma de proceder tiene una desventaja fundamental. Cuando la distribución hipotética que queremos contrastar es discreta, los intervalos pasan a ser puntos y el contraste de especicación se corresponde desde un primer momento con el de la multinomial que hemos estudiado en la sección anterior. Sin embargo, cuando nuestra F0 es continua, el método se basa en una discretización de su recorrido, es decir, la transformamos en una multinomial con toda la pérdida de información que eso conlleva. Podemos dar innitos ejemplos de funciones con distribuciones notadamente diferentes que dan lugar a multinomiales idénticas si los intervalos se eligen de la manera adecuada (pero lo cual no es nada adecuado para nuestro propósito). Así, aceptar la hipótesis nula del contraste multinomial no garantiza en la realidad que podamos aceptar la hipótesis nula del contraste original. La solución a este defecto pasa por elegir el número de intervalos y la posición de sus fronteras de la forma adecuada. Al nal del capítulo comentaremos con detenimiento el proceso del contraste aplicado a distribuciones continuas y intentaremos dar ciertas pautas a seguir. 16 CAPÍTULO 2. CONTRASTES BASADOS EN LA FUNCIÓN DE DISTRIBUCIÓN ˆ Fn(x) = 1 n#{Xi∈Sn:Xi≤x}=1 n n X i=1 I(Xi≤x) 1 Es inmediato ver que ˆ Fn es una función no decreciente y escalonada, con cada escalón extendiéndose entre xi−1 y xi y lim x→−∞F(x) = 0 , lim x→∞F(x) = 1 . También cumple que es continua por la derecha y tiene límite por la izquierda de todo punto x , con lo que en particular ˆ Fn es una función de distribución. Ejemplo de función de distribución empírica Figura 2.1: Los segmentos negros se corresponden con los escalones de la función de distribución empírica obtenida a partir de 20 datos simulados de una normal estándar. Su función de distribución real se representa con una línea roja punteada Ahora nos centraremos en el comportamiento puntual de ˆ Fn . Dado que nˆ Fn(x) es la función que cuenta el número de observaciones de la muestra Sn menores o iguales que x , y que cada dato tiene una probabilidad F(x) de ser menor o igual que x , es obvio que nˆ Fn(x) sigue una Binomial de parámetros n y F(x) (de hecho al denirla la hemos expresado como la suma de variables Bernouilli). De aquí deducimos las siguientes propiedades: Lema 2.2. 1. E(ˆ Fn(x)) = F(x)∀x, ∀n 2. Por la ley fuerte de los grandes números, cuando n→ ∞ ˆ Fn(x)a.s. −→ F(x)∀x 2 3. Por el Teorema Central de Límite, cuando n→ ∞ √n(ˆ Fn(x)−F(x)) d −→ N(0, F(x)(1 −F(x)) ∀x 1 I(Xi≤x) es el indicador del evento (Xi≤x) , que jada x se comporta como una Bernouilli con p=F(x) 2 a.s. signica almost sure en inglés y se usa para denotar que la convergencia es casi segura. 2.2. EL PROCESO EMPÍRICO 17 Como corolario inmediato de la propiedad (2) , podemos armar que la FDE es un estimador consistente de F(x) en cada punto x . Estas propiedades nos dan información sobre el comportamiento asintótico de ˆ Fn(x) , pero hacen referencia a convergencias puntuales. Sin embargo, la propiedad (2) puede ser extendida para la convergencia uniforme de ˆ Fn(x) a F(x) . Esto es lo que arma el Teorema de Givenko-Cantelli: Teorema 2.3. Givenko-Cantelli. Con la notación anterior, y cuando n→ ∞ , sup x|(ˆ Fn(x)−F(x))|a.s. −→ 0 Estos resultados aseguran que ˆ Fn se acerca cada vez más a F cuando la muestra es grande. Por ello, si queremos realizar el contraste H0:F(x) = F0(x) , tiene sentido medir la diferencia entre ˆ Fn(x) y F0(x) . Esta es justamente la idea subyacente de los test que veremos en esta sección, cuyos estadísticos serán de la forma: Tn=c(n)d(ˆ Fn, F0), (2.1) siendo c(n) un factor de escala y d(., .) una distancia o función de divergencia 3 entre funciones. Nótese que ˆ Fn es una función de x al igual que F0 , por lo que tiene sentido medir esta diferencia. Por denición, se cumple que d(F, F0)=0 4 ⇔H0 es cierta. Además, notemos que tal y como ˆ Fn es un estimador de de F , d(ˆ Fn, F0) es un estimador de d(F, F0) , con lo que siempre trabajaremos con medidas cuyos estimadores sean consistentes. Esto implica que cuando n se hace grande, d(ˆ Fn, F0) bajo la nula tiende a 0 en probabilidad. Una vez que elegimos d podemos estudiar las propiedades de Tn , para lo que suele ser crucial el comportamiento asintótico de √n(ˆ Fn(x)−F(x)) . Sin embargo, los resultados de convergencia puntual no suelen ser sucientes para la mayor parte de los Tn . 2.2. El proceso empírico A pesar de que la FDE constituye la base sobre la cual se articula la familia de test que estudiamos en este capítulo, suele ser más conveniente trabajar directamente con el 3 Una función de divergencia es más débil que una distancia en el sentido de que no es necesariamente simétrica ni cumple la desigualdad triangular. 4 Esta igualdad puede ser entendida como la igualdad estricta de funciones (que podría ser contrastada con el estadístico de Kolmogorov-Smirnov), o, más en general, como la igualdad de las clases de F y F0 en el espacio de funciones Ln (a la que se limita los estadísticos de Anderson-Darling). Como estamos en un contexto probabilístico nos llega con que se cumpla la igualdad para casi todo punto, con lo que entenderemos de esta última forma la igualdad de la hipótesis nula y no nos preocuparemos más por esta cuestión 18 CAPÍTULO 2. CONTRASTES BASADOS EN LA FUNCIÓN DE DISTRIBUCIÓN proceso empírico, que se denirá a continuación. Denición 2.4. Dada una muestra de n observaciones de una variable aleatoria X con función de distribución F(x) , denimos el proceso empírico como la función aleatoria Bn(x) = √n(ˆ Fn(x)−F(x)) . Por el lema 2.8, conocemos el comportamiento asintótico de Bn para cada x . Pero esto no será suciente, y precisamos estudiar el comportamiento asintótico funcional del proceso empírico. Un primer paso para enfrentarnos a este problema es considerar el vector k -dimensional, donde cada coordenada corresponde a Bn evaluada en diferentes puntos del dominio de F , y estudiar sus propiedades. Ello da lugar al siguiente resultado: Proposición 2.5. Para todo x1, x2, ..., xk valores en el dominio de F(x) , cuando n→ ∞ se tiene: (Bn(x1), ..., Bn(xk)) d −→ (B(x1), ..., B(xk)) ∈N(0k,Σ) , donde Σ=(σij) , con σij =Cov{B(xi), B(xj)}=F(xi∧xj)−F(xi)F(xj) 5 Demostración. Es consecuencia directa del Teorema Central del Límite Multivariante, solo hay que tener en cuenta que Bn(x) = √n(ˆ Fn(x)−E[I(Xi≤x)]) y que ˆ Fn(x) = 1 nPn i=1 I(Xi≤x) . Veamos que Cov{B(xl), B(xm)}=F(xl∧xm)−F(xl)F(xm) . Por el TCLM: Cov{B(xl), B(xm)}=Cov{I(X1≤xl), I(X1≤xm)}= E[I(X1≤xl)·I(X1≤xm)] −E[I(X1≤xl)] ·E[I(X1≤xm)] . Usando que I(X1≤xl)·I(X1≤xm) = I(X1≤min{xl, xm}) , obtenemos el resultado. Conforme tomamos más puntos y aumentamos la dimensión del vector, este se vuelve una mejor aproximación de la función Bn . Pero para llegar a un TCL funcional no vale solo con dejar que k tiena a innito, sino que hacen falta más condiciones. Sin embargo, para los resultado de este trabajo será suciente con pensar un TCL funcional como el límite de un TCL multivariante. Decimos entonces que el proceso empírico Bn converge débilmente al proceso límite B , y lo denotamos por: Bnw −→ B La función límite B es lo que se conoce como proceso Gaussiano, es decir, un conjunto de variables aleatorias indexadas (en este caso por los x ), tales que toda colección nita 5 x∧y=min{x, y} 2.3. EL ESTADÍSTICO DE KOLMOGOROV-SMIRNOV 19 de dichas variables tiene una distribución normal multivariante. En particular, se trata de un proceso Gaussiano de media cero y covarianza la dada por la proposición 2.5. Para nalizar esta sección presentamos un último teorema de gran utilidad: el teorema de Mann-Wald: Teorema 2.6. Sea g una función continua. Si Bnw −→ B , entonces g(Bn)w −→ g(B) cuando n→ ∞ . 2.3. El estadístico de Kolmogorov-Smirnov En las anteriores secciones hemos recopilado los resultados más importantes referentes a la FDE y el proceso empírico. Ahora podemos introducir ya el primer estadístico de este capítulo y uno de los primeros que se emplearon para realizar un contraste de bondad de ajuste: el de Kolmogorov-Smirnov, a menudo abreviado como KS. Para contrastar la hipótesis nula H0:F=F0 frente H1:F6=F0 , el estadístico toma la forma: Dn=√n·sup x|ˆ Fn(x)−F0(x)|, (2.2) con lo que Dn es de la forma de la expresión 2.1, con d(ˆ Fn, F0) = sup x|ˆ Fn(x)−F0(x)| y c(n) = √n . Además, bajo la hipótesis nula Dn=sup x|Bn(x)| y así podremos aplicar los resultados de la sección anterior. Nótese que Dn es directamente proporcional al supremo de la diferencia entre la función de la hipótesis G y la FDE. A menudo esta diferencia también se suele escribir como ˆ Fn(G−1(p)) −p(p=G(x)) . En muchos manuales el estadístico se presenta sin el factor √n , aunque nosotros preferimos nuestra notación por poner de relieve su relación clara con el proceso empírico. Smirnov introdujo dos estadísticos íntimamente relacionados con Dn , y que se denen como sigue: D+ n=√n·sup x(ˆ Fn(x)−F0(x)) = sup x(Bn(x)) (2.3) D− n=√n·sup x(F0(x)−ˆ Fn(x)) = sup x(−Bn(x)) (2.4) Está claro que D− n mide la mayor desviación positiva entre la FDE y la F0 , mientras que D− n mide la mayor desviación negativa. Así, se emplean para los contrastes de hipótesis nula H0:F=F0 y alternativas H1:F > GF0 y H1:F < F0 respectivamente, siendo > 20 CAPÍTULO 2. CONTRASTES BASADOS EN LA FUNCIÓN DE DISTRIBUCIÓN la relación de orden estocástico 6 . Justamente, combinando ambos estadísticos, obtenemos KS: Dn=max{D+ n, D− n} . Normalmente, encontrar el supremo de una función no diferenciable requiere la evaluación de la función en muchos puntos. Sin embargo, si F es una función continua, podemos suponer que: x1< x2< ... < xn . Así, dado que ˆ Fn es una función escalonada y G monótona creciente, las expresiones (2.3) y (2.4) trivialmente se simplican a: D+ n=√n·max 1≤i≤n(i n−F0(X(i))) (2.5) D− n=√n·max 1≤i≤n(F0(X(i))−i−1 n)) (2.6) Para hallar Dn solo hace falta evaluar F0 en los n valores de la muestra. Además, esta expresión de los estadísticos permite entrever una propiedad muy importante sobre su comportamiento: su distribución nula no depende de F0 . Proposición 2.7. Dada una muestra x1, x2, ..., xn de una variable aleatoria X con una función de distribución F continua, las distribuciones nulas de los estadísticos Dn, D+ n, D− n no dependen de F . Demostración. Por las expresiones 2.5 y 2.6, si se cumple H0:F=F0 deducimos que D+ n y D− n solo dependen de las variables U(1) =F(X(1)), U(2) =F(X(2)), ..., U(n)=F(X(n)) , que son los estadísticos de orden de las variables aleatorias U1=F(X1), U2=F(X2), ..., Un= F(Xn) . Pero si F es continua, Ui=F(Xi) sigue una uniforme de intervalo [0,1] independientemente de F con lo que tenemos el resultado. Escribiendo Dn=max{D+ n, D− n}= max 1≤i≤n(i n−F(X(i)), F(X(i))−i−1 n) tenemos el resultado también para Dn . Ahora, si queremos emplear el estadístico KS para inferencia necesitamos conocer su distribución bajo la nula. Como se acaba de ver que esta es libre, este problema suele ser abordado suponiendo que F(x) es una uniforme en (0,1) y a partir de ahí se halla la distribución de cada Dn . Aunque estas se han computado hasta un n grande, el proceso es bastante tedioso y suele ser más común emplear la distribución asintótica de los Dn . Gracias a los resultados vistos en las dos primeras secciones de este capítulo podemos hallarla muy fácilmente: Lema 2.8. Cuando n→ ∞ y si H0:F=F0 es cierta, Dn=sup x|Bn(x)|d −→ sup x|B(x)|=: D , donde B es el proceso límite de Bn . Demostración. Considerando el proceso límite Bnw −→ B es consecuencia directa del teorema 2.6. 6 Dadas X e Y dos variables aleatorias con funciones de distribución F y G respectivamente, decimos que X es estocásticamente menor que Y y representamos por F > G , si para todo ∀z , Pr{X < z}< Pr{Y < z} . 2.3. EL ESTADÍSTICO DE KOLMOGOROV-SMIRNOV 21 Existe incluso una expresión analítica de la función de distribución de D que se ha comprobado sucientemente buena para n > 35 , aunque la demostración no se incluirá en este trabajo: FD(d)=1−2∞ X j=1 (−1)j+1exp(−2j2d2) (2.7) Entonces, si n > 35 el test rechazaría la hipótesis nula H0:F=F0 con un nivel de signicación α , si y solo si Dn> F−1 D(1 −α) . Como FD viene expresado en serie de potencias, habría que aproximar el valor FD(1 −α) para realizar el contraste. 2.3.1. El estadístico KS para hipótesis compuestas Hasta ahora hemos estudiado el estadístico KS en el contexto de los test simples. Sin embargo, análogamente a los estadísticos basados en la multinomial, este puede ser adaptado para realizar también contrastes de hipótesis compuestas. Aunque veremos que no hay una forma general de proceder para ello. Supongamos que la distribución de la hipótesis, F0 , es conocida hasta un vector de parámetros βt k= (β1, ..., βk) de la que es dependiente. Denotamos entonces F0(x) = F0(x;β) para explicitar dicha dependencia. La manera de proceder será sustituir la hipótesis inicial por H0:F(x) = F0(x;ˆ βn) , siendo ˆ βn un estimador del que se precisará que cumpla varias condiciones más adelante. Todos los estadísticos de este capítulo están basados en el proceso empírico, con lo que es necesario estudiar cómo se comporta cuando la distribución F0 no está completamente especicada. Denotamos Bn(x;β) = √n(ˆ Fn(x)−F0(x;β)) y B(x;β) su proceso límite bajo la hipótesis nula. Como sabemos, este tiene media cero y su función de covarianza cumple: Cov{B(x;β), B(y;β)}=F0(x∧y;β)−F0(x;β)F0(y;β) . Notemos que a pesar de escribir la función de distribución F0(x) de la forma F0(x;β) , simplemente es una cuestión de notación y la teoría de la sección anterior se aplica normalmente. El problema aparece cuando β es desconocido y ha de ser sustituido a la hora de calcular los procesos empíricos por su estimador ˆ βn , con lo que naturalmente los resultados vistos hasta ahora dejan de ser válidos. El siguiente teorema muestra la nueva distribución asintótica de ˆ Bn(x;ˆ βn) = √n(ˆ Fn(x)−F0(x;ˆ βn)) , para un tipo muy concreto de estimadores: Teorema 2.9. Si ˆ βn es un estimador lineal localmente asintótico 7 , bajo la hipótesis nula el 7 Se dice del estimador ˆ βn de β que verica la siguiente expresión: ˆ βn−β=1 nPn i=1 ψ(Xi;β)+op(n−1/2) , donde ψt= (ψ1, ..., ψp) es una función vectorial continuamente diferenciable de Rp en Rp con media cero y E{ψ(X;β)ψt(X;β)} es nito y no singular. 22 CAPÍTULO 2. CONTRASTES BASADOS EN LA FUNCIÓN DE DISTRIBUCIÓN proceso empírico estimado ˆ Bn(x) converge débilmente a un proceso Gaussiano ˆ B de media cero y función de covarianza: Cov{ˆ B(x),ˆ B(y)}= =F0(x∧y;β)−F0(x;β)F0(y;β)−ψt(x;β)h(y;β)−ψt(y;β)h(x;β) + ht(x;β)Σψh(y;β) , donde h(x;β) = ∂F0(x;β)/∂β , Ψ(x;β) = Rx −∞ ψ(z;β)dF0(z;β) y Σψ=V ar(ψ(x;β)) . Hay dos consecuencias muy importantes de la convergencia de ˆ Bn a ˆ B . La primera es que nos permite hallar la distribución nula del estadístico KS. Así, bajo H0 , cuando n→ ∞ : ˆ Dn=Dn(ˆ βn) = √n·sup x|ˆ Fn(x)−F0(x;ˆ βn)|=sup x|ˆ Bn|d −→ sup x|ˆ B(x)| (2.8) La segunda es que la distribución límite depende del β desconocido, además de la distribución F0 . Por tanto, ya no hay una forma general de realizar el contraste. Sin embargo, si F0 es una distribución invariante a cambios de escala y posición, 8 ˆ D se simplica de forma que deja de depender de β (si bien continúa dependiendo de F0 ). Un caso particular de este tipo de distribuciones es la normal, que ha sido estudiado por Lilliefors. Este propuso aplicar el estadístico KS a las variables estandarizadas Zi=Xi−X √S (siendo X y S la media y la cuasivarianza muestral). Para ello calcularíamos FDE tras la estandarización, Fn , y el estadísticos toma la expresión: ˆ Dn=√n·sup z|Fn(z)−φ(z)| , donde φ es la distribución de una normal estándar. Lilliefors fue el primero en tabular las distribuciones exactas de los ˆ Dn y luego se han empleado diferentes métodos computacionales para hallar su distribucón asintótica. Los puntos críticos del estadístico para los p -valores más empleados se pueden consultar en el manual Nonparametric Statistical Inference de Jean Dickinson Gibbons y Shubhabrata Chakabort. 2.3.2. El estadístico KS aplicado al caso discreto Hasta ahora, todos los resultados expuestos sobre las propiedades del estadístico KS se han basado en que su distribución no depende de la de la variable aleatoria de la muestra, y por tanto requieren como hipótesis la continuidad de la función G . Sin embargo, la teoría presentada en las dos primeras secciones, incluyendo el teorema de Givenko-Cantelli 2.3, no precisan de esta suposición, con lo que existe una motivación para emplear el estadístico en el caso de G discreta. También puede interesar agrupar la muestra en k clases, pero ahora querer emplear el KS y no los estadísticos del primer capítulo. Como esto es básicamente una discretización de la función de distribución F , ambos problemas son equivalentes. 8 Se dice que una distribución es invariante a transformaciones de escala y posición si admite función de densidad f0 y esta cumple: g(x;µ, σ) = g((x−µ)/σ; 0,1) 2.3. EL ESTADÍSTICO DE KOLMOGOROV-SMIRNOV 23 Empleando la notación del primer capítulo, supongamos que tenemos k clases, cada una con una probabilidad esperada p0 i dada por la función F0 , y Ni representa el número de observaciones de cada clase de la muestra. Entonces tenemos que, dado j∈ {1, ..., k} : ˆ Fn(xj) = j X i=1 Ni/n, G(xj) = j X i=1 p0 i, (2.9) con lo que el estadístico KS toma la forma: Dn=max i≤j≤k j X i=1 Ni−np0 i √n=max i≤j≤k j X i=1 Oi−Ei √n (2.10) De esta forma, la distribución de Dn depende de cómo fueron establecidas las clases a la hora de agrupar los datos. La distribución exacta de √nDn ha sido tabulada (ver Pettitt and Stephens ). Es especialmente signicativa la similitud entre el KS para el caso discreto y el estadístico de Pearson X2 . Justamente, para terminar esta sección, haremos una pequeña comparación entre ambos estadísticos. Una primera diferencia obvia es que el estadístico de Pearson requiere que los datos estén agrupados en clases, al contrario que el KS. Por tanto, para distribuciones continuas el KS tiene un uso más completo de la muestra, inuyendo cada dato de forma propia, y su uso supone no enfrentarse al problema del la elección del número de clases y de las fronteras de cada una. Además, la distribución exacta de los Dn es conocida y ha sido tabulada para todo n , mientras que a la hora de trabajar con los estadísticos de divergencia solo podemos basarnos en sus distribuciones asintóticas a la χ2 , que supone una buena aproximación siempre y cuando tengamos sucientes datos y las frecuencias esperadas sean lo bastante grandes. Sin embargo, los estadísticos de divergencia tienen la ventaja de tener una distribución asintótica conocida cuando hay presentes parámetros desconocidos (recordamos que cada parámetro resta un grado de libertad de la χ2 límite), mientras que ˆ Dn tiene una distribución diferente a Dn , para la cual no contamos con una fórmula general. Para nalizar, cabe decir que el test del KS tiene mayor poder que el de la chi-cuadrado de Pearson, tanto para datos categorizados como para distribuciones continuas. Con lo que se puede concluir que, exceptuando los casos de hipótesis compuestas para las que no se conozca la distribución límite del ˆ Dn , el estadístico KS es el más adecuado de los dos. 24 CAPÍTULO 2. CONTRASTES BASADOS EN LA FUNCIÓN DE DISTRIBUCIÓN 2.4. Estadísticos de Anderson-Darling Hasta ahora hemos estudiado con detenimiento el estadístico KS y su comportamiento a la hora de contrastar hipótesis simples o compuestas, y para distribuciones continuas y discretas. Sin embargo, este solo es un ejemplo del conjunto de estadísticos construibles siguiendo la fórmula 2.1. Otra familia muy importante de estadísticos es la introducida por Anderson-Darling: Tn=ZS w(G(x))B2 n(x)dG(x), (2.11) donde w(.) actúa como una función de peso. Cuando w(u)=1 para todo 0≤u≤1 , el estadístico suele ser conocido como el estadístico de Cramér-von Mises, que abreviaremos por CvM. Además, aunque hemos presentado una colección de estadísticos indexada por una función de peso, en la práctica solo hay un único estadístico de Anderson-Darling que es corrientemente empleado (aparte del CvM): el que tiene w(x) = 1/(u(1 −u)) como función de peso, y que por tanto será al que nos reramos cuando hablemos de estadístico de Anderson-Darling y abreviaremos por AD. La razón por la cual se emplea esa función de peso particular es que tiene un efecto estabilizador en la varianza: pw(u)Bn(u) tiene varianza constante igual a 1 . A pesar de que calcular el valor de los estadísticos pueda parecer complicado, para CvM y AD existen expresiones explícitas simples. Denotémolos de forma respectiva por Wn y An para una muestra de tamaño n . Sean Ui=G(Xi) y U(i) la estadística de orden i -ésimo de las variables U1, ..., Un . Entonces: An=−n−1 n n X i=1 (2i−1)(log U(i)+log(1 −U(n+1−i))) (2.12) Wn=1 12n+ n X i=1 U(i)−2i−1 2n2 (2.13) Cuando el contraste involucra una hipótesis compuesta: H0:F(x) = G(x;β) se procede de manera similar al test KS. Primero, calcularíamos el estimador de β , ˆ βn , que suponemos asintóticamente lineal, y sustituimos G(x;β) por G(x;ˆ βn) . Menos es conocido sobre el comportamiento exacto y asintótico de los estadísticos CvM y AD que sobre el del KS. Para todo estadístico de Anderson-Darling Tn ( ˆ Tn cuando trabajamos con hipótesis compuestas), los resultados vistos en la sección anterior siguen siendo válidos, a saber: la distribución de Tn bajo la nula no depende de G y usando 2.4. ESTADÍSTICOS DE ANDERSON-DARLING 25 de nuevo la convergencia débil del proceso empírico Bn podemos extender el lema 2.8 al siguiente teorema: Teorema 2.10. Si R1 0w(t)B2(t)dt < ∞ y ˆ βn es un estimador asintóticamente lineal de β , entonces, bajo la hipótesis nula simple y compuesta respectivamente: Tnd →Z1 0 w(t)B2(t)dt (2.14) ˆ Tnd →Z1 0 w(t)ˆ B2(t)dt, (2.15) cuando n→ ∞ . Ninguna de estas fórmulas es verdaderamente útil para hallar directamente los puntos críticos de los estadísticos. En el caso de la hipótesis compuesta recordemos que ˆ B depende de G y de β , con lo que en general no se pueden tabular los p-valores. Aunque sí que es posible hacerlo cuando G es invariante a cambios de locación y escala, como en el caso de que G sea una normal o una exponencial. En el caso de la hipótesis simple se han podido obtener varios resultados interesantes. A partir de sus funciones características, Anderson y Darling hallaron expresiones de las distribuciones asintótitcas de Wn y An , W y A respectivamente, como sumas innitas de variables aleatorias: W=∞ X j=1 1 j2π2Z2 j (2.16) A=∞ X j=1 1 j(j+ 1)Z2 j, (2.17) donde Z1, Z2, ... son variables aleatorias i.i.d siguiendo una normal estándar. En cuanto a las distribuciones exactas de los estadísticos Wn y An , a pesar de haber sido largamente estudiadas durante las últimas décadas, los resultados no son tan satisfactorios como los del KS. El estadístico CvM solo ha sido tabulado para n= 1, ..., 7 , y el AD para un escueto n= 1 . Sin embargo, se ha visto que la distribución asintótica supone una buena aproximación para valores de n tan pequeños como n= 8 . El test funcionaría de manera similar al del estadístico KS. 2.4.1. Comentarios nales Utilizando técnicas de simulación, se ha llegado a la conclusión de que tanto el estadístico CvM como el AD tienen una muy buena potencia en comparación con muchos otros 32 CAPÍTULO 3. CONTRASTES BASADOS EN LA FUNCIÓN DE DENSIDAD construir medidas del ajuste global. Queremos ver cómo se comporta nuestro estimador en la recta real y no sólo en un punto indeterminado. Es decir, buscamos una distancia funcional, y precisamente emplearemos el cuadrado de la inducida por el producto escalar de L2 y que da lugar a un medidor muy empleado en estadística paramétrica: el Error Cuadrático Integrado. Este viene dado por la siguiente fórmula: ECI(ˆ fn(x;h)) = Z[ˆ fn(x;h)−f(x)]2dx. (3.8) Este medidor es de utilidad si solo nos preocupa la muestra que tenemos. Para estudiar el comportamiento de nuestro estimador no para un ejemplo concreto, sino para una muestra cualquiera de nuestra densidad f , tendríamos que calcular su esperanza, lo que da lugar al Error Cuadrático Integrado Medio (ECIM), ECIM(ˆ fn(x;h)) = E[ECI(ˆ fn(x;h))] = EZ(ˆ fn(x;h)−f(x)2dx (3.9) Haciendo un cambio en el orden de integración obtenemos una fórmula en fuunción del ECM: ECIM(ˆ fn(x;h)) = ZE(ˆ fn(x;h)−f(x))2=ZECM(ˆ fn(x;h)) (3.10) Y aplicando 3.7 obtenemos una expresión muy util del ECIM: ECIM(ˆ fn(x;h)) = n−1Z[(K2 h∗f)(x)−(Kh∗f)2(x)]dx+ +Z[(Kh∗f)(x)−f(x)]2dx (3.11) Podemos notar aquí una consecuencia muy importante del sesgo de fn(x;h) : el ECI no es cero ni tiene a cero conforme n se hace grande; lo cual supondrá un problema más tarde para nuestro estadístico de contraste. Para terminar esta sección, daremos una expresión alternativa de ECIM que proviene de emplear la extensión de f(x−hz) como serie de Taylor y realizar unos cálculos no muy complicados que pueden ser consultados en el manual Kernel Smoothing de Wand y Jones. Así, siempre que tenga sentido, tenemos la igualdad siguiente: ECIM(ˆ fn(x;h)) = (nh)−1ZK(x)2dx+1 4h4Zz2K(z)dz2Zf00(x)2dx+o[(nh)−1+h4]. 5 (3.12) 3.2. CONTRASTE DE HIPÓTESIS 33 A la suma de dos primeros términos de esta expresión se la denomina ECIM asintótico (ECIMA): ECIMA(ˆ fn(x;h)) = (nh)−1ZK(x)2dx +1 4h4Zz2K(z)dz2Zf00(x)2dx, (3.13) y comportamiento en función de K y h es clave en la búsqueda de estimadores de densidad kernel óptimos. 3.2. Contraste de hipótesis 3.2.1. Estadísticos de contraste Hasta aquí lo que hemos hecho ha sido presentar la fórmula general para el estimador kernel de densidad y estudiar sus propiedades más importantes. Ahora emplearemos los resultados anteriores para motivar el uso de los dos estadísticos propuestos por Bickel y Rosenblatt y analizarlos para poder entender sus limitaciones. Como ya se ha comentado en la introducción de este capítulo, el estudio de los estimadores de densidad ha sido históricamente independiente del de los contrastes de especicación, a pesar de su potencial claro como herramienta para este tipo de test. Fueron Bickel y Rosenblatt en su artículo de 1973 quienes por primera vez proponen el uso de estadísticos basados en estimadores de densidad kernel y estudian sus propiedades asintóticas. Lo más interesante es la intuición que tienen de traducir el test de la χ2 de Pearson (recordamos que estaba basado en el estadístico: X2=Pn i=1 (Oi−Ei)2 Ei ) a este nuevo contexto. Así, donde Pearson usa observaciones Oi y valores esperados Ei , Bickel y Rosenblatt los sustituirán por el estimador de densidad ˆ fn(x;h) y la densidad hipotética f0(x) . Y como estamos trabajando con variables aleatorias continuas, el sumatorio se traduce en una integral, y en vez de dividir por los valores esperados, añaden una función de peso w(x) a la integral. Es decir, la fórmula general del primer estadístico de Bickel y Rossenblet es: D0 n=Z(ˆ fn(x;h)−f0(x))2w(x)dx (3.14) En este trabajo tomaremos siempre w(x) = 1 y trabajaremos entonces con el estadístico: 5 Se emplea la notación an=o(bn) cuando n→ ∞ , si y solo si lim n→∞|an/bn|= 0 . Se escribe an=O(bn) cuando n→ ∞ , si y solo si lim n→∞|an/bn|<∞ . 34 CAPÍTULO 3. CONTRASTES BASADOS EN LA FUNCIÓN DE DENSIDAD T0 n=Z(ˆ fn(x;h)−f0(x))2dx (3.15) Decidimos esto así por dos motivos. El primero es para evitar extendernos en la cuestión de la elección de w(x) y su efecto en el test, aunque todos los resultados que se presenten más adelante pueden ser generalizados para estadísticos de la forma 3.14. Y el segundo es para poder basarnos en el análisis paramétrico anterior de nuestro estimador de densidad a la hora de estudiar las propiedades de nuestro estadístico T0 n , que bajo la hipótesis nula es exactamente el ICE de nuestro estimador. Así, todos los resultados vistos en la sección anterior se aplican. Naturalmente, esto implica que T0 n arrastra también la problemática del sesgo de ˆ fn(x;h) . No debemos olvidar que lo que buscamos es una distancia entre funciones, d(., .) , tal que H0:f=f0⇔d(f, f0)=0 . En este sentido, d(ˆ fn, f0) es un estimador de d(f, f0) , y por tanto nos interesa que tenga buenas propiedades, como que sea insesgado o consistente. La fórmula 3.11 nos da una expresión de la esperanza del estadístico bajo la hipótesis nula: E(T0 n) = n−1Z[(K2 h∗f0)(x)−(Kh∗f0)2(x)]dx +Z[(Kh∗f0)(x)−f0](x)]2dx (3.16) Vemos claramente que la esperanza de T0 n es estrictamente mayor que 0 , y que, de hecho, no tiene a 0 cuando n se hace más grande. Es decir, T0 n no es un estimador insesgado ni asintóticamente insesgado de d(f, f0) bajo la hipótesis nula, lo cual es consecuencia del sesgo de ˆ fn . Existen dos formas de solventar este defecto. La primera y más obvia es sustituir el estadístico T0 n por uno que tenga mejores propiedades; en particular, que sea asintóticamente insesgado. Como el problema siempre es la diferencia entre el producto de convolución de nuestra función hipotética y la propia función, la solución más inmediata consiste en sustituir f0 por Kh∗f0 en la denición de nuestra T0 n , lo que da lugar al nuevo estadístico Tn con fórmula: Tn=Z[ˆ fn(x;h)−(Kh∗f0(x))]2dx, (3.17) que también podemos escribir como: Tn=Z[ˆ fn(x;h)−EH0(ˆ fn(x;h))]2dx, (3.18) donde EH0(ˆ fn(x;h)) indica la esperanza de nuestro estimador de densidad bajo la hipótesis nula. Este estadístico es propuesto de una forma más general por Bickel y Rosenblatt 3.2. CONTRASTE DE HIPÓTESIS 35 al incluir también una función de peso, pero como ya hemos indicado preferiremos tomar siempre como peso la unidad y así obviar el problema de su elección. La segunda solución consiste en hacer tender nuestro ancho de banda a 0 conforme aumenta el tamaño de la muestra. Efectivamente, si en la expresión 3.12 sustituimos la f por nuestra f0 tendremos la otra expresión de la media de T0 n bajo la hipótesis nula, siempre y cuando se cumplan ciertas condiciones de regularidad. E(T0 n)=(nh)−1ZK(x)2dx +1 4h4Zz2K(z)dz2Zf00(x)2dx +o[(nh)−1+h4]. (3.19) Lo interesante de esta expresión es que depende de nuestra h de una forma muy simple y nos permite ver que si conseguimos hacer tender tanto (nh)−1 como h a 0 , tendremos que T0 n sí será asintóticamente insesgado. Además, este expresión nos ayudará a entender el motivo de algunas de las restricciones que precisaremos para poder asegurar la convergencia de ˆ Tn que estudiaremos ahora. Lo que podemos adelantar es que necesitaremos tomar una sucesión de anchos de banda hn (por comodidad seguiremos empleando h para denotar h=hn ) tal que lim n→∞h= 0 y lim n→∞(nh)−1= 0 . 3.2.2. Comportamiento asintótico Hasta ahora hemos presentado la fórmula general de un estimador de densidad kernel, hemos empleado herrramientas del análisis paramétrico para llegar a la expresión de los dos estadísticos que estudiaremos en este capítulo y hemos obtenido ciertas restricciones que nos aseguran un mejor comportamiento de al menos uno de ellos. Hemos intentado que todo ello fuese de forma uida y, al mismo tiempo, todo lo riguroso que hemos sabido. Pero este enfoque tiene que ser abandonado en este punto. Queda por ver cuáles son las distribuciones asintóticas de nuestros estadísticos de contraste bajo la nula, de forma que nos sirvan para hacer inferencia. El problema es que las demostraciones toman una complejidad bastante elevada, tanto por su tamaño como por la base teórica que requieren, y por ello serán pasadas por alto. Solo daremos entonces las condiciones que necesitaremos para asegurar sus convergencias (muchas de ellas están estrechamente relacionadas con la fórmula 3.19 como se podrá ver), y sus propias distribuciones límite. Añadamos entonces las siguientes restricciones a las ya dadas: 1. Nuestro ancho de banda es una sucesión h=hn tal que lim n→∞h= 0 y n−1=o(h) . 2. K cumple o que (a) es nula fuera de un intervalo [−A, A] y absolutamente continua 36 CAPÍTULO 3. CONTRASTES BASADOS EN LA FUNCIÓN DE DENSIDAD en [−A, A] , o bien que (b) es absolutamente continua en la recta real y R|K00|k<∞ para k= 1,2 . 3. La densidad f es continua, positiva y acotada. Y f1/2 es absolutamente continua y su derivada 1 2f0/f1/2 está acotada. Y además: Z[|z|3/2≥3] |z|3/2[log log|z|]1/2[|K0(z)] + |K(z)|]dz < ∞ (3.20) 4. La segunda derivada f00 de f existe y está acotada. 5. Existe Rz2K(z)dz Una vez expuestas estas restricciones podemos comenzar por el estadístico Tn . Denimos µ(K) = 1 nh RK2(z)dz y σ(K)2= 2[R(RK(u+v)K(v)dv)2du]Rf2 0(x)dx . Se cumple entonces el siguiente teorema: Teorema 3.1. Supngamos que se cumplen las restricciones 1-3, entonces, bajo la hipótesis nula y cuando n→ ∞ : n√hTn−µ(K) σ(K) d −→ N(0,1) (3.21) Para presentar este resultado hemos hecho uso de solo las tres primeras restricciones que hemos presentado, lo cual tiene sentido ya que Tn es un estadístico con mejores propiedades que T0 n . Nótese además en 3.19, que necesitamos las restricciones 4 y 5 para que ˆ Tn tenga esperanza. Sin embargo, el estadístico T0 n es probablemente de mayor interés que Tn , y gracias a las restricciones 4 y 5, podemos asegurar también su convergencia a una normal de la que indicaremos la media y la varianza en el siguiente teorema. Teorema 3.2. Supongamos que se cumplen 1-5. Entonces, bajo la hipótesis nula y cuando n→ ∞ : n√hT0 n−µ(k) σ(k) d −→ N(0,1) (3.22) Es decir, siempre que se cumplan las restricciones 1-5, los dos estadísticos tienen el mismo comportamiento asintótico. El contraste de la hipótesis simple H0:f=f0 , con un nivel de signicación α lo haríamos de la siguiente forma: hallamos el valor de Tn a partir de la muestra, calculamos d(α) = µ(K)+(n√h)−1σ(K)φ(1 −α) y, si Tn> d(α) , rechazamos la hipótesis nula. En caso contrario la aceptamos. 3.3. COMENTARIOS FINALES 37 3.2.3. Hipótesis compuesta Si el uso del estimador de densidad kernel para el contraste de especicación de una hipótesis simple es de mediados de la década de los ochenta, su estudio para el contraste para una hipótesis compuesta H0:f∈ {fθ:θ∈Θ} es aún más reciente. En efecto, Bickel y Rossenblet no estudian el comportamiento de los estadísticos T0 n y Tn cuando sustituimos nuestra densidad f0 por una fˆ θ con ˆ θ un estimador de θ obtenido a partir de la muestra. Ya hemos visto en el primer capítulo (a raíz de la sonada controversia entre Pearson y Fisher) que las distribuciones límites no tienen por qué ser en general las mismas, pues la inclusión del parámetro estimado dependiente de la muestra puede tener un efecto sobre su comportamiento asintótico. De nuevo, vamos a evitar entrar en demostraciones demasiado laboriosas e iremos directamente a los resultados que nos interesan. Solo vamos a presentar un teorema sacado de Fan (1994), quien propone emplear el estadístico ˆ Tn que resulta de sustituir f0 por fˆ θ en la fórmula de Tn , y donde ˆ θ lo hallamos empleando el método de máxima verosimilitud. Es decir: ˆ Tn=Z(ˆ fn(x;h)−(Kh∗fˆ θ)(x))2dx (3.23) Bajo ciertas hipótesis de regularidad que no introduciremos (se puede consultar en Fan (1994)) y bajo la hipótesis nula, tenemos que, cuando n→ ∞ : n√hˆ Tn−µ(k) σ(K, ˆ fn)) d −→ N(0,1), (3.24) donde ahora σ(K, ˆ fn)2= 2 "ZZK(u+v)K(v)dv2 du#Zˆ fn(x;h)2(x)dx El test resultante rechaza la hipótesis nula H0:f∈fθ:θ∈Θ cuando ˆ Tn> d(α) con d(α) = µ(K)+(n√h)−1σ(K)φ−1(1 −α) . 3.3. Comentarios nales Para nalizar este capítulo vamos a comentar ciertas cuestiones importantes que hay que tener cuenta a la hora de implementar los test en la práctica. En primer lugar, apenas se ha hablado de la inuencia de la función kernel elegida para el contraste, aunque sí que hemos adelantado que el comportamiento del estimador 38 CAPÍTULO 3. CONTRASTES BASADOS EN LA FUNCIÓN DE DENSIDAD es bastante similar independientemente de la K escogida. Se ha demostrado que el kernel de Epanechnikov, que denotamos por K∗ , es el que minimiza el ECIMA y suele ser usado como referencia para estudiar el comportamiento de otros kernels propuestos. Decimos que el kernel K tiene una eciencia el 90 % si el estimador de densidad óptimo (el que emplea K∗ como kernel) puede conseguir el mismo ECIMA utilizando el 95 % de los datos que usando K . Un estudio de la eciencia de distintas funciones muestra que los kernels clásicos, entre los que se encuentra la normal o la uniforme, tienen un rendimiento superior al 90 % . Por ello la elección suele atender a criterios computacionales, y kernels como el uniforme (que es constante a trozos) o incluso el de Epanechnikov (cuya derivada es discontinua) suelen ser rechazados por otros más regulares. En segundo lugar, queda el problema de la elección del ancho de banda. Estudios por simulación de los test basados en estimadores de densidad kernel han mostrado una gran dependencia del smoothing parameter , por lo que es una cuestión clave. Sin embargo, también es una cuestión de muchísima complejidad. Se han propuesto numerosos criterios para denir una expresión de h que resulte en buenas propiedades para nuestro estimador (por suepuesto, todas estas expresiones cumplen que h tiende a cero conforme el tamaño de la muestra se hace grande). Sin embargo, su éxito en la reducción del ECIMA depende en gran medida de la función de densidad que siga la muestra, con lo que es difícil establecer un criterio general. Además, como estas cuestiones nos interesan desde el punto de vista de los contrastes de especicación, existe la dicultad añadida del hecho de que muchos de estos criterios hagan depender a h no solo de n , sino de los propios valores de la muestra. Esto crearía una nueva dependencia de nuestros estadísticos sobre la muestra, y no tenemos la seguridad de la validez de los teoremas presentados sobre su comportamiento asintótico. De esta forma, la elección óptima del parámetro es una cuestión que continúa hoy en día sin tener una respuesta clara. Capítulo 4 Ilustración en bases de datos simulados Para nalizar el trabajo, vamos a llevar a la práctica los test de contraste vistos en los tres capítulos anteriores. Usaremos para ello bases de datos simulados, lo que nos permitirá poder interpretar el resultado de los test sabiendo en todo caso si nuestra hipótesis nula es cierta o no. Es decir, sabremos si el test es acertado o está cometiendo algún tipo de error. Es necesario aclarar que en ningún caso pretendemos hacer un estudio pormenorizado y sistemático de los estadísticos de contraste. No buscamos concluir qué estadístico es el más adecuado, qué número de intervalos es el óptimo para el test de Pearson o cuál es el ancho de banda idóneo para nuestro estimador de densidad. Para ver estudios verdaderamente rigurosos a partir de simulación nos remitimos a la bibliografía. En este trabajo hemos comentado muchas de las conclusiones de estos estudios (así como cuestiones que están lejos de estar concluidas), y no pretendemos verdaderamente añadir ni modicar nada. Nuestro objetivo es otro, y, además, es doble. En primer lugar, y especialmente, buscamos ejemplicar todos los métodos de los que hemos hablado, de forma que su comprensión teórica se complete con un caso aplicado a una muestra concreta. De poco sirve comprender la teoría que lleva a la identicación del proceso empírico en un proceso Gaussiano, si no sabemos cómo realizar un contraste empleando el estadístico de Kolmogorov-Smirnov. Este trabajo trata en último lugar sobre contrastes, y aunque nos hayamos centrado en los aspectos más teóricos, pues nos interesaba sobretodo hacer una revisión de las familias de estadísticos más importantes y de las herramientas matemáticas en las que se basan, no podemos olvidar el problema no tan obvio de cómo aplicar la teoría en un caso práctico. Por otro lado, este capítulo nos servirá para seguir introduciendo conceptos y técnicas de la estadística computacional tan fundamentales como el método Monte Carlo o el 39 40 CAPÍTULO 4. ILUSTRACIÓN EN BASES DE DATOS SIMULADOS Bootstrap paramétrico. Nuestro acercamiento a estos métodos será sobretodo intuitivo y a partir de ejemplos. Los iremos motivando a través de las diferentes cuestiones que nos irán surgiendo sobre las propiedades de nuestros estadísticos al trabajar con nuestros datos simulados. Un enfoque más riguroso ocuparía demasiado espacio y haría que nos desviáramos excesivamente del tema que nos ocupa. De ahora en adelante daré una descripción de los métodos implementados en R para el contraste e incluiré los resultados obtenidos seguidos de las conclusiones que podemos sacar. El código lo presentaremos al nal del trabajo en un apéndice. 4.1. Hipótesis simple Vamos a simular en R un muestra Sn de una normal estándar tamaño n que de momento no jaremos pues nos interesará variarla en ocasiones. Nuestra hipótesis nula será correcta, H0:Sn∈N(0,1) frente a H1:Sn/∈N(0,1) , con lo que podremos estudiar el error de tipo I de nuestros test. 4.1.1. Contraste usando los estadísticos de divergencia Comenzamos por los test del primer capítulo. La primera cuestión es cómo discretizar la normal. Como estamos en el caso de la hipótesis simple, podemos aplicar el criterio de equiprobabilidad. Además vamos a tomar siempre n≥50 , con lo que podremos tomar k= 10 intervalos equiprobables y de forma que las observaciones esperadas en cada intervalo sean siempre mayores o iguales que 5 , como recomendábamos para que la χ2 sea una buena aproximación de la distribución exacta de nuestros estadísticos. Vamos a trabajar con tres estadísticos: el X2 de Pearson, el G2 del test de razón de verosimilitudes y el 2nIλ con λ= 2/3 , que es el sugerido como óptimo por Read y Cressie de entre todos los estadísticos de la familia de divergencia y que denominamos PD. Como tienen expresiones analíticas podemos calcularlos con facilidad (existen funciones de paquetes de R pero solo las usaremos para comparar). Usaremos siempre una seed determinada de forma que nuestros resultados puedan ser más tarde corroborados. Fijaremos siempre set.seed (3000) y tomamos n= 50 . Para hallar el p-valor utilizamos la distribución asintótica χ2 9 : para el estadístico X2 lo calcularíamos como 1−Fχ2 9(X2) siendo Fχ2 9(X2) la función de distribución de la χ2 9 . Para G2 y PD se haría de forma similar. Los resultados los recogemos en la tabla 4.1. Tomando cualquier nivel de signicación clásico y aplicando cualquiera de los test aceptaríamos la hipótesis nula. Nótese que los tres valores son muy similares (como hemos probado). Este sería el procedimiento que deberíamos para realizar el contraste. Por 4.1. HIPÓTESIS SIMPLE 41 Estadístico X2G2PD Valor 12,000 11,860 11,751 p-valor 0,213 0,221 0,228 Cuadro 4.1: Tabla con el valor de los estadísticos supuesto, en R ya existen funciones incluidas en paquetes que realizan estos test, y que recomendamos utilizar. Se puede comprobar sin embargo que los resultados obtenidos son exactamente los mismos. Podemos sin embargo intentar ir un paso más allá. Este resultado concreto parecería indicar que los tres test funcionan correctamente a la hora de reconocer la función hipotética, pero una sola muestra no parece suciente para poder sacar conclusiones sobre el error de tipo I. Una idea interesante sería repetir este proceso un número muy grande de veces y comparar cuántas veces la hipótesis nula es rechazada por cada tipo de test respecto del total. Necesariamente tenemos que olvidarnos de jar la seed , pues queremos muestras que puedan ser diferentes sacadas de la normal estándar. Lo haremos para 100000 repeticiones y haciendo variar la n . Surge aquí un problema con el estadístico de razón de verosimilitudes: no tiene valor cuando algún intervalo no tiene valores observados. Cuando n es grande este inconveniente desaparece pues siempre habrá alguna observación en alguno de los 10 intervalos, pero cuando n= 50 este problema se da de forma ocasional y si usamos varios miles de muestras va a haber sin duda casos en los que no esté denido. Lo ideal sería poder elegir intervalos equiprobables no vacíos, pero hacerlos depender de la muestra es arriesgado como ya hemos comentado. Para no entrar en mayores complicaciones y como este problema solo existe si n es bajo, vamos a calcular de igual forma cuántas veces el test da positivo de todas las veces en las que esté denido, e indicaremos en una nota a pié de página cuántas veces no se ha podido denir. Lo haremos para n∈ {50,100,500,1000} . Los resultados están recogidos en la tabla 4.2. Como era de esperar, el número de veces que los estadísticos han aceptado el test es similar, y parece tender a un 95 % del total. Esto es se explica porque el nivel de signicación que hemos usado es α= 0,05 . Los test dan positivo siempre que el estadístico de contraste es inferior al cuantil de orden 95 de la χ2 9 , con lo que estamos comprobando que este valor es muy cercano al cuantil de orden 95 de las distribuciones exactas de nuestros estadísticos según cada n . Esto es lógico, pues por la teoría asintótica vista sabemos que 48 CAPÍTULO 4. ILUSTRACIÓN EN BASES DE DATOS SIMULADOS Para analizar el error de tipo I de nuestros estadísticos usamos Monte Carlo: simulamos muchas muestras y hallamos el error de tipo empírico. Nótese que este método es probablemente el más costoso computacionalmente de todos, no tanto por el número de veces que tenemos que realizar aproximaciones a integrales, como por la necesidad de hallar el estimador kernel de cada muestra. Como nos interesa que sea una función y no un vector, hemos empleado el comando kdensity del paquete kdensity , pero esto conlleva un coste alto. Vamos a repetir el método 500 veces, para lo que necesitamos ya bastante potencia dependiendo de qué n tomemos. Lo haremos de nuevo para n= 50,100,500,1000 . Los resultados los recogemos en la tabla 4.9. Estadístico T0T n=50 0,926 0,994 n=100 0,95 0,994 n=500 0,974 0,98 n=1000 0,97 0,972 Cuadro 4.9: Errores de tipo I empíricos Los resultados son ligeramente insatisfactorios, pues aunque el error empírico es bastante cercano al teórico, no parece que conforme aumente el tamaño de la muestra este vaya acercándose a 0,95 , como es esperable teniendo en cuenta la convergencia en distribución de los estadísticos. También hemos realizado el test de KS para testear normalidad y siempre hemos obtenido p-valores muy bajos, con lo que rechazaríamos en todo caso que los estadísticos pudieran ser aproximados correctamente por su distribución asintótica. Por tener una idea visual de cómo de cerca se encuentra la distribución exacta de la hipotética, hemos computado la distribución empírica y la hemos representado grácamente junto con la nula para n= 1000 . Este comportamiento ligeramente alejado del teórico quizás se deba a las limitaciones con las que nos encontramos para poder realizar los test, como el tamaño modesto de la muestra o el número de repeticiones. Seguramente también tenga que ver la inuencia del ancho de banda, que ya habíamos adelantado que juega un papel determinante en el comportamiento de los estadísticos, pero para el que no contamos con criterios para optimizarlo. Teniendo en cuenta estas cuestiones, los resultados que hemos obtenido pueden darse por satisfactorios dentro de las aspiraciones de nuestro trabajo. 4.2. HIPÓTESIS COMPUESTA 49 (a) n= 100 (b) n= 1000 Figura 4.3: Ambas distribuciones parecen bastante cercanas a la normal. Sin embargo, teniendo en cuenta que estamos ante un valor de n grande, el test de kolmogorv-smirnov es muy sensible a desviaciones de lo esperado, y por tanto los p-valores son muy pequeños. Por otro lado, el test sí que ha aceptado ambos estadísticos siguen la misma distribución exacta con un p-valor muy alto (p=0.82) 4.2. Hipótesis compuesta Vamos a ejemplicar ahora el contraste suponiendo que tenemos una hipótesis compuesta. De nuevo generaremos datos de una normal estándar (elegimos de nuevo la normal porque ello nos permitirá emplear el test de Lilliefors del segundo capítulo), y consideraremos la hipótesis nula correcta, en este caso: Sn∈ {N(µ, σ) : µ∈R, σ ∈R+} . De esta forma podremos estimar de nuevo el error empírico de tipo I de nuestros test. Todos los test de nuestros capítulo precisan la estimación de los parámetros de la normal a partir de la muestra para poder llevarse a cabo. Vamos a emplear el método de razón de verosimilitudes para realizar esta estimación, ya que es el más universalmente extendido y el único para el que tenemos la garantía de todos los resultados asintóticos de nuestros estadísticos para la hipótesis compuesta. Además, contamos con la función mle del paquete stats4 para poder realizar la estimación. Recordamos sin embargo que este no es el único estimador que existe, ni el único con buenas propiedades (recuérdese la familia de estimadores BAN a la que pertenecían los estimadores asociados a los estadísticos de divergencia), pero varios de los teoremas que hemos presentado aseguran la convergencia solo para el estimador máximo verosímil. De esta forma, jemos la semilla set.seed(3000) con n= 100 y generamos la muestra. Utilizamos la función mle y obtenemos ˆµ=−0,033 y ˆσ= 0,924 , efectivamente los parámetros estimados están muy cerca de los reales. Ahora, en vez de proceder como en la sección anterior y separar los contraste por capítulos, como muchas de las cuestiones relativas a los contrastes ya han sido tratadas vamos a hacer una sección más compacta para poder 50 CAPÍTULO 4. ILUSTRACIÓN EN BASES DE DATOS SIMULADOS ya centrarnos en la comparación entre los test, aunque comentando igualmente las dudas que puedan surgir sobre los métodos aplicados al caso compuesto. Comenzando por los estadísticos de divergencia, la principal diferencia en el método respecto a la anterior sección es que ya no podemos tomar intervalos equiprobables (o por lo menos no sin que esta elección dependa de la muestra). El único criterio que podíamos seguir entonces era situar las fronteras de los intervalos de forma aleatoria. Como necesitamos que haya por lo menos un dato por intervalo, usaremos la distribución uniforme para obtener nueve puntos entre el máximo y el mínimo de nuestra muestra. Obviamente esto hace depender nuestros intervalos de los datos, pero utilizar un método aleatorio cualquiera probablemente inhabilitaría totalmente nuestro test ya que clasicaría todos los datos en un mismo intervalo. Tomaremos de nuevo k= 10 y, salvo estimación de las probabilidades de cada y intervalo y la sustitución de la χ2 9 por la nueva χ2 7 , el contraste se realiza de exacta igual forma. Siguiendo con los test del segundo capítulo nos encontramos con dos opciones. La primera consiste en aprovechar que nuestra hipótesis nula supone la normalidad de la muestra para aplicar el test de Lilliefors. Para ello estandarizamos nuestros datos empleando la media y cuasivarianza muestrales, y calculamos el estadístico de KS usando como función de distribución hipotética la de la normal estándar. Para este método contamos con tablas donde se recogen los valores críticos asociados a cada nivel de signicación para cada n . El p-valor se calcula como: 0,886/√n , para α= 0,05 y cada valor de n . La segunda posibilidad consistiría en realizar el mismo procedimiento pero empleando en vez los estadísticos de CvM y AD. En teoría esto sería posible debido a que como la normal es invariante a cambio de locación y escala, si se estandariza todos los estadísticos siguen un distribución ja (la de la normal estándar). El problema es que en la bibliografía seguida no se muestran tablas con sus p-valores más importantes, ni tampoco hemos conseguido una fuente able en internet. La solución que podemos tomar es emplear MonteCarlo para hallar el cuantil empírico de orden 95 y utilizarlo luego para realizar el contraste. Hemos obtenido 0,127 y 0,759 para CvM y AD respectivamente. Finalmente, del tercer capítulo solo tenemos un test, el propuesto por Fan (1994). Su cálculo no tiene ninguna dicultad añadida a la primera sección una vez que tenemos los estimadores de la media y la varianza, después podemos emplearlos para especicar la función de densidad en el producto de convolución y hallar el valor del estadístico. También es necesario modicar ligeramente la fórmula de la varianza asintótica. Ahora podemos intentar estudiar el error empírico de nuestros test y contrastar si nuestros estadísiticos se acercan a la distribución asintótica que les corresponde. Para los test del segundo capítulo no tenemos distribuciones asintóticas y dado que el punto crítico 4.2. HIPÓTESIS COMPUESTA 51 Test ˆ X2ˆ G2ˆ PD Lilliefors CvM AD T Valor 2,54 2,46 2,50 0,526 0,0424 0,269 0,00115 Resultado contraste Positivo Positivo Positivo Positivo Positivo Positivo Positivo Cuadro 4.10: Tabla con los resultados de los test para la muestra lo hemos simulado nosotros para los estadísticos de Anderson-Darling, solo tendría interés comprobar que la fórmula dada por Lilliefords efectivamente corresponde al cuantil de orden 95. De cualquier forma haremos el test para los tres, y así veremos que los cuantiles calculados por MonteCarlos son correctos. Para el resto de test tenemos que estimar los parámetros, lo cual conlleva un proceso de optimización que es bastante costoso. Solo vamos a correr los test para n= 100 con mil repeticiones (tabla 4.11). Test ˆ X2ˆ G2ˆ PD Lilliefors CvM AD T Error tipo I emp. 0,88 0,96 4 0,91 0,954 0,946 0,947 0,980 p-valor del ks.test 7,07e−6 0,126 7,8e−4 No aplica No aplica No aplica 8,88e−16 Cuadro 4.11: Tabla con el error empírico y el resultado del test de KS para medir la adecuación de los estadístico a su comportamiento asintótico Los resultados son quizás los esperables. Los estadísticos de divergencia parecen tender efectivamente a su comportamiento asintótico (el estadístico G2 no ha estado denido para más de la mitad de las muestras, con lo que su valor está muy sesgado en perjuicio de aquellas que tienen intervalos vacíos), aunque se encuentren aún lejos de tal distribución. Probando a aumentar el valor de n parece que efectivamente su comportamiento se acerca cada vez más al de una χ2 7 , si bien el test que hemos diseñado tiene muchísima variabilidad, seguramente por causa del método aleatorio elegido para la elección de las fronteras. A falta de criterios más claros y de una mayor fuerza computacional, estos test, tal y como los hemos implementado, no parecen idóneos para el contraste de la hipótesis compuesta. Los test de la segunda familia son los que tienen un comportamiento más regular. Obviamente, como nosotros mismos hemos estimado los cuantiles de CvM y AD, era esperable que el error empírico fuese muy cercano a 0,95 . También se ha comprobado que la fórmula 4 El estadístico no está denido para 696 muestras 52 CAPÍTULO 4. ILUSTRACIÓN EN BASES DE DATOS SIMULADOS de Lilliefors para hallar el p valor da una buena aproximación. Siguiendo estos resultados, parece natural que efectivamente estos test sean los más populares para medir la normalidad de entre todos los que hemos propuesto. Su desventaja es que están limitados a dicha distribución, o, adaptados, a otras distribuciones invariantes a cambios de escala y locación. El test de Fan, como ya se podría anticipar a juzgar por la sección anterior, no tiene el comportamiento esperado. A pesar de que su error empírico no está demasiado lejos del esperable, su esperanza empírica es claramente negativa para todas las veces que hemos corrido el test, y su distribución está muy lejos de ser una normal estándar, dejando pocas esperanzas a que se normalice con una aumento de la n . La causa de este comportamiento errático probablemente se encuentre en el criterio de la elección de la h , que nosotros hemos establecido de forma bastante arbitraria (aunque cumpliendo con las restricciones asintóticas). Sería interesante entonces probar a establecer métodos de selección que tomen los anchos de banda que mejor se ajusten a la muestra, y estudiar su efecto en el comportamiento asintótico. Hemos probado con la regla del pulgar y el resultado no ha sido mejor, de hecho no suele ser recomendado en el contexto de los estimadores de densidad. El método que sí se suele sugerir es el de validación cruzada, pero precisa resolver un problema de optimización e implementarlo se haría muy pesado. Por tanto no hemos ido más allá y hemos decidido terminar aquí el análisis de nuestros test. Conclusión Empezábamos la introducción de este trabajo defendiendo la elección de las tres familias de test que hemos incluido sobre todos los tipos existentes de contrastes de especicación. Alegábamos que lo que las destaca por encima de las demás es su importancia histórica y conceptual, dos criterios que consideramos fundamentales para un trabajo que quiere funcionar como una introducción teórica a las técnicas más fundamentales con las que medir la divergencia entre una muestra y una distribución teórica. Hemos visto, sin embargo, que su implementación en la práctica está lejos de ser trivial a pesar de lo intuitivo de su losofía. Por supuesto, siempre se pueden emplear paquetes estadísticos que incluyan test ya diseñados que realicen el contraste en R o cualquier otro lenguaje de programación. Pero si intentamos adentrarnos en el código sobre el que se basan estos test, vemos que algunos de ellos, y en especial para la hipótesis simples, tienen que sortear problemas como la elección de los intervalos para el test de la X2 de Pearson o el bandwidth para los estimadores de densidad kernel. Estas cuestiones que, en general, no tienen una respuesta convincente que surja de la teoría de la probabilidad, y que son abordadas utilizando técnicas computacionales. Intentar integrar en un solo papel la base teórica de cada una de las familias, junto con una descripción de los métodos computacionales que se usan para estudiarlas, e incluir los resultados que se han obtenido hasta la fecha sería una empresa tremenda. Este trabajo solo intenta ser una primera aproximación a ese ámbito mucho más grande que son los contrastes de especicación A pesar de ello, hemos obtenido resultados bastante satisfactorios. Aún con lo rudimentario del código que hemos usado y de las limitaciones computacionales, hemos visto que para la hipótesis simple nuestros test tienen un buen comportamiento. Incluso, hemos obtenido resultados correctos cuando hemos aplicado los test del segundo capítulo al contraste de una hipótesis compuesta, y alentadores cuando usamos los estadísticos de divergencia. En todo caso, el método seguido dista de poder ser considerado riguroso, y siempre nos hemos limitado a la normal estándar y a una hipótesis nula correcta. Quedaría por ver, por ejemplo, qué ocurre cuando realizamos el contraste para una distribución discreta, o cuando nuestra hipótesis nula es diferente a la distribución que usamos para 53 54 CAPÍTULO 4. ILUSTRACIÓN EN BASES DE DATOS SIMULADOS simular nuestras muestras. Ya hemos discutido los motivos por los cuales probablemente los estadísticos de Bickel y Rosenblatt no hayan rendido como era esperable. Lo cierto es que de entre las tres familias de contrastes, la basada en la densidad es la más reciente y menos estudiada. Además, la relativa complejidad de sus estadísticos hace que sea con diferencia la menos popular de las tres, a pesar del enorme interés que parece suscitar la estimación de la densidad en ámbitos como la econometría o las nanzas. Sin embargo, la popularidad no ha sido un factor tan determinante a la hora de elegir los estadísticos. Si no, deberíamos haber hablado, por ejemplo, del estadístico de Shapiro-Wilk, que se encuentra entre los más populares para contrastar la normalidad de un conjunto de datos, y que se basa en la comparación de dos estimadores alternativos de la varianza de la muestra. Además, nos hemos dejado otras grandes familias de contrastes de especicación, que por lo menos querríamos introducir a continuación. Cuando en el capítulo segundo introducíamos la fórmula general de un estadístico basado en la FDE usábamos la expresión: Tn=c(n)d(ˆ Fn, F0), con d(., .) una función de divergencia (aunque al nal siempre hemos utilizado distancias). Esta expresión es muy general, y nosotros siempre hemos hecho depender al estadístico de la FDE y la F0 a través del proceso empírico. De esta forma teníamos que B(x)=0∀x⇔ F(x) = F0(x)∀x , donde usábamos B para denotar el proceso empírico. Pero existen más maneras de caracterizar una distribución que a través de su función de distribución (o de su densidad, si la admite). También la función cuantil dada por F−1(x) , o la función característica φF(x) = Ef[exp(ixY )] sirven para caracterizar una distribución. De esta forma, podríamos considerar otras opciones para B como: B(x) = F−1(x)−F−1 0(x) B(x) = φF(x)−φF0(x) Ambas posibilidades siguen cumpliendo la condición B(x) = 0 ∀x⇔F(x) = F0(x)∀x , con lo que nos exponen dos otras grandes familias de contrastes de especicación: las basadas en el estimador de la función cuantil y el estimador de la función característica respectivamente. Abarcarlas en este trabajo era casi imposible por limitaciones de tiempo y espacio, con lo que nos conformamos con este pequeño esboce. Aquí terminamos nuestro trabajo. A continuación presentamos el código que hemos usado en el capítulo cuarto. A pesar de ser rudimentario, y sin duda mejorable, los test de CONCLUSIÓN 55 los paquetes de R han conrmado que funciona correctamente. Lo hemos presentado tal y como ha quedado luego de cada análisis, comenzando con la denición de los parámetros del test de Monte Carlo, y terminando por sus resultados. En medio del cuerpo se genera la muestra y se calculan los estadísticos. Hemos incluido comentarios para que sea más legible. Para más información sobre las diversas cuestiones que hemos tratado, o para profundizar más en los contrastes de especicación, recomendamos leer la bibliografía que presentamos al nal del trabajo. En ella, destaca especialmente el manual Comparing Distributions de Oliver Thas, que ha servido como la principal referencia para este trabajo. 56 CONCLUSIÓN Código de R .1. Hipótesis simple .1.1. Primer capítulo #DEFINIMOS EL TEST QUE NOS PERMITIRA REPETIR VARIAS VECES EL METODO \begin{ l s t l i s t i n g }[ language= R ] n=1000 rep =1000 tp=1: rep tr =1: rep tc =1: rep vp=1: rep vr=1: rep vc=1: rep for ( j in 1: rep ){ #GENERAMOS DATOS DE UNA NORMAL ESTANDAR mu=0 sig=1 data = rnorm (n ,mu, sig ) #Generamos intervalos equiprobables con prob 1 / k k=10 x= seq (0 ,1 ,1 / k) fron= qnorm (x ,mu, sig ) p=1 / k vp= rep (p,10) #Creamos la multinomial 57 64 CÓDIGO DE R #LA FUNCION: con VECTORIZADA conv < − function (s) sapply (s , con) #YA PODEMOS HALLAR FACILMENTE EL VALOR DEL ESTADISTICO int2 < − function (x) (kde(x) − conv(x)) ∗∗ 2 T2=integrate ( int2 , − Inf , Inf ) $ value #NOS INTERESA AHORA PODER REALIZAR EL TEST. TENEMOS QUE HALLAR LOS PARAMETROS. #CALCULAMOS LA MEDIA Y LA VARIANZA ASINTOTICA k2 < − function (x) dnorm (x ,0 ,1) ∗∗ 2 ik2=integrate (k2, − Inf , Inf ) $ value mu=ik2 / (n ∗ h) #PARA HALLAR LA VARIANZA TENGO QUE EVALUAR UNA INTEGRAL DOBLE. #DEFINO LA FUNCION INTEGRANDO kk < − function (u , v) dnorm (v+u,0 ,1) ∗ dnorm (v ,0 ,1) #DEFINO LA FUNCION DE LA PRIMERA INTEGRAL (COMO ikk ES UNA INTEGRAL #NO PUEDO APLICARLE integrate () , NECESITO VECTORIZARLA ANTES. DEFINO ikkv ) ikk < − function ( s ) integrate ( function (v , u=s ) dnorm (v+u,0 ,1) ∗ dnorm (v ,0 ,1) , − Inf , Inf ) $ value ikkv < − function (s) sapply (s , ikk ) #LA ELEVO AL CUADRADO ikkv2 < − function (x) ikkv (x) ∗∗ 2 #YA PUEDO INTEGRARLA NORMALMENTE Y HALLAR LA VARIANZA DEL ESTADISTICO sig= sqrt (2 ∗ integrate ( ikkv2 , − Inf , Inf ) $ value ∗ ik2 ) #PODEMOS YA REALIZAR LOS TEST EMPLEANDO LA DISTRIBUCION ASINTOTICA: #TEST1: 1 − pnorm (n ∗ sqrt (h) ∗ (T1 − mu) / sig ,0 ,1) #TEST2: 1 − pnorm (n ∗ sqrt (h) ∗ (T2 − mu) / sig ,0 ,1) .2. HIPÓTESIS COMPUESTA 65 t1 [ j]=1 − pnorm (n ∗ sqrt (h) ∗ (T1 − mu) / sig ,0 ,1) >0.05 t2 [ j]=1 − pnorm (n ∗ sqrt (h) ∗ (T2 − mu) / sig ,0 ,1) >0.05 v1 [ j ]=n ∗ sqrt (h) ∗ (T1 − mu) / sig v2 [ j ]=n ∗ sqrt (h) ∗ (T2 − mu) / sig } #RESULTADOS DE LOS ANaLISIS DE LOS ESTADISTICOS: #ERRORES EMPIRICOS sum ( t1 ) /rep sum ( t2 ) /rep #ks . test ks . test (v1 , pnorm ,0 ,1) ks . test (v2 , pnorm ,0 ,1) .2. Hipótesis compuesta .2.1. Primer capítulo library ( stats4 ) rep =1000 tp=1: rep tr =1: rep tc =1: rep vp=1: rep vr=1: rep vc=1: rep for ( j in 1: rep ){ #GENERAMOS LA MUESTRA ALEATORIA n=100 data = rnorm (n,0 ,1) q =2 #CALCULAMOS EL ESTIMADOR MAXIMO VEROSIMIL DE LA MEDIA Y LA VARIANZA CON EL COMANDO mle #DEFINIMOS LA FUNCION LOG − VER: lv < − function (mu, sig ){ r= dnorm ( data ,mu, sig ) 66 CÓDIGO DE R # − sum ( log ( r )) } #USAMOS mle PARA HALLAR LOS PARaMETROS QUE LA MINIMIZAN mle( minuslogl=lv , start = list (mu=1, sig =1)) #OBTENEMOS NaNs DEBIDO A QUE EN EL PROCESO DE OPTIMIZACION LA FUNCION INTENTA #ASIGNAR VALORES NEGATIVOS A sig . NO INFLUYE EN EL RESULTADO, LOS ESTIMADORES SON #CORRECTOS. #ACCEDEMOS A ELLOS DE LA MANERA SIGUIENTE mu=mle( minuslogl=lv , start = list (mu=1, sig =1))@coef [ 1] sig=mle( minuslogl=lv , start = list (mu=1, sig =1))@coef [ 2] f0 < − function (x) dnorm (x ,mu, sig ) F0 < − function (x) pnorm (x ,mu, sig ) # −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− #CAPiTULO 1 #GENERAMOS LA MULTINOMIAL k=10 fron=1:k+1 fron[1]= − Inf fron [ k+1]=Inf fron [ 2: k]= sort ( runif (k − 1, min ( data ) , max ( data ))) mult=1:k for ( i in 1:k){ mult [ i ]= sum ( fron [ i+1]> data&data >fron [ i ]) } #CALCULAMOS LAS PROBABILIDADES ASOCIADAS A CADA INTERVALO p=1:k for ( i in 1:k){ p[ i ]=integrate ( f0 , fron [ i ] , fron [ i +1]) $ value .2. HIPÓTESIS COMPUESTA 67 } #CALCULAMOS EL ESTADISTICO DE PEARSON x2=0 for ( i in 1:k){ x2=x2+((mult [ i ] − n ∗ p[ i ]) ∗∗ 2) / (n ∗ p[ i ]) } #CALCULAMOS EL ESTADISTICO DE MAX.VER. g2=0 for ( i in 1:k){ g2=g2+2 ∗ mult [ i ] ∗ ( log (mult [ i ] / n) − log (p[ i ]))} #CALCULAMOS EL ESTADISTICO DE CRESSIE AND READ CON LAMBDA=2 / 3 lamb=2 / 3 pd=0 for ( i in 1:k){ pd=pd+mult [ i ] ∗ (( mult [ i ] / (n ∗ p[ i ]) ) ∗∗ lamb − 1)} pd=2 / (lamb ∗ (lamb+1)) ∗ pd #HE COMPROBADO TOMANDO LAMBDA=1 QUE DA EL MISMO RESULTADO QUE X2 vp [ j ]=x2 vr [ j ]=g2 vc [ j ]=pd tp [ j]=1 − pchisq (x2 ,k − 1 − q )>0.05 tr [ j]=1 − pchisq (g2 ,k − 1 − q )>0.05 tc [ j]=1 − pchisq (pd ,k − 1 − q )>0.05 } #RESULTADOS DEL TEST: #TP: sum (tp ) /rep #TR: 68 CÓDIGO DE R sum ( tr , na . rm =TRUE) / ( rep − sum ( is . na ( tr ))) sum ( is . na ( tr )) #TC: sum ( tc ) /rep #TESTEO QUE EFECTIVAMENTE LOS ESTADISTICOS SIGUEN UN CHISQ ks . test (vp , pchisq ,k − 1 − q ) ks . test (vr , pchisq ,k − 1 − q ) ks . test (vc , pchisq ,k − 1 − q ) .2.2. Capítulo 2 library (" goftest ") #COMO NO CONTAMOS CON TABLAS CON LA DISTRIBUCION DE AD Y CvM PARA TESTEAR NORMALIDAD #VAMOS A HALLAR EL CUANTIL DE ORDEN 95 POR MONTECARLO rep2=10000 n=100 vw=1:rep2 va=1: rep for ( j in 1: rep2 ){ #GENERAMOS LOS DATOS. NOS VALE CUALQUIER TIPO DE NORMAL data = rnorm (n,3 ,9) #ESTANDARIZAMOS Y ORDENAMOS data =( data − mean ( data )) /sd ( data ) data = sort ( data ) #HALLAMOS EL VALOR DE LOS ESTADISTICOS sum =0 for ( i in 1:n){ sum = sum +(2 ∗ i − 1) ∗ ( log ( pnorm ( data [ i ] ,0 ,1))+ log (1 − pnorm ( data [n − i +1] ,0 ,1))) } A= − n − (1 / n) ∗ sum va [ j ]=A .2. HIPÓTESIS COMPUESTA 69 W=0 for ( i in 1:n){ W=W+( pnorm ( data [ i ] ,0 ,1) − (2 ∗ i − 1) / (2 ∗ n)) ∗∗ 2 } vw[ j ]=W } #HALLAMOS EL CUANTIL EMPIRICO qam= sort (va )[ rep2 ∗ 0.95] qwm= sort (vw)[ rep2 ∗ 0.95] qam qwm # −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− rep =1000 #VECTORES DONDE GUARDAREMOS LOS RESULTADOS DE LOS TEST tk=1: rep tw=1: rep ta=1: rep for ( j in 1: rep ){ #GENERAMOS DATOS DE UNA NORMAL ESTANDAR Y LOS ORDENAMOS data = rnorm (n,0 ,1) data = sort (( data − mean ( data )) /sd ( data )) #CONSTRUIMOS LA FUNCION DE DISTRIBUCION EMPIRICA Y LA VISUALIZAMOS CON LA REAL fde=ecdf ( data ) x= seq ( − 4,4,0.01) y= pnorm (x ,0 ,1) plot ( fde ) lines (x ,y , " l " , lty =3, col ="red") #CONSTRUIMOS EL ESTADISTICO DE KS (SIN MULTIPLICARLO POR SQRT(N) DE MOMENTO) D =0 for ( i in 1:n){ dif1= abs ( i / n − pnorm ( data [ i ] ,0 ,1)) 70 CÓDIGO DE R dif2= abs (( i − 1) / n − pnorm ( data [ i ] ,0 ,1)) dif= max ( dif1 , dif2 ) if ( dif> D ) D =dif } #PROBAMOS AHORA CON LOS ESTADISTICOS DE ANDERSON − DARLING. sum =0 for ( i in 1:n){ sum = sum +(2 ∗ i − 1) ∗ ( log ( pnorm ( data [ i ] ,0 ,1))+ log (1 − pnorm ( data [n − i +1] ,0 ,1))) } A= − n − (1 / n) ∗ sum W=0 for ( i in 1:n){ W=W+( pnorm ( data [ i ] ,0 ,1) − (2 ∗ i − 1) / (2 ∗ n)) ∗∗ 2 } #GUARDAMOS LOS RESULTADOS tk [ j ]= D <(0.886 /sqrt (n)) tw [ j ]=W<qwm ta [ j ]=A<qam } sum ( tk ) /rep sum (tw) /rep sum ( ta ) /rep #LLAMAMOS A LOS PAQUETES QUE USAREMOS library ( kdensity ) library ( stats4 ) #DEFINIMOS LOS PARAMETROS DEL TEST rep =100 n=500 v2=1: rep t2=1: rep for ( j in 1: rep ){ #GENERAMOS LA MUESTRA ALEATORIA .2. HIPÓTESIS COMPUESTA 71 data = rnorm (n,0 ,1) q =2 #CALCULAMOS EL ESTIMADOR MAXIMO VEROSIMIL DE LA MEDIA Y LA VARIANZA CON EL COMANDO mle #DEFINIMOS LA FUNCION LOG − VER: lv < − function (mu, sig ){ r= dnorm ( data ,mu, sig ) # − sum ( log ( r )) } #USAMOS mle PARA HALLAR LOS PARAMETROS QUE LA MINIMIZAN mle( minuslogl=lv , start = list (mu=1, sig =1)) #OBTENEMOS NaNs DEBIDO A QUE EN EL PROCESO DE OPTIMIZACION LA FUNCION INTENTA #ASIGNAR VALORES NEGATIVOS A sig . NO INFLUYE EN EL RESULTADO, LOS ESTIMADORES SON #CORRECTOS. #ACCEDEMOS A ELLOS DE LA MANERA SIGUIENTE mu=mle( minuslogl=lv , start = list (mu=1, sig =1))@coef [ 1] sig=mle( minuslogl=lv , start = list (mu=1, sig =1))@coef [ 2] f0 < − function (x) dnorm (x ,mu, sig ) F0 < − function (x) pnorm (x ,mu, sig ) # −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− #CONSTRUIMOS EL ESTIMADOR DE DENSIDAD KERNEL h=1.06 ∗ sd ( data ) ∗ n ∗∗ ( − 1 / 5) kde=kdensity ( data , kernel="gaussian" ,bw=h) #SEGUNDO ESTADISTICO DE BICKEL Y ROSENBLATT #TENEMOS QUE DEFINIR LA FUNCION RESULTANTE DEL PRODUCTO DE CONVOLUCION con < − function ( s ) integrate ( function (y , x=s ) ( dnorm ((x − y) / h,0 ,1) / h) ∗ dnorm (y ,mu, sig ) , − Inf , Inf ) $ value #COMO con SE DEFINE CON EL COMANDO INTEGER, QUE NO ES UNA FUCCION VECTORIZADA, #NO PODRIAMOS APLICAR DE NUEVO INTEGER PARA HALLAR EL VALOR DEL ESTADISTICO. 72 CÓDIGO DE R #SIN EMBARGO ESTE INCONVENIENTE PUEDE ARREGLARSE USANDO sapply Y DEFINIMOS ASI #LA FUNCION: con VECTORIZADA conv < − function (s) sapply (s , con) #YA PODEMOS HALLAR FACILMENTE EL VALOR DEL ESTADISTICO int2 < − function (x) (kde(x) − conv(x)) ∗∗ 2 T2=integrate ( int2 , − Inf , Inf ) $ value #NOS INTERESA AHORA PODER REALIZAR EL TEST. TENEMOS QUE HALLAR LOS PARAMETROS. #CALCULAMOS LA MEDIA Y LA VARIANZA ASINTOTICA k2 < − function (x) dnorm (x ,0 ,1) ∗∗ 2 ik2=integrate (k2, − Inf , Inf ) $ value mut=ik2 / (n ∗ h) #PARA HALLAR LA VARIANZA TENGO QUE EVALUAR UNA INTEGRAL DOBLE. #DEFINO LA FUNCION INTEGRANDO kk < − function (u , v) dnorm (v+u,0 ,1) ∗ dnorm (v ,0 ,1) #DEFINO LA FUNCION DE LA PRIMERA INTEGRAL (COMO ikk ES UNA INTEGRAL #NO PUEDO APLICARLE integrate () , NECESITO VECTORIZARLA ANTES. DEFINO ikkv ) ikk < − function ( s ) integrate ( function (v , u=s ) dnorm (v+u,0 ,1) ∗ dnorm (v ,0 ,1) , − Inf , Inf ) $ value ikkv < − function (s) sapply (s , ikk ) #LA ELEVO AL CUADRADO ikkv2 < − function (x) ikkv (x) ∗∗ 2 #HALLO LA INTEGRAL DEL CUADRADO DE LA FUNCION ESTIMADOR ikd2=integrate ( function (x) kde(x) ∗∗ 2, − Inf , Inf ) $ value #YA PUEDO INTEGRARLA NORMALMENTE Y HALLAR LA VARIANZA DEL ESTADISTICO sigt= sqrt (2 ∗ integrate ( ikkv2 , − Inf , Inf ) $ value ∗ ikd2 ) #PODEMOS YA REALIZAR EL TEST EMPLEANDO LA DISTRIBUCION ASINTOTICA: #TEST2: 1 − pnorm (n ∗ sqrt (h) ∗ (T2 − mut) / sigt ,0 ,1) .2. HIPÓTESIS COMPUESTA 73 t2 [ j]=1 − pnorm (n ∗ sqrt (h) ∗ (T2 − mut) / sigt ,0 ,1) >0.05 v2 [ j ]=n ∗ sqrt (h) ∗ (T2 − mut) / sigt } #RESULTADOS DEL TEST sum ( t2 ) /rep ks . test (v2 , pnorm ,0 ,1)