🕷️ Crawler Inspector

# ¿Cómo se estima una densidad? #### José R. Berrendero - [Preparación de los datos](https://verso.mat.uam.es/~joser.berrendero/blog/kernel.html#preparaci%C3%B3n-de-los-datos) - [Una alternativa al histograma a partir de la definición de densidad](https://verso.mat.uam.es/~joser.berrendero/blog/kernel.html#una-alternativa-al-histograma-a-partir-de-la-definici%C3%B3n-de-densidad) - [Estimadores del núcleo](https://verso.mat.uam.es/~joser.berrendero/blog/kernel.html#estimadores-del-n%C3%BAcleo) - [Cómo se representan en ggplot() estos estimadores](https://verso.mat.uam.es/~joser.berrendero/blog/kernel.html#c%C3%B3mo-se-representan-en-ggplot-estos-estimadores) - [Selección del parámetro de suavizado](https://verso.mat.uam.es/~joser.berrendero/blog/kernel.html#selecci%C3%B3n-del-par%C3%A1metro-de-suavizado) - [Valores numéricos](https://verso.mat.uam.es/~joser.berrendero/blog/kernel.html#valores-num%C3%A9ricos) ### Preparación de los datos Para empezar bajamos y preparamos los datos que nos van a servir de ejemplo. Los datos corresponden a longitud, peso y contenido de mercurio en tejidos de una muestra de peces capturados en varias estaciones de dos ríos (Lumber y Wacamaw). ``` enlace <- "http://matematicas.uam.es/~joser.berrendero/datos/mercurio.txt" mercurio <- read.table(enlace, header = TRUE) %>% mutate(RIO = recode(RIO, "0" = "Lumber", "1" = "Wacamaw")) %>% mutate(RIO = as.factor(RIO)) %>% mutate(ESTACION = as.factor(ESTACION)) glimpse(mercurio) ``` ``` ## Rows: 171 ## Columns: 5 ## $ RIO <fct> Lumber, Lumber, Lumber, Lumber, Lumber, Lumber, Lumber, Lu... ## $ ESTACION <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1... ## $ LONG <dbl> 47.0, 48.7, 55.7, 45.2, 44.7, 43.8, 38.5, 45.8, 44.0, 40.4... ## $ PESO <dbl> 1616, 1862, 2855, 1199, 1320, 1225, 870, 1455, 1220, 1033,... ## $ CONC <dbl> 1.60, 1.50, 1.70, 0.73, 0.56, 0.51, 0.48, 0.95, 1.40, 0.50... ``` El histograma de los datos se obtendría así. He comentado para que se entiendan mejor las opciones: ``` ggplot(mercurio) + geom_histogram(aes(x = CONC, y = ..density..), # y = ..density.. normaliza para área = 1 bins = 10, # número de intervalos col = 'black', # color del borde fill = 'olivedrab4') + # color del interior labs(x = "Concentración (ppm)", y = NULL) # cambia etiquetas de los ejes ```  ### Una alternativa al histograma a partir de la definición de densidad Sabemos que las áreas bajo la función de densidad corresponden a las probabilidades de que la variable tome valores en el correspondiente intervalo. Por lo tanto, si \\(h\\approx 0\\) se tiene la siguiente aproximación: \\\[ \\mbox{P}(x-h\\leq X\\leq x+h) = \\int\_{x-h}^{x+h} f(t) dt \\approx 2h f(x). \\\] Usando esta observación construimos un estimador de la función de densidad tomando \\(h\>0\\) pequeño y sustituyendo la probabilidad anterior por la proporción de datos que caen dentro del intervalo \\((x-h, x+h)\\). El estimador resultante es \\\[ \\hat{f}(x) = \\frac{1}{2h} \\frac{\\\#\\{i:\\, \|x-X\_i\| \< h\\}}{n}. \\\] En la expresión anterior \\(\\\#\\) denota el número de elementos de un conjunto. La elección de \\(h\\) sería análoga a la elección del número de intervalos en un histograma. Más adelante comentamos con más detalle en qué se debe basar. ### Estimadores del núcleo Si definimos la función \\(K(x) = (1/2) \\mathbb{I}\_{\\{\|x\|\\leq 1\\}}\\), donde \\(\\mathbb{I}\_A\\) denota la función indicatriz sobre \\(A\\) (que vale 1 en \\(A\\) y 0 en el complementario de \\(A\\)), entonces podemos reescribir el estimador anterior como: \\\[ \\hat{f}(x) = \\frac{1}{nh} \\sum\_{i=1}^nK\\left(\\frac{x-X\_i}{h}\\right)=\\frac{1}{n} \\sum\_{i=1}^n \\frac{1}{h}K\\left(\\frac{x-X\_i}{h}\\right). \\\] Fijado un valor de \\(x\\), contamos a cuántos intervalos centrados en \\(X\_i\\) y de radio \\(h\\) pertenece \\(x\\). Asignamos un valor \\(1/(2h)\\) si esto ocurre y cero en caso contrario, y entonces promediamos los valores. Nótese que la función \\(K\\) que hemos definido corresponde a la función de densidad de una v.a. uniforme en el intervalo \\(\[-1,1\]\\). Si sustituimos la función indicatriz por una función más suave obtendremos estimadores que se asemejan más a una función de densidad típica. Así se obtienen los estimadores del núcleo de la densidad: \\\[ \\hat{f}(x) = \\frac{1}{nh} \\sum\_{i=1}^nK\\left(\\frac{x-X\_i}{h}\\right), \\\] donde \\(h\>0\\) es el **parámetro de suavizado** y \\(K\\) es el núcleo (que suele ser una función de densidad simétrica). Si, por ejemplo, \\(K(x)=e^{-x^2/2}/\\sqrt{2\\pi}\\) es un núcleo gaussiano (la función de densidad de una v.a. normal estándar) entonces los valores \\\[ \\frac{1}{h}K\\left(\\frac{x-x\_i}{h}\\right) \\\] corresponderán a las funciones de densidad de v.a. con distribución normal de media \\(x\_i\\) y desviación típica \\(h\\). El estimador del núcleo es entonces el promedio de los valores de estas densidades en cada punto \\(x\\). En el gráfico siguiente (tomado de [wikipedia](https://en.wikipedia.org/wiki/Kernel_density_estimation)), las curvas rojas corresponden a \\\[ \\frac{1}{nh}K\\left(\\frac{x-x\_i}{h}\\right),\\ \\ i = 1,\\ldots,n \\\] La suma de estas curvas nos da el estimador del núcleo, \\(\\hat{f}(x)\\).  ### Cómo se representan en ggplot() estos estimadores En `ggplot2` usamos el comando `geom_density` para representar los estimadores del núcleo. El comando `geom_density` usa por defecto un núcleo gaussiano. ``` ggplot(mercurio) + geom_density(aes(x = CONC), size = 1.2, # grosor de línea fill = 'olivedrab4') + labs(x = "Concentración (ppm)", y = " ") ```  Mediante el argumento `kernel` se pueden usar otros núcleos diferentes al gaussiano. Por ejemplo, el núcleo rectangular \\(K(x) = (1/2) \\mathbb{I}\_{\\{\|x\|\\leq 1\\}}\\) que da el estimador que vimos anteriormente proporciona un resultado menos suave, más parecido al histograma: ``` ggplot(mercurio) + geom_density(aes(x = CONC), size = 1.2, fill = 'olivedrab4', kernel = 'rectangular') + labs(x = "Concentración (ppm)", y = " ") ```  Las opciones de elección de núcleo son `kernel = c("gaussian", "epanechnikov", "rectangular", "triangular", "biweight", "cosine", "optcosine")`. Como técnica descriptiva, una ventaja de los estimadores del núcleo respecto a los histogramas es que son más apropiados para superponerlos en el mismo gráfico. A continuación se comparan los estimadores correspondientes a los peces de cada uno de los dos ríos: ``` ggplot(mercurio) + geom_density(aes(x = CONC, linetype = RIO), size = 1.2) + labs(x = "Concentración (ppm)", y = " ") ```  Una manera elegante de comparar estimadores es mediante un *gráfico de cordillera* (usando [el paquete `ggridges`](https://cran.r-project.org/web/packages/ggridges/vignettes/introduction.html)): ``` ggplot(mercurio) + ggridges::geom_density_ridges(aes(x = CONC, y = RIO, fill = RIO), size = 1.2) + labs(x = "Concentración (ppm)", y = "Río") ``` ``` ## Picking joint bandwidth of 0.273 ```  ### Selección del parámetro de suavizado La selección del parámetro de suavizado \\(h\\) es un problema difícil e importante ya que el valor de este parámetro determina drásticamente las propiedades del estimador. Su papel es análogo al del número de clases de un histograma. Hay que elegir un valor de \\(h\\) tal que el sesgo y la varianza del estimador estén bien equilibrados. Compárense por ejemplo las situaciones de la siguiente figura:  La densidad que se desea estimar es la normal estándar (línea discontinua) y se han aplicado las fórmulas anteriores a un conjunto de muestras de esta población. En el panel superior se ha aplicado un valor de \\(h\\) demasiado grande. El resultado es que los estimadores se parecen bastante entre sí (tienen poca varianza) pero presentan desviaciones sistemáticas de la verdadera densidad (tienen bastante sesgo). Por ejemplo, en la parte central siempre infraestiman. El panel central describe lo que ocurre para valores demasiado pequeños de \\(h\\). El promedio de los estimadores corresponde más o menos a la densidad verdadera (sesgo pequeño) pero cada estimador es muy diferente del resto (varianza muy alta). Para obtener un buen resultado debemos elegir un valor de \\(h\\) que produzca un equilibrio adecuado entre sesgo y varianza (la situación del panel inferior, en la que se ha usado un método sofisticado y automático de selección de \\(h\\)). Existe teoría (basada en aproximaciones al error cuadrático medio integrado de los estimadores) según la cual un valor equilibrado de \\(h=h\_n\\) tiene que ser del orden de \\(n^{-1/5}\\), donde \\(n\\) es el número de datos. Esto significa que debe decrecer con el tamaño muestral pero bastante lentamente. También parece lógico que el valor de \\(h\\) dependa de la escala. Concretamente, el valor que usa `geom_density()` es \\(h\_n = 0.9 \\hat{\\sigma} n^{-1/5}\\), donde \\(\\hat{\\sigma}=\\min\\{s,\\mbox{RI}/1.34\\}\\) (con \\(s\\) la desviación típica y \\(\\mbox{RI}\\) el rango intercuartílico de los datos). El valor \\(\\hat{\\sigma}\\) es una estimación de la desviación típica de la población con un cierto grado de robustez, gracias a la intervención del rango intercuartílico. El valor por defecto de \\(h\\) se puede cambiar usando el argumento `adjust`. Por ejemplo, para usar un valor de \\(h\\) que sea la cuarta parte del que se usa por defecto: ``` ggplot(mercurio) + geom_density(aes(x = CONC), size = 1.2, adjust = 0.25) + labs(x = "Concentración (ppm)", y = " ") ```  Claramente, este es un valor de \\(h\\) demasiado pequeño, con muchas variaciones espurias, que probablemente no corresponden a la densidad original. ### Valores numéricos Hemos visto cómo se representan los estimadores del núcleo con `geom_density()`. Si lo que queremos son los valores numéricos del estimador \\(\\hat{f}(x)\\) se pueden obtener a través del comando `density()`, que subyace a todas las representaciones gráficas anteriores. Por ejemplo, para obtener el valor del estimador en cinco puntos equiespaciados en el intervalo \\(\[1,2\]\\) habría que usar: ``` estimador_nucleo <- density(mercurio$CONC, n = 5, from = 1, to = 2) estimador_nucleo$x # puntos donde se calcula la densidad ``` ``` ## [1] 1.00 1.25 1.50 1.75 2.00 ``` ``` estimador_nucleo$y # densidad en esos puntos ``` ``` ## [1] 0.5214496 0.4341464 0.3918171 0.2951467 0.1816800 ``` [Enlace a *Caminos aleatorios*](http://caminosaleatorios.wordpress.com/) [Enlace a mi página web de la UAM](https://verso.mat.uam.es/~joser.berrendero/index.html)