La estimación no paramétrica de la densidad

Traballo Fin de Grao La estimación no paramétrica de la densidad Laura Cotos Fernández-Arruty Julio,2022 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao La estimación no paramétrica de la densidad Laura Cotos Fernández-Arruty Julio, 2022 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Trabajo propuesto Área de Conocimiento: Estadística e Investigación Operativa Título: La estimación no paramétrica de la densidad Breve descripción del contenido La densidad es una función característica de muchísimo interés en la Estadística que fue y sigue siendo ampliamente estudiada desde diferentes perspectivas. Históricamente se propusieron diferentes modelos paramétricos cuyo buen funcionamiento recae sobre la vericación de determinadas condiciones que o bien no siempre se cumplen o bien no pueden ser testadas en la práctica. Por esta razón surgieron los métodos no paramétricos: tipo núcleo, splines, wavelets... que son menos restrictivos en las condiciones requeridas y ofertan una mayor exibilidad en su aplicación. En este trabajo haremos una revisión sobre la estimación de la densidad en general, centrándonos en los procedimientos no paramétricos y más particularmente en los estimadores tipo núcleo. A título orientativo el trabajo se puede organizar en las siguientes secciones: Revisión de la función de densidad y su estimación. Métodos no paramétricos para la densidad: la estimación tipo núcleo. El problema de selección de la ventana. Estudio de simulación y/o análisis de datos reales. Recomendaciones Otras observaciones iii Índice Resumen vii 1 Introducción 1 1.1 Conceptosesenciales .................................. 1 1.2 El problema de estimación de la densidad . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Casodeestudio .................................... 5 2 El estimador tipo núcleo de la densidad 11 2.1 Criteriosdeerror.................................... 13 3 Métodos de selección del parámetro ventana 21 3.1 Regla del pulgar de Silverman . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.2 Reglaplug-in ...................................... 22 3.3 Validacióncruzada................................... 24 4 Estudio de simulación 27 5 Aplicación a datos reales 37 6 Conclusiones y trabajo futuro 39 I Código para el estudio de simulación 41 II Código para la aplicación a datos reales 53 Referencias 57 v Resumen La estimación de la densidad es uno de los campos de mayor interés dentro de la Estadística, tanto por su aspecto propio como por las múltiples aplicaciones que de él derivan. En este trabajo nos centramos en el estimador tipo núcleo de la densidad, que es uno de los métodos no paramétricos más empleados en la literatura. Comenzaremos por introducir aquellos conceptos esenciales para su construcción, para luego detallar y derivar sus principales características. A pesar de la exibilidad que coneren los métodos no paramétricos, la cuestión más esencial para el cálculo del estimador tipo núcleo de la densidad es la elección del denominado parámetro ventana. La elección de este parámetro ventana no es inmediata, y a lo largo de los años se han desarrollado diferentes métodos basados en ideas diversas. En este trabajo nos centraremos en los tres más empleados: el conocido como método de Silverman, el método de Sheather y Jones y el método de validación cruzada. Uno de los objetivos es comparar el comportamiento de estos selectores tanto entre sí como con respecto al óptimo teórico. Para ello realizamos un completo estudio de simulación empleando el paquete estadístico . En dicho estudio se analiza el comportamiento de los procedimientos en muestras nitas con varios tamaños muestrales y cuatro modelos teóricos que cubren una amplia gama de las posibles características presentes en una función de densidad. Por último, ilustraremos el uso del estimador tipo núcleo y los selectores estudiados con un conjunto de datos reales del ámbito biomédico en el que se recogen diversas medidas útiles en la detección de la cardiomegalia (condición humana que se caracteriza por un tamaño del corazón anormalmente grande). vii 6 1. Introducción Figura 1.2: ilustración de una radiografía de tórax obtenida de (Sogancioglu y cols., 2020) con máximo diámetro horizontal torácico 266.72 mm y máximo diámetro horizontal cardíaco 146.86 mm. en Estados Unidos, están disponibles en la base de datos SCR (van Ginneken, Stegmann, y Loog, 2006a) y se describen con detalle en (Shiraishi y cols., 2000) y (van Ginneken, Stegmann, y Loog, 2006b). mdhc mdht ctr Mínimo 76.3 223.3 0.273 Q1 101.3 265.6 0.357 Mediana 110.2 283.1 0.388 Media 111.1 285.1 0.392 Desviación típica 13.251 24.997 0.052 Q3 120.4 302.9 0.427 Máximo 151.2 349.3 0.538 Tabla 1.1: resumen numérico de la base de datos del estudio. En la Tabla 1.1 hemos realizado un resumen numérico de las tres variables anteriormente mencionadas, mdhc , mdht y ctr . Observamos que la variable mdhc tiene valores más pequeños que mdht , mientras que la variable ctr tiene valores entre 0 y 1 , pues es el cociente entre las otras dos. En la Figura 1.3 se muestran los diagramas de caja para cada una de las variables de interés. En los diagramas de caja podemos observar, igual que en la tabla anterior, que el rango de valores de la variable mdhc es menor que el de la variable mdht . También se aprecia la distinción entre las escalas de las dos primeras variables con la variable ctr . 1.3. Caso de estudio 7 Figura 1.3: diagramas de caja correspondientes a las variables de interés, mdhc , mdht y ctr de izquierda a derecha. Aunque los diagramas de caja nos aportan cierta información sobre las muestras, parece necesario evaluar los histogramas de las mismas, que pueden verse en la Figura 1.4. En los histogramas de las tres muestras observamos que se tratan de distribuciones unimodales y con bastante simetría en torno al punto medio, si bien es cierto que el histograma correspondiente a mdht pueda presentar cierta asimetría hacia la izquierda aunque no parece muy marcada. Para ilustrar con este ejemplo práctico las limitaciones del histograma, podemos ver en la Figura 1.5 dos nuevos histogramas para las mismas muestras, y resulta menor que en los mostrados en la Figura 1.4. En este caso se ve un histograma mucho menos rugoso, pero que parece mantener las propiedades de simetría y unimodalidad comentados anteriormente. Este trabajo n de grado se organiza en 6 capítulos, siendo el primero esta Introducción que acabamos de hacer. En el Capítulo 2 denimos el estimador tipo núcleo de la densidad de manera formal junto con varios criterios de error que utilizaremos en otros capítulos. En el Capítulo 3 mostramos tres de los procedimientos de selección del parámetro ventana más conocidos para el estimador tipo núcleo. 8 1. Introducción Figura 1.4: histogramas correspondientes a las variables mdhc , mdht y ctr de izquierda a derecha con 30 intervalos. El Capítulo 4 es un estudio de simulación en el que comparamos los selectores introducidos en el Capítulo 3 entre sí y estudiamos su comportamiento. En el Capítulo 5 estimamos la densidad de un conjunto de datos reales y aplicamos los selectores estudiados. Por último, el Capítulo 6 consta de las conclusiones del presente trabajo y de una breve exposición de trabajos futuros que se podrían realizar a partir de este. 1.3. Caso de estudio 9 Figura 1.5: histogramas correspondientes a las variables mdhc , mdht y ctr de izquierda a derecha con diez intervalos. 10 1. Introducción Capítulo 2 El estimador tipo núcleo de la densidad En este capítulo vamos a introducir en el estimador tipo núcleo de la densidad y a estudiar tanto a nivel teórico como práctico sus principales características. Supongamos que disponemos de una muestra aleatoria simple, X1, ..., Xn , de una variable aleatoria X con función de densidad f . El estimador tipo núcleo de la función de la densidad se dene según (Parzen, 1962) y (Rosenblatt, 1956) como: ˆ fnh(x) = 1 nh n X i=1 Kx−Xi h, (2.1) donde h es un número real positivo que se denomina ancho de banda, parámetro ventana o parámetro de suavización, y K es lo que se conoce como función núcleo, una función unimodal simétrica respecto al cero. El estimador tipo núcleo puede verse como una suma de protuberancias de las observaciones, donde la función tipo núcleo K determina la forma de las protuberancias y h determina el ancho de las mismas. Esta interpretación se ilustra en la Figura 2.1. Como acabamos de ver, el parámetro ventana, h , representa la amplitud de la función núcleo centrado en las observaciones {xi}n i=0 , esto es algo así como en qué rango de valores de la muestra nos jamos para hacer la estimación en un punto dado, aquellas observaciones de la muestra más cercanas al punto de interés aportarán más peso a la estimación que las más lejanas. Como veremos con más detalle en el Capítulo 3, su elección es compleja y tiene mucha inuencia sobre el resultado de la estimación. Si la ventana toma un valor muy pequeño, tendremos una estimación con mucho ruido, lo que signica que se centra demasiado en cada dato particular y la estimación es rugosa (apareciendo lo que se conoce como el fenómeno de infrasuavización). Por el contrario, si la ventana toma un valor grande, lo que ocurre es que suaviza la estimación 11 12 2. El estimador tipo núcleo de la densidad Figura 2.1: ilustración de la construcción del estimador tipo núcleo de la densidad con una muestra de tamaño n= 10 escogida aleatoriamente de la varibale mdht de los datos de estudio. en exceso, perdiendo información necesaria y generando un fenómeno de sobresuavización. En la Figura 2.2 hemos representado la densidad de la variable mdht de los datos presentados en el Capítulo 1, con dos elecciones de parámetro ventana distintas ( h= 1 en la izquierda y h= 50 en la derecha), así podemos observar ese gran efecto del parámetro ventana que comentamos anteriormente. En la Figura 2.2 (izquierda), observamos una estimación rugosa con gran cantidad de picos, mientras que en la Figura 2.2 (derecha), observamos que la curva está achatada. Otro de los elementos presentes en el estimador tipo núcleo denido en (2.1) es la función núcleo, K , que recordemos es una función real, no negativa e integrable con RRK(x)dx= 1 , unimodal y simétrica respecto al origen. Como ejemplos de esta función núcleo podríamos citar los más habituales: el núcleo Normal que se corresponde con la función de densidad Normal estándar, o el núcleo Epanechnikov, denido por: K(x) = 3 4(1 −x2)1{|x|<1}, donde 1{·} denota la función indicadora de conjuntos. Se pueden consultar más funciones tipo 2.1. Criterios de error 13 Figura 2.2: estimaciones tipo núcleo de la variable mdht con distintos anchos de banda, h= 1 (izquierda) y h= 50 (derecha). núcleo en una tabla del apéndice B de (Wand y Jones, 1994). Es sabido que variaciones en la función núcleo no tienen una inuencia tan importante en el resultado de la estimación como el parámetro ventana. Ejemplo de ello es la Figura 2.3, donde podemos observar dos estimaciones tipo núcleo de la densidad de la muestra de estudio con el núcleo Normal (izquierda) y con el de Epanechnikov (derecha), entre las que apenas se aprecian variaciones. 2.1. Criterios de error Para poder valorar de manera justa cualquier procedimiento de estimación, es necesario denir mecanismos que midan la bondad de ajuste de los mismos, son los denominados criterios de error, que miden la discrepancia entre la estimación y el valor teórico o real. Una posibilidad habitual es emplear como medida de discrepancia el Error Cuadrático Medio ( ECM ). Dado un parámetro θ y una estimación del mismo, ˆ θ , se dene en (Wand y Jones, 1994) el ECM como: ECM(ˆ θ) = Eˆ θ−θ2. El ECM posee una expresión más fácilmente interpretable que depende del concepto de sesgo y varianza. 14 2. El estimador tipo núcleo de la densidad Figura 2.3: estimaciones tipo núcleo de la densidad para la variable mdht con distintas funciones núcleo: núcleo Normal (izquierda) y núcleo Epanechnikov (derecha). Denición 2.1. Denominamos sesgo de un estimador, ˆ θ , para un parámetro poblacional θ , a Sesgo(ˆ θ) = E[ˆ θ]−θ, (2.2) siempre y cuando exista E[ˆ θ] . Diremos que el estimador es insesgado si su sesgo es nulo. Nótese que esta cualidad de insesgadez es una propiedad deseable para cualquier estimador. Vamos ahora entonces a reescribir la expresión del ECM en función del sesgo y varianza para un estimador cualquiera ˆ θ de un valor teórico θ : ECM(ˆ θ) = Eˆ θ−θ2=Eh(ˆ θ−E[ˆ θ] + E[ˆ θ]−θ)2i =Eh(ˆ θ−E[ˆ θ])2+ 2(ˆ θ−E[ˆ θ])(E[ˆ θ]−θ)+(E[ˆ θ]−θ)2i = (a) Eh(ˆ θ−E[ˆ θ])2i+ 2Eh(ˆ θ−E[ˆ θ])(E[ˆ θ]−θ)i+Eh(E[ˆ θ]−θ)2i = (b) Var hˆ θi+ 2Ehˆ θ−E[ˆ θ]iEhE[ˆ θ]−θi+E[ˆ θ]−θ2= (c) Var hˆ θi+Sesgo2(ˆ θ), (2.3) donde en (a) utilizamos la linealidad de la esperanza, en (b) utilizamos la denición de varianza y en (c) la denición de sesgo de (3.1), teniendo en cuenta que, de nuevo, por la linealidad de la esperanza, Ehˆ θ−E[ˆ θ]i= 0 . Nuestro objetivo ahora es obtener la expresión del ECM del estimador tipo núcleo de la densidad para poder estudiar su comportamiento. Haciendo uso de la descomposición anterior, necesitaremos obtener su media y su varianza. 2.1. Criterios de error 15 Utilizando la siguiente notación, Kh(·) = 1 hK(· h) , podemos reescribir (2.1) de manera más compacta: ˆ fnh(x) = 1 n n X i=1 Kh(x−Xi). (2.4) En primer lugar procederemos al cálculo de la media para un punto x∈Sop(f)⊂R cualquiera, donde Sop(f) denota el soporte de la función de densidad. Ehˆ fnh(x)i=E"1 n n X i=1 Kh(x−Xi)#= (a) 1 nnE[Kh(x−Xi)] = (b)ZR Kh(x−y)f(y)dy, (2.5) donde en (a) hemos utilizado la linealidad de esperanza y en (b) aplicamos la denición de esperanza de una variable aleatoria continua, en particular de una transformación de una variable. Podemos escribir entonces el sesgo del estimador en un punto como: Sesgo hˆ fnh(x)i=Ehˆ fnh(x)i−f(x) = ZR Kh(x−y)f(y)dy−f(x)=(Kh⊛f)(x)−f(x), (2.6) donde ⊛ denota la convolución entre funciones, esto es, dadas g1 y g2 dos funciones reales de variable real, se dene (g1⊛g2)(x) = RRg1(x−y)g2(y)dy. Una vez obtenida la expresión explícita del sesgo, calcularemos su varianza: Var hˆ fnh(x)i=Var "1 n n X i=1 Kh(x−Xi)#= (a) 1 n2nVar [Kh(x−X)] = 1 nVar [Kh(x−X)] = (b) 1 nEK2 h(x−X)−E[Kh(x−X)]2 = (c) 1 nZR K2 h(x−y)f(y)dy−1 nZR Kh(x−y)f(y)dy2 =1 nh(K2 h⊛f)(x)−(Kh⊛f)2(x)i, (2.7) donde en (a) aplicamos propiedades de la varianza, en (b) usamos la denición de varianza y en (c) la denición de la esperanza. Combinando el sesgo al cuadrado con la varianza tal y como se indica en (2.3) obtenemos el error cuadrático medio del estimador en un punto dado: ECM hˆ fnh(x)i=1 nh(K2 h⊛f)(x)−(Kh⊛f)2(x)i | {z } Var + Sesgo2 z }| { (Kh⊛f)(x)−f(x)2. (2.8) Como se ha ido mencionando, el ECM es un criterio puntual, pues evalúa la bondad de ajuste de nuestro estimador para cada x∈Sop(f) . Es de interés denir medidas de error global, 22 3. Métodos de selección del parámetro ventana Partiendo de la hipótesis de que f sigue una distribución Normal, se tiene que R(f′′) = ZR f′′(x)2dx=3 8 1 √πσ5, (3.2) donde σ denota la desviación típica de la distribución. Esta expresión sigue sin ser calculable en la práctica porque depende de la desviación típica teórica, σ . Para solventar esto, basta con sustituir ese valor por una estimación del mismo obtenida a partir de la muestra, que denotaremos en general por ˆσ. Así, sustituyendo esto y el valor de R(f′′) dado en (3.2) en la expresión de la ventana óptima (3.1), se obtiene la regla del pulgar: hSilv =8√πR(K) 3µ2(K)2nˆσ. (3.3) Este estimador, ˆσ , puede ser la desviación típica, el rango intercuartílico, o el mínimo entre ambos. Esta última es la elección más recomendada, pues evita el sobresuavizado y funciona bien con las densidades unimodales y bimodales. 3.2. Regla plug-in Ahora pasamos a introducir el método comunmente conocido como el método de Sheather y Jones, pues fue introducido en (Sheather y Jones, 1991). Se trata de un método plug-in, estas son una familia de métodos que se basan en la idea de pinchar ( plugging-in ) estimaciones de las cantidades desconocidas que aparecen en la fórmula de la ventana óptima dada en (3.1). Para aligerar las expresiones que necesitamos manejar, vamos a comenzar introduciendo la siguiente notación: Rf(r)=ZR f(r)(x)2dx, R f(r)= (−1)rZR f(2r)(x)f(x)dx y ψ(r)=ZR f(r)(x)f(x)dx, (3.4) siendo (r) la notación empleada para la derivada de orden r . Este método de selección del parámetro ventana se basa, como la regla del pulgar de Silverman anteriormente detallada, en la ventana óptima que minimiza el ECMIA , pero en lugar de suponer que f sigue una distribución Normal, lo que hace es sustituir R(f′′) por un estimador. Para ello tenemos en cuenta que: R(f′′) = ZR f′′(x)2dx=ZR f(4)(x)f(x)dx=ψ(4). (3.5) 3.2. Regla plug-in 23 Así, siguiendo las indicaciones de (Wand y Jones, 1994), podemos reescribir la expresión de la ventana óptima como: hECMIA =R(K) µ2(K)2ψ(4)n1 5 . (3.6) En este punto, sustituimos ψ(4) por un estimador, en este caso por su estimador tipo núcleo 3 , b ψ(4) g=1 nPn i=1 ˆ f(4) g(Xi) = 1 n2Pn i=1 Pn j=1 L(4) g(Xi−Xj) con L función núcleo y g parámetro ventana, obteniendo: ˆ h="R(K) µ2(K)2b ψ(4) gn#1 5 . (3.7) El problema con esto es que no es una regla completamente automática, pues como podemos observar, la ventana óptima depende ahora de la elección de una segunda ventana, g , denominada ventana piloto. Por tanto, podemos, de manera análoga a lo realizado en el Capítulo 2 para el estimador tipo núcleo de la densidad, obtener las expresiones de los distintos errores y, a partir de sus versiones asintóticas, la expresión de la ventana óptima para la estimación de ψ(4) . No realizamos esa labor en este trabajo puesto que no aporta gran innovación de desarrollo con respecto a lo ya visto, más allá de un incremento de complejidad en la notación. Sin embargo, puede verse en (Wand y Jones, 1994) que: gECMA ="k!L(r)(0) −µk(L)ψ(r+k)n#1 r+k+1 , (3.8) siendo L una función núcleo sobre la que se necesita incrementar el orden de las condiciones de regularidad en función del valor de k . De esta manera, particularizando para L=K tendríamos, gECMA ="2K(4)(0) −µ2(K)ψ(6)n#1 7 , (3.9) que depende de ψ(6) , cantidad de nuevo desconocida. Con un procedimiento recursivo, esta cantidad se podría volver a estimar, pinchando por ejemplo un estimador tipo núcleo, y su ventana óptima dependería de ψ(8) y así sucesivamente, pues la ventana óptima para estimar ψ(r) con métodos tipo núcleo va a depender de ψ(r+2) . Llegados a este punto, la propuesta de (Wand y Jones, 1994) consiste en parar la segunda iteración y reemplazar la estimación tipo núcleo de ψ(6) por una estimación paramétrica, esto es, estimar ψ(6) =RRf(6)(x)f(x)dx asumiendo que f sigue una distribución Normal (con parámetros estimados a partir de la muestra, como ya se hacía en la regla del pulgar). Así, se obtendría un valor para la ventana piloto g , que se sustituye en (3.7), y se obtiene hSJ : 3 Los estimadores tipo núcleo de las derivadas de la función de densidad se construyen de modo muy intuitivo. Aunque exceda los objetivos del presente trabajo, puede verse detallado en (Wand y Jones, 1994). 24 3. Métodos de selección del parámetro ventana hSJ ="R(K) µ2(K)2ˆ ψ4(g)n#1 5 . (3.10) Esta regla plug-in en dos pasos es una de las más empleadas en la literatura por su buen comportamiento con una amplia variedad de modelos. 3.3. Validación cruzada Para nalizar con los métodos de selección de ventana que exponemos en este trabajo, vamos a detallar el procedimiento de validación cruzada ( V C ). Este selector se basa fundamentalmente en el desarrollo del ECI denido en (2.9): ECI hˆ fnhi=ZRˆ fnh(x)−f(x)2dx =ZRˆ fnh(x)2−2ˆ fnh(x)f(x) + f(x)2dx =ZR ˆ fnh(x)2dx−ZR 2ˆ fnh(x)f(x)dx+ZR f(x)2dx =ZR ˆ fnh(x)2dx−2ZR ˆ fnh(x)f(x)dx+ZR f(x)2dx, (3.11) donde observamos que el último término no depende de h . Por tanto, minimizar (3.11) es equivalente a minimizar la siguiente expresión: ZR ˆ fnh(x)2dx−2ZR ˆ fnh(x)f(x)dx. (3.12) En (3.12), el primero de los sumandos se puede calcular a partir de la muestra, mientras que este no es el caso del segundo, puesto que depende de f . Para solventar esto, hay que tener en cuenta que RRˆ fnh(x)f(x)dx=Ehˆ fnh(X)i , es decir, es la media de la variable aleatoria transformada por el estimador tipo núcleo. Por tanto, el problema se reduce a estimar una media. La primera idea sería \ Ehˆ fnh(x)i=1 nPn i=1 ˆ fnh(Xi) , sin embargo, esta propuesta presenta un problema, y es que estamos usando los datos en dos ocasiones, tanto en la construcción del estimador de la densidad, ˆ fnh , como en su evaluación, lo que produce imprecisiones en la estimación. La forma habitual para solventar este problema es lo que se conoce en inglés como leave one out , es decir, dejar un dato fuera. Esto signica que para evaluar en Xi , el correspondiente estimador de la densidad lo calcularemos con toda la muestra salvo el dato i-ésimo. Formalmente escribimos: \ Ehˆ fnh(x)i=1 n n X i=1 ˆ f−i nh(Xi), (3.13) 3.3. Validación cruzada 25 donde ˆ f−i nh =1 n−1Pj=iKh(x−Xj) . Así, se dene la función de validación cruzada como: V C(h) = Zˆ fnh(x)2dx−21 n n X i=1 ˆ f−i nh(Xi), (3.14) y el selector de ventana es el obtenido con la minimización de esta función: hV C =arg m´ın h∈R+V C(h). (3.15) Se sabe que este selector de ventana, a diferencia de la regla del pulgar, tiende a infrasuavizar las estimaciones, por lo que suele ser empleado en la literatura en estimaciones que, por conocimiento experto o nociones previas, se intuye que proceden de modelos más complejos. En el capítulo siguiente exponemos los resultados del estudio de simulación comparativo que hemos llevado a cabo para analizar el comportamiento de los tres métodos presentados. 26 3. Métodos de selección del parámetro ventana Capítulo 4 Estudio de simulación Nuestro objetivo en este estudio de simulación es comparar diferentes selectores de ventana entre sí, de manera que podamos valorar el error cometido por cada uno de ellos en la estimación de la densidad. Vamos a centrarnos en los selectores que hemos visto en el Capítulo 3, es decir, la propuesta de Silverman, la de Sheather y Jones, y la de validación cruzada. Cuando se están evaluando metodologías que atañen a la función de densidad, un documento muy empleado en la literatura es (Marron y Wand, 1992). En este trabajo se especican 15 mixturas de normales, cuyas densidades cubren sino todas, la mayor parte de las posibles características de un modelo para la densidad (simetría, asimetría, unimodalidad, multimodal, densidades suaves, densidades con picos...). Estos modelos suelen por tanto, ser empleados para la evaluación de métodos relativos a la densidad. En este trabajo vamos a utilizar cuatro de esos modelos: la densidad Normal estándar, la unimodal con curtosis, la bimodal asimétrica y la comunmente conocida como garra, que corresponden a la primera, cuarta, octava y décima en el artículo original, números que mantendremos para evitar confusiones en la nomenclatura de los modelos. En , podemos encontrar estos modelos en la librería nor1mix . En la Tabla 4.1 se resumen los cuatro modelos que nosotros vamos a utilizar. El Modelo 1 se corresponde con la densidad Normal estándar, una densidad unimodal y simétrica respecto del origen. El Modelo 4 y Modelo 8 son mixturas de dos distribuciones normales, la principal diferencia entre ellas reside en que el Modelo 4 tiene una única moda, mientras que el 8 consta de dos modas. El Modelo 10 por su parte, se corresponde a una densidad con cinco modas. En la Figura 4.1 presentamos las correspondientes funciones de densidad, que consideramos cubren un abanico suciente de características para abordar el objetivo de nuestro estudio de simulación. 27 28 4. Estudio de simulación Nomenclatura Expresión Modelo 1 (Normal estándar) N(0,1) Modelo 4 (Unimodal con curtosis) 2 3N(0,1) + 1 3N0,(1 10)2 Modelo 8 (Bimodal asimétrica) 3 4N(0,1) + 1 4N3 2,(1 3)2 Modelo 10 (Garra) 1 2N(0,1) + P4 i=0 1 10N(i 2−1,(1 10)2) Tabla 4.1: expresiones analíticas de los Modelos 1, 4, 8 y 10 extraídos de (Marron y Wand, 1992). (a) Modelo 1: Normal estándar (b) Modelo 4: unimodal con curtosis (c) Modelo 8: bimodal asimétrica (d) Modelo 10: garra Figura 4.1: Funciones de densidad de los Modelos 1, 4, 8 y 10 de (Marron y Wand, 1992). Para realizar el estudio de simulación hemos elaborado un código en , el cual se puede consultar en el Anexo I del presente trabajo. Notemos que, por simplicidad, en todo el estudio se usa la función núcleo Normal, aunque un procedimiento análogo podría seguirse con el núcleo de Epanechnikov o cualquier otro. A pesar de que a nivel teórico el núcleo de Epanechnikov es el óptimo, hemos decidido emplear el Normal porque es el más utilizado en la literatura, y la diferencia en términos de eciencia entre funciones núcleo es muy pequeña. 29 El estudio se realiza para cuatro tamaños muestrales, n= 50,100,200 y 500 , para los que hemos generado aleatoriamente M= 1000 muestras usando procedimientos de Monte Carlo. Para cada una de las muestras calculamos el selector teórico óptimo visto en (2.20) y los tres selectores de ventana detallados en el Capítulo 3. Hemos programado las funciones h.Silverman y h.SJ para los selectores de Silverman y de Sheather y Jones respectivamente, mientras que el procedimiento de validación cruzada se ha implementado usando la función bw.bcv del paquete stats de . A continuación, para estas muestras, analizamos el comportamiento de los tres selectores de ventana calculando distintas medidas de error. Denotaremos, respectivamente, por e1 y por e2 a la media y a la desviación típica del error cuadrático integrado, denido en (2.9), e1=1 M M X k=1 ECI(ˆ hk), e2=1 M M X k=1 (ECI(ˆ hk)−e1)2, donde ECI(ˆ hk) denota el ECI del estimador tipo núcleo de la densidad calculado con la ventana ˆ h que corresponde en cada caso y para cada muestra k∈ {1, ..., M}. Además, no solo nos interesa la estimación en sí misma, sino también cuánto de lejos o cerca está cada selector de la ventana teórica óptima. En este sentido calculamos también la media y la desviación típica de la diferencia entre la ventana óptima teórica y el selector que corresponda en cada caso: e3=1 M M X k=1 (hECMIA −ˆ hk) y e4=1 M M X k=1 (hECMIA −ˆ hk)−e32, donde recordamos que hECMIA es la ventana óptima teórica denida en (2.20), que varía para cada modelo, pero no con cada muestra (aunque sí con el tamaño muestal). Las medidas e1 y e2 sirven para para valorar el comportamiento de los selectores en cuanto al ajuste de la estimación con respecto a la densidad teórica, mientras que las medidas e3 y e4 valoran la diferencia entre el parámetro escogido por cada método de selección y el óptimo teórico que persiguen, tal y como detallamos en el Capítulo 3. Como se puede ver en la Tabla 4.2 para el Modelo 1, en donde se han marcado en gris los menores valores para cada escenario y criterio de error, los errores observados disminuyen al aumentar el tamaño muestral para todos los procedimientos, tanto en media como en desviación típica. Este comportamiento es el razonable y esperable para cualquier procedimiento que funciona correctamente, pues al aumentar el tamaño muestral, estamos incrementando la información. En este modelo sencillo, los tres selectores funcionan de manera similar, aunque destaca el método de validación cruzada, que incluso con tamaños muestrales pequeños, tiene valores muy 30 4. Estudio de simulación Modelo 1 Teórica Silverman SJ VC n= 50 e1 0.18808 0.22016 0.24471 0.19905 e2 0.15382 0.16760 0.19259 0.15212 e3 -3.16497 3.33129 -3.48689 e4 - 6.31035 8.47856 5.42367 n= 100 e1 0.10906 0.11864 0.12820 0.11437 e2 0.07998 0.08402 0.09123 0.08219 e3 - 1.84759 1.76136 -2.79653 e4 - 3.54184 5.33694 3.25326 n= 200 e1 0.06839 0.07206 0.07714 0.07084 e2 0.04839 0.04903 0.05285 0.04933 e3 -1.08172 1.23446 -2.25542 e4 - 2.18961 3.68541 2.18375 n= 500 e1 0.03567 0.03676 0.03834 0.03687 e2 0.02331 0.02350 0.02448 0.02423 e3 -0.51010 0.51675 -1.65295 e4 -1.22339 2.18754 1.46697 Tabla 4.2: media y desviación típica del ECMI y de la diferencia entre la ventana óptima y los selectores medidos, e1 a e4 , para el Modelo 1 (mutiplicado por 102 ). cerca de los valores teóricos óptimos que se persiguen, a excepción del mayor tamaño muestral ( n= 500 ) en donde el selector de Silverman es el que presenta un mejor comportamiento, tal y como cabía esperar para el modelo Normal estándar. En cuanto a las medidas e3 y e4 , se puede observar como los selectores propuestos por Silverman y Sheater y Jones obtienen en media valores inferiores que la ventana óptima, infrasuavizando la densidad de los datos, frente al selector por validación cruzada que sobresuaviza (en media) la estimación de la densidad de los datos. El selector por validación cruzada es el que mayor diferencia en valor absoluto alcanza en media, mientras que en desviación típica es el selector de Sheather y Jones, que sigue al de validación cruzada en media. La menor diferencia tanto en media como en desviación típica, todo en valor absoluto, la consigue el selector de Silverman. Cabe destacar que al aumentar el tamaño muestral tanto en media como en desviación típica se reducen sus valores. En la Tabla 4.3 se muestran las medidas de error correspondientes al Modelo 4. En este caso, el selector con menor error medio, e1 , es el dado por Sheater y Jones. Podemos ver como Silverman y validación cruzada distan bastante del óptimo, aunque este efecto parece corregirse 31 Modelo 4 Teórica Silverman SJ VC n= 50 e1 1.70010 2.92887 2.28989 4.44436 e2 0.76464 1.15899 1.13525 0.64362 e3 - -16.26108 -9.78686 -34.42257 e4 - 7.65648 6.64838 6.14960 n= 100 e1 0.99215 2.32190 1.42027 3.91228 e2 0.42711 0.90617 0.74401 0.70285 e3 - -13.70533 -6.32907 -29.16260 e4 - 5.01927 3.39886 6.03745 n= 200 e1 0.56979 1.85411 0.81109 3.04998 e2 0.22821 0.67412 0.41692 1.19725 e3 - -11.93858 -4.19317 -21.84400 e4 - 3.28612 1.73525 9.53687 n= 500 e1 0.29255 1.32571 0.37893 0.64565 e2 0.12046 0.44687 0.19592 0.93246 e3 - -9.88094 -2.43869 -3.62009 e4 - 1.97104 0.78283 6.87838 Tabla 4.3: media y desviación típica del ECMI y de la diferencia entre la ventana óptima y los selectores medidos, e1 a e4 , para el Modelo 4 (multiplicado por 102 ). en la validación cruzada para tamaños muestrales grandes. En términos de variabilidad, e2 , tanto Silverman como Sheather y Jones parecen más estables que la validación cruzada. Con respecto a la medida e3 , vemos que en este modelo es siempre negativa, lo que indica que los selectores tienden a tomar valores de ventana mayores que la óptima. Globalmente, en lo que respecta al valor de la ventana el selector de Sheather y Jones tiene un mejor comportamiento En la Tabla 4.4 tenemos, análogamente a las anteriores, un resumen para el Modelo 8. En este caso, de los tres selectores propuestos, el selector con menor error cuadrático medio, e1 , es el de Sheather y Jones, salvo para n= 50 , que sería el de Silverman. En desviación típica, e2 , para tamaños muestrales pequeños el menor error se alcanza con el selector de validación cruzada, mientras que para tamaños muestrales más grandes son los tres selectores muy similares, siendo en este caso el menor valor el de Silverman. En cuanto a las diferencias entre ventanas con respecto a la ventana teórica, medidas e3 y e4 , podemos observar como los tres selectores sobresuavizan la estimación de la densidad, pues la diferencia entre la ventana óptima teórica y los selectores es negativa, lo que indica que estos últimos son mayores que la óptima. El selector que menor diferencia tiene en valor absoluto es el de Sheather y Jones, 38 5. Aplicación a datos reales Figura 5.1: estimaciones tipo núcleo de la densidad para los selectores Silverman (trazo negro y continuo), Sheather y Jones (trazo azul y discontinuo) y validación cruzada (trazo rojo y de puntos). Arriba, de izquierda a derecha, las variables mdhc , mdht y abajo la variable ctr . mdhc , mdht y ctr , utilizando el núcleo Normal y los selectores de ventana de Silverman, de Sheater y Jones y validación cruzada. Cabe destacar que, dado que los valores ventana son muy similares las estimaciones son prácticamente idénticas, dichos valores obtenidos por los distintos selectores en cada muestra se resumen en la Tabla 5.1. Si analizamos con un poco de detalle los valores de la Tabla 5.1 vemos algunos datos curiosos: por ejemplo, para ninguna de las variables la ventana de validación cruzada es la que toma menor valor, cuando sabemos por lo comentado en el Capítulo 4 y en la literatura, que es un selector que tiende a infrasuavizar y por consiguiente sus valores tienden a ser pequeños. Si bien es cierto, la diferencia entre unos selectores y otros es tan pequeña que no puede considerarse relevante. Capítulo 6 Conclusiones y trabajo futuro A lo largo de este trabajo hemos revisado el estimador tipo núcleo de la densidad y sus propiedades, tanto a nivel teórico como a través de un completo estudio de simulación. El estimador tipo núcleo tiene dos elementos clave, la función núcleo y el parámetro ventana. Hemos visto como el primero de los elementos tiene un efecto escaso en el resultado de la estimación. Aunque a nivel teórico el núcleo de Epanechnikov resulta ser óptimo, la diferencia en la aplicación con otras funciones núcleo es muy pequeña, por lo que en la práctica, y por cuestiones de facilidad de manejo, la función núcleo más empleada sea probablemente la Normal. Sin embargo, el segundo de los elementos, el parámetro ventana, si hemos visto que presenta una gran inuencia en el resultado de la estimación. La selección de este parámetro ha sido un problema por tanto de gran interés en la literatura estadística y que ha derivado en numerosos procedimientos. En el presente trabajo hemos analizado el comportamiento de tres de los más relevantes: el método del pulgar de Silverman, el método plug-in de Sheather y Jones y el método de validación cruzada. En el Capítulo 4 realizamos un estudio de simulación para valorar el comportamiento de estos selectores. A raíz de los resultados obtenidos, y como ya hemos ido detallando y justicando a lo largo del trabajo, el selector que presenta un mejor comportamiento en términos generales es el de Sheather y Jones. Por último, después de aplicar los métodos de estimación de la densidad observamos que en nuestros datos los tres selectores son muy parecidos en todas las variables de estudio, el máximo diámetro horizontal cardíaco, el máximo diámetro horizontal torácico y el ratio cardio-torácico. Como propuestas de trabajo futuro o posibles mejoras de la presente memoria, se podría ampliar el estudio de simulación mediante el estudio de otros selectores del parámetro ventana que existen en la literatura, como son el método de validación cruzada insesgado o el selector 39 40 6. Conclusiones y trabajo futuro de ventana bootstrap, basado, como su nombre indica, en procedimientos de remuestreo. . Otra opción de ampliación sería utilizar los doce modelos restantes de (Marron y Wand, 1992), o incluso introducir nuevos modelos que no sean mixturas de Normales, y así estudiar el comportamiento de los selectores a nivel todavía más general, o en casos especialmente patológicos. Otra cuestión que encontramos relevante a raíz de la aplicación a datos reales, aunque excedería con creces los objetivos del presente trabajo sería la estimación multivariante de la densidad. Esto es, en lugar de estar interesados en la función de densidad de una variable aleatoria (real de variable real) podríamos estar interesados en la función de densidad de un vector aleatorio d-dimensional. En el caso de nuestros datos podría ser interesante estimar conjuntamente la densidad del máximo diámetro horizontal torácico y del máximo diámetro horizontal cardíaco. Este problema supone aumentar la complejidad del proceso de estimación con respecto al que abordamos en este trabajo, y además, el parámetro ventana pasa a ser una matriz de ventanas de dimensión dxd , con el consiguiente incremento de complejidad en el problema de su selección. Anexo I Código para el estudio de simulación # Selector ventana de Silverman h.Silverman<-function(x,K){ n<-length(x) sd<-sqrt(var(x)) ri<-diff(quantile(x, probs = c(0.25, 0.75)))/(qnorm(0.75)-qnorm(0.25)) # el rango intercuartilico muestral tambien se puede obtener con IQR(x) sigma.hat<-min(sd,ri) if(K=="Gauss"){ mu_2<- 1 Rk<-1/(2*sqrt(pi)) } else if (K=="Epa"){ mu_2<-1/5 Rk<-3/5 } mu_2<-mu_2^2 h.Silverman<-((8*sqrt(pi)*Rk)/(3*mu_2*n))^(1/5)*sigma.hat return(h.Silverman) } #Selector de ventana de Sheather y Jones h.SJ<-function(x,K){ n<-length(x) if(K=="Gauss"){ K6<-function (x) {exp((-x^2)/2)*(x^6-15*x^4+45*x^2-15)/(sqrt(2*pi))} K4<-function (x) {exp((-x^2)/2)*(x^4-6*x^2+3)/(sqrt(2*pi))} 41 42 I. Código para el estudio de simulación mu_2<-1 Rk<-1/(2*sqrt(pi)) } else if (K=="Epa"){ K6<-function (x) {0} K4<-function (x) {0} # la derivada segunda ya es una constante mu_2<-1/5 Rk<- 3/5 #mu_2 y Rk sacados de la tabla B.2 de Kernel Smoothing } #1 º estimamos fi8 utilizando la regla del pulgar sd<-sd(x) #sqrt(var(x)) ri<-diff(quantile(x, probs = c(0.25, 0.75)))/(qnorm(0.75)-qnorm(0.25)) sigma.hat<-min(sd,ri) fi8<-105/(32*sqrt(pi)*sigma.hat^9) #ahora calculamos la ventana óptima de fi6 (que es la que depende de fi8) g1<-((2*K6(0))/(-mu_2*fi8*n))^(1/9) #ahora puedo calcular la estimación de fi6 con la ventana óptima g1 a<-outer(x,x,"-")/g1 #matriz nxn con (Xi-Xj)/g1 Keval<-matrix(K6(a),nc=n,nrow=n) #evalúo la matriz en la # derivada correspondiente de K fi6<-(1/(n^2*g1^7))*sum(Keval) #calculamos ahora la ventana optima g2 de fi4 g2<-((2*K4(0))/(-mu_2*fi6*n))^(1/7) #calculamos fi4 a<-outer(x,x,"-")/g2 #matriz nxn con (Xi-Xj)/g1 Keval2<-matrix(K4(a),nc=n,nrow=n) #evalúo la matriz en la # derivada correspondiente de K fi4<-(1/(n^2*g2^5))*sum(Keval2) h.SJ<-(Rk/((mu_2)^2*fi4*n))^(1/5) # fifthroot <- function(x)sign(x)*abs(x)^(1/5) # h.SJ <- fifthroot(h.SJ) return(h.SJ) } #Estimador tipo núcleo f.kernel<-function(x1,x0,h,K){ n<-length(x1) 43 a<-outer(x0,x1,FUN="-")/h # me da una matriz dxn con d la longitud de x0 if(K=="Gauss"){Kf=K_gauss} else if (K=="Epa"){Kf=K_epa} Keval<-matrix(Kf(a),nc=n,nrow=length(x0)) # funcion nucleo de a yy=rowSums(Keval)/(n*h) # promedio los nucleos por filas return(data.frame('xrej'=x0,'y'=yy)) } K_gauss <- function(x){1/sqrt(2*pi)*exp(-x^2/2)} #con -infty<x<infty # Estudio de simulación library(nor1mix) n<-50 # n=50,100,200,500 M<-1000 #numero de muestras con cada tamaño muestral n.rej<-100 #longitud de la rejilla #creamos unas matrices para guardar los datos de las ventanas y los errores mat.err <- matrix(NA, ncol=4, nrow=M) colnames(mat.err) <- c("Teorica", "Silverman","SJ","CV") mat.e1 <- matrix(NA, ncol=4, nrow=1) colnames(mat.e1) <- c("Teorica", "Silverman","SJ","CV") mat.e2 <- matrix(NA, ncol=4, nrow=1) colnames(mat.e2) <- c("Teorica", "Silverman","SJ","CV") mat.h <- matrix(NA, ncol=4, nrow=M) colnames(mat.h) <- c("Teorica", "Silverman","SJ","CV") mat.dif <- matrix(NA,ncol=3,nrow=M) colnames(mat.dif) <-c("Teorica-Silverman","Teorica-SJ","Teorica-CV") mat.e3 <- matrix(NA, ncol=3, nrow=1) colnames(mat.e3) <- c("Silverman","SJ","CV") mat.e4 <- matrix(NA, ncol=3, nrow=1) colnames(mat.e4) <- c("Silverman","SJ","CV") #hacemos una función con la ventana teórica dada por el ECMIA ya que esta #solo va a depender del tamaño muestral y del modelo teórico h.teorica<- function (n,modelo) { Rk <- 1/(2*sqrt(pi)) #estoy suponiendo que K es el núcleo gaussiano mu2 <- 1 #estoy suponiendo que K es el núcleo gaussiano 44 I. Código para el estudio de simulación if (modelo=="MW.nm1"){ #N(0,1) Rf <- 3/(8*sqrt(pi)) #sacado de Wand & Jones h <- (Rk/(n*mu2^2*Rf))^0.2 } else { if (modelo=="MW.nm4") { f4 <- function (x) {(2/(3*sqrt(2*pi)) * exp(-x^2/2) * (-1+x^2)+1/(3*sqrt(pi*2)) * exp(-50*x^2) * (-10^3+10^5*x^2))^2} Rf <- integrate(f4,-Inf,Inf) h <- (Rk/(n*mu2^2*Rf$value))^0.2 } else { if(modelo=="MW.nm8") { f8 <- function(x) {(0.75* 1/sqrt(2*pi) * exp(-x^2/2)*(-1+x^2) + 0.25*1/sqrt(2*pi) * exp(-9/2*(x-1.5)^2) * (-27+243*(x-1.5)^2))^2} Rf <- integrate(f8,-Inf,Inf) h <- (Rk/(n*mu2^2*Rf$value))^0.2 } else{ if(modelo=="MW.nm10") { f10 <- function(x) { mu <- c(-1,-0.5,0,0.5,1) sig <- 0.1 res <- 0.5 * (x^2-1) * 1/sqrt(2*pi) * exp(-x^2/2) + 0.1/sig^3 * (((x-mu[1])/sig)^2-1) * 1/sqrt(2*pi) * exp(-0.5*((x-mu[1])/sig)^2) + 0.1/sig^3 * (((x-mu[2])/sig)^2-1) * 1/sqrt(2*pi) * exp(-0.5*((x-mu[2])/sig)^2) + 0.1/sig^3 * (((x-mu[3])/sig)^2-1) * 1/sqrt(2*pi) * exp(-0.5*((x-mu[3])/sig)^2) + 0.1/sig^3 * (((x-mu[4])/sig)^2-1) * 1/sqrt(2*pi) * exp(-0.5*((x-mu[4])/sig)^2) + 0.1/sig^3 * (((x-mu[5])/sig)^2-1) * 1/sqrt(2*pi) * exp(-0.5*((x-mu[5])/sig)^2) return(res^2) } Rf <- integrate(f10,-Inf,Inf) h <- (Rk/(n*mu2^2*Rf$value))^0.2 } } } } return(h) } 45 ############################ MODELO 1 #calculamos entonces la densidad teórica: h.teor <- h.teorica(n,modelo="MW.nm1") rejilla <- seq(qnorMix(0.01,MW.nm1),qnorMix(0.99,MW.nm1),len=n.rej) d.teor <- dnorMix(rejilla,obj= nor1mix::MW.nm1) #densidad set.seed(12345) #fijamos la semilla for (m in 1:M){ #generamos una muestra x <- nor1mix::rnorMix(n, obj = nor1mix::MW.nm1) # Aplicamos cada uno de los selectores: mat.h[m,1] <- h.teor mat.h[m,2] <- h.Silverman(x,K="Gauss") mat.h[m,3] <- h.SJ(x,K="Gauss") mat.h[m,4] <- bw.bcv(x) # Calculamos la estimacion con cada uno de los selectores d.teor.est <- f.kernel(x,x0=rejilla,h=mat.h[m,1],K="Gauss") d.silverman <- f.kernel(x,x0=rejilla,h=mat.h[m,2],K="Gauss") d.sj <- f.kernel(x,x0=rejilla,h=mat.h[m,3],K="Gauss") d.cv <- f.kernel(x,x0=rejilla,h=mat.h[m,4],K="Gauss") # Calculamos medida de error ECI (de cada muestra) mat.err[m,1] <- mean((d.teor.est[,2]-d.teor)^2) mat.err[m,2] <- mean((d.silverman[,2]-d.teor)^2) mat.err[m,3] <- mean((d.sj[,2]-d.teor)^2) mat.err[m,4] <- mean((d.cv[,2]-d.teor)^2) # Calculamos el error de las ventanas (de cada muestra) mat.dif[,1] <- mat.h[,1]-mat.h[,2] mat.dif[,2] <- mat.h[,1]-mat.h[,3] mat.dif[,3] <- mat.h[,1]-mat.h[,4] } # Calculamos la media (ECMI) y sd de los errores (de las M muestras) 46 I. Código para el estudio de simulación mat.e1[1] <- mean(mat.err[,1]) mat.e1[2] <- mean(mat.err[,2]) mat.e1[3] <- mean(mat.err[,3]) mat.e1[4]<- mean(mat.err[,4]) mat.e2[1] <- sd(mat.err[,1]) mat.e2[2] <- sd(mat.err[,2]) mat.e2[3] <- sd(mat.err[,3]) mat.e2[4] <- sd(mat.err[,4]) # Calculamos media y sd de los errores (diferencias de ventanas) mat.e3[1] <- mean(mat.dif[,1]) # media de h.silverman - h.teorica mat.e3[2] <- mean(mat.dif[,2]) # media de h.sj - h.teorica mat.e3[3] <- mean(mat.dif[,3]) mat.e4[1] <- sd(mat.dif[,1]) mat.e4[2] <- sd(mat.dif[,2]) mat.e4[3] <- sd(mat.dif[,3]) #boxplot de los errores(ECI) boxplot(mat.err*100) ######################## MODELO 4 #calculamos entonces la densidad teórica: h.teor <- h.teorica(n,modelo="MW.nm4") rejilla <- seq(qnorMix(0.01,MW.nm4),qnorMix(0.99,MW.nm4),len=n.rej) d.teor <- dnorMix(rejilla,obj= nor1mix::MW.nm4) set.seed(12345) #fijamos la semilla for (m in 1:M){ # generamos una muestra x <- nor1mix::rnorMix(n, obj = nor1mix::MW.nm4) # Aplicamos cada uno de los selectores: mat.h[m,1] <- h.teor mat.h[m,2] <- h.Silverman(x,K="Gauss") mat.h[m,3] <- h.SJ(x,K="Gauss") 47 mat.h[m,4] <- bw.bcv(x) # Calculamos la estimacion con cada uno de los selectores d.teor.est <- f.kernel(x,x0=rejilla,h=mat.h[m,1],K="Gauss") d.silverman <- f.kernel(x,x0=rejilla,h=mat.h[m,2],K="Gauss") d.sj <- f.kernel(x,x0=rejilla,h=mat.h[m,3],K="Gauss") d.cv <- f.kernel(x,x0=rejilla,h=mat.h[m,4],K="Gauss") # Calculamos medida de error ECI (de cada muestra) mat.err[m,1] <- mean((d.teor.est[,2]-d.teor)^2) mat.err[m,2] <- mean((d.silverman[,2]-d.teor)^2) mat.err[m,3] <- mean((d.sj[,2]-d.teor)^2) mat.err[m,4] <- mean((d.cv[,2]-d.teor)^2) # Calculamos tambien el error de las ventanas (de cada muestra) mat.dif[,1] <- mat.h[,1]-mat.h[,2] mat.dif[,2] <- mat.h[,1]-mat.h[,3] mat.dif[,3] <- mat.h[,1]-mat.h[,4] } # Calculamos la media (ECMI) y sd de los errores (de las M muestras) mat.e1[1] <- mean(mat.err[,1]) mat.e1[2] <- mean(mat.err[,2]) mat.e1[3] <- mean(mat.err[,3]) mat.e1[4]<- mean(mat.err[,4]) mat.e2[1] <- sd(mat.err[,1]) mat.e2[2] <- sd(mat.err[,2]) mat.e2[3] <- sd(mat.err[,3]) mat.e2[4] <- sd(mat.err[,4]) # Calculamos media y sd de los errores (diferencias de ventanas) mat.e3[1] <- mean(mat.dif[,1]) # media de h.silverman - h.teorica mat.e3[2] <- mean(mat.dif[,2]) # media de h.sj - h.teorica mat.e3[3] <- mean(mat.dif[,3]) mat.e4[1] <- sd(mat.dif[,1]) mat.e4[2] <- sd(mat.dif[,2]) 54 II. Código para la aplicación a datos reales a<-outer(x0,x,FUN="-")/h # me da una matriz dxn con d la longitud de x0 if(K=="Gauss"){Kf=K_gauss} else if (K=="Epa"){Kf=K_epa} Keval<-matrix(Kf(a),nc=n,nrow=length(x0)) # funcion nucleo de a yy=rowSums(Keval)/(n*h) # promedio los nucleos por filas return(data.frame('xrej'=x0,'y'=yy)) } # Representamos la estimación de x1 con la densidad real f<-f.kernel(x1,x0=seq(200,400,len=100),h=h.Silverman(x1,K="Gauss"),K="Gauss") plot(f$xrej,f$y,type="l",lwd=2,ylab="",xlab="") lines(density(x1),col="red") # Programamos ahora las ventanas # Regla del pulgar de Silverman h.Silverman<-function(x,K){ n<-length(x) sd<-sqrt(var(x)) ri<-diff(quantile(x, probs = c(0.25, 0.75)))/(qnorm(0.75)-qnorm(0.25)) # el rango intercuartilico muestral tambien se puede obtener con IQR(x) sigma.hat<-min(sd,ri) if(K=="Gauss"){ mu_2<- 1 Rk<-1/(2*sqrt(pi)) } else if (K=="Epa"){ mu_2<-1/5 Rk<-3/5 } mu_2<-mu_2^2 h.Silverman<-((8*sqrt(pi)*Rk)/(3*mu_2*n))^(1/5)*sigma.hat return(h.Silverman) } # Sheather & Jones h.SJ<-function(x,K){ n<-length(x) if(K=="Gauss"){ K6<-function (x) {exp((-x^2)/2)*(x^6-15*x^4+45*x^2-15)/(sqrt(2*pi))} 55 K4<-function (x) {exp((-x^2)/2)*(x^4-6*x^2+3)/(sqrt(2*pi))} mu_2<-1 Rk<-1/(2*sqrt(pi)) } else if (K=="Epa"){ K6<-function (x) {0} K4<-function (x) {0} # la derivada segunda ya es una constante mu_2<-1/5 Rk<- 3/5 } #1 º estimamos fi8 utilizando la regla del pulgar sd<-sd(x) ri<-diff(quantile(x, probs = c(0.25, 0.75)))/(qnorm(0.75)-qnorm(0.25)) sigma.hat<-min(sd,ri) fi8<-105/(32*sqrt(pi)*sigma.hat^9) #ahora calculamos la ventana óptima de fi6 (que es la que depende de fi8) g1<-((2*K6(0))/(-mu_2*fi8*n))^(1/9) #estimación de fi6 con la ventana óptima g1 a<-outer(x,x,"-")/g1 #matriz nxn con (Xi-Xj)/g1 Keval<-matrix(K6(a),nc=n,nrow=n) fi6<-(1/(n^2*g1^7))*sum(Keval) #ventana optima g2 de fi4 g2<-((2*K4(0))/(-mu_2*fi6*n))^(1/7) #calculamos fi4 a<-outer(x,x,"-")/g2 #matriz nxn con (Xi-Xj)/g1 Keval2<-matrix(K4(a),nc=n,nrow=n) fi4<-(1/(n^2*g2^5))*sum(Keval2) h.SJ<-(Rk/((mu_2)^2*fi4*n))^(1/5) return(h.SJ) } 56 II. Código para la aplicación a datos reales Referencias Ibarrola, R. V., y Pérez, A. G. (2006). Principios de inferencia estadística . Universidad Nacional de Educación a Distancia. Marron, J. S., y Wand, M. P. (1992). Exact mean integrated squared error. The Annals of Statistics , 20 (2), 712736. Parzen, E. (1962). On Estimation of a Probability Density Function and Mode. The Annals of Mathematical Statistics , 33 (3), 1065  1076. Pearson, K., y Henrici, O. M. F. E. (1895). X. contributions to the mathematical theory of evolution. skew variation in homogeneous material. Philosophical Transactions of the Royal Society of London. (A.) , 186 , 343-414. Rosenblatt, M. (1956). Remarks on Some Nonparametric Estimates of a Density Function. The Annals of Mathematical Statistics , 27 (3), 832  837. Sheather, S. J., y Jones, M. C. (1991). A reliable data-based bandwidth selection method for kernel density estimation. Journal of the Royal Statistical Society. Series B (Methodological) , 53 (3), 683690. Shiraishi, J., Katsuragawa, S., Ikezoe, J., Matsumoto, T., Kobayashi, T., Komatsu, K.-i., ... Doi, K. (2000). Development of a digital image database for chest radiographs with and without a lung nodule: receiver operating characteristic analysis of radiologists' detection of pulmonary nodules. American Journal of Roentgenology , 174 (1), 7174. Silverman, B. W. (1986). Density estimation for statistics and data analysis . Chapman Hall, London. Sogancioglu, E., Murphy, K., Calli, E., Scholten, E. T., Schalekamp, S., y Van Ginneken, B. (2020). Cardiomegaly detection on chest radiographs: Segmentation versus classication. IEEE Access , 8 , 94631-94642. van Ginneken, B., Stegmann, M. B., y Loog, M. (2006a). Scr database: Segmentation in chest radiographs. Descargado 2022-03-01, de https://www.isi.uu.nl/Research/Databases/ SCR/index.php van Ginneken, B., Stegmann, M. B., y Loog, M. (2006b). Segmentation of anatomical structures in chest radiographs using supervised methods: a comparative study on a public database. 57 58 Referencias Medical Image Analysis , 10 (1), 1940. Vélez Ibarrola, R., y Hernández Morales, V. (1995). Cálculo de probabilidades 2. Universidad Nacional de Educación a Distancia. Wand, M. P., y Jones, M. C. (1994). Kernel smoothing . CRC press.