1. Teoría

1.1 Definición matemática

Dado un conjunto de datos \(x=(x_1,x_2,..., x_n)\), la función de distribución de densidad \(f(x)\) puede aproximarse utilizando una Estimación de Densiada de Kernel tal que:

\[\widehat f(x)=\frac{1}{n}\sum_{i=1}^{n}K_{h}(x-x_{i})=\frac{1}{nh}\sum_{i=1}^{n}K\frac{(x-x_{i})}{h}\]

  • n: Es el número de datos (observaciones). Cada uno de ellos es el centro sobre el que se coloca un kernel.
  • h: es el ancho de banda. Controla cuanto se expande la influencia de cada observación. Si se emplea como kernel una distribución normal, equivale a la desviación típica. Este es el valor más determinante a la hora de ajusta KDE, ya que condiciona el nivel de sobreajuste.
  • K: es el Kernel, una función que define la forma y la distribución de la influencia (peso) que se asocian a cada observación.

Kernel

SESGO

Analizaremos el error de estimación del KDE, donde asumimos que X1,··· ,Xn son IID de una función de densidad desconocida de la funcion p. En el problema de estimación de la densidad, el parámetro de interés es p, la verdadera función de densidad.

Para simplificar el problema, supongamos que nos enfocamos en un punto dado X0 y queremos analizar la calidad de nuestro estimador \[\widehat{p}_n(x_o)\] Primero analizamos el sesgo de KDE: \[E(\widehat{p}_n(x_0))-p(x_o)=E(\frac{1}{nh}\sum_{i=1}^{n}K(\frac{X_i-x_o}{h}))-p(x_o)\] \[=\frac{1}{h}E(\sum_{i=1}^{n}K(\frac{X_i-x_o}{h}))-p(x_o)\] \[=\frac{1}{h}\int K(\frac{x-x_o}{h})p(x)dx-p(x_o)\] Ahora realizamos un cambio de variable \[y=\frac{x-x_o}{h}\] \[E(\widehat{p}_n(x_0))-p(x_o)=\int K(\frac{x-x_o}{h})p(x)\frac{dx}{h}-p(x_o)\] \[=\int K(y)p(x_o+hy)dy-p(x_o)\] Dado que \[ x=x_o+hy\] Ahora por expansion de Taylor dado que h es pequeño, \[p(x_o+hy)=p(x_o)-hy.p'(x_o)+\frac{1}{2}h^2y^2p''(x_o)+o(h^2)\] Teniendo en cuenta que \[o(h^2)\] significa que es un termino de orden mas pequeño en comparacion con \[h^2\] cuando h tiende a o nuevamente con el sesgo obtenemos: \[E(\widehat{p}_n(x_0))-p(x_o)=\int K(y)p(x_o+hy)dy-p(x_o)\] \[=\int K(y)[p(x_o)+hy.p'(x_o)+\frac{1}{2}h^2y^2p''(x_o)+o(h^2)]dy-p(x_o)\]

\[=\int K(y)p(x_o)dy+\int K(y)hy.p'(x_o)dy+\int K(y)\frac{1}{2}h^2y^2p''(x_o)dy+o(h^2)-p(x_o)\] \[=p(x_o)\int K(y)dy+hp'(x_o)\int yK(y)dy+\frac{1}{2}h^2p''(x_o)\int y^2K(y)dy+o(h^2)-p(x_o)\] \[=p(x_o)+\frac{1}{2}h^2p''(x_o)\int y^2K(y)dy+o(h^2)-p(x_o)\] \[=\frac{1}{2}h^2p''(x_o)\int y^2K(y)dy+o(h^2)\] \[=\frac{1}{2}h^2p''(x_o)u_k+o(h^2)\] donde \[u_k=\int y^2K(y)dy\] Es decir es sesgo va ser igual a: \[sesgo(\widehat{p}_n(x_0))=\frac{1}{2}h^2p''(x_o)u_k+o(h^2)\]

VARIANZA

Para el análisis de varianza, podemos obtener un límite superior utilizando un cálculo sencillo:

\[Var_{\hat p_n(x_0)} = var \left(\frac{1}{nh}\sum_{i=1}^{n}K \left(\frac{X_i-x_0}{h}\right)\right)\] \[ = \frac{1}{nh^2} Var\left(K \left(\frac{X_i-x_0}{h}\right)\right)\] \[ \leq \frac{1}{nh^2} E\left(K^2 \left(\frac{X_i-x_0}{h}\right)\right)\] \[ = \frac{1}{nh^2} \int K^2 \left(\frac{X_i-x_0}{h}\right)p(x)dx\] \[\text { Usando } y=\frac{x-x_0}{h} \text{ y } dy=dx/h :\] \[ = \frac{1}{nh} \int K^2(y)p(x_0+hy)dy \] \[ = \frac{1}{nh} \int K^2(y)\left[p(x_0+hyp'(x_0)+o(h)\right])dy\] \[ = \frac{1}{nh} \left(p(x_0) \int K^2(y)dy+o(h)\right)dy\] \[ = \frac{1}{nh} p(x_0) \int K^2(y)dy+o\left(\frac{1}{nh}\right)\] \[ = \frac{1}{nh} p(x_0) \sigma ^2_ {K}+o\left(\frac{1}{nh}\right)\] donde \(\sigma^2_{K}=\int K^2 (y)dy\).Por lo tanto, la varianza se reduce a una tasa \(O \left (\frac{1}{nh}\right)\) cuando \(n \rightarrow \infty\) y \(h \rightarrow 0\). Un hecho interesante de la varianza es que en el punto donde el valor de la densidad es grande,la varianza también lo es.

Error cuadrático medio

Una vez conocido el sesgo y la varianza, se puede obtener el Error cuadrático medio (Mean Squared Error MSE) de KDE:

\[MSE(\hat p_n(x_0)) = bias^2(\hat p_n(x_0)) + Var\] \[ = \frac{1}{4}|p''(x_0)|^2\mu^2_K+\frac{1}{nh}p(x_0)\sigma^2_K+o(h^4)+o \left(\frac{1}{nh}\right)\] \[ = O(h^4)+O\left(\frac{1}{nh}\right)\] Los términos \[\frac{1}{4}|p''(x_0)|^2\mu^2_K+\frac{1}{nh}p(x_0)\sigma^2_K\] se conocen como el error cuadrático medio asintótico (Asymptotic Mean Square Error AMSE). En el KDE, suavizar el ancho de banda \(h\) es algo que se puede elegir. Por lo tanto, el ancho de banda \[h\] que minimiza el AMSE es

\[h_{opt}(x_0)^2_K =\left(\frac{4}{n}\cdot \frac{p(x_0)}{|p''(x_0)|^2} \frac{\sigma^2_K}{\mu^2_K}\right)^\frac{1}{5}=C_1 \cdot n^{\frac{-1}{5}}\]

Con esta elección de suavizar el ancho de banda conduce a un MSE de tasa:

\[MSE_{opt}(\hat p_n (x_0)) =O\left(n^{\frac {-4}{5}}\right)\]

1.2 Teoremas

1.2.1 Teorema

Para \(h > 0\),

\(E [\widehat{J}(h)] = E [J(h)]\)

Tambien,

\(\widehat{J}(h)\approx\frac{1}{hn^2}\sum_i\sum_jK^*\left( \frac{X_i-X_j}{h}\right)+\frac{2}{nh}K(0)\)

donde \(K^*(x)=K^{(2)}(x)-2K(x)\) y \(K^{(2)}(z)=\int{K(z-y)}K(y)dy\).

En particular, si \(K\) es \(N(0,1)\) un nucleo gaussiano, despues \(K^{(2)}(z)\) es la densidad \(N(0,2)\)

A continuacion, elegimos el ancho de banda \(h_n\) que minimiza \(\widehat{J}(h).^3\) Una justificacion para este metodo viene dado por el Teorema de Stone.

1.2.2 Teorema de Stone

Supongamos que \(f\) está acotada.

Sea \(\widehat{f}(h)\) el estimador de Kernel con ancho de banda \(h\) y sea \(h_n\) el ancho de banda elegido por validación cruzada.

Entonces,

\({\frac{\int({f(x)-\widehat{f}h_n(x)})^2dx}{inf_h\int({f(x)-\widehat{f}h_n(x)})^2dx}}\underrightarrow{P}1\)

1.3 Regresión no paramétrica

Considere los pares de puntos \((x_1, Y_1), … , (x_n, Y_n)\) relacionado como \(Y_i= r(x_i) + \epsilon_i\)

  • En donde \(E(\epsilon_i) = 0\). Las \(x_i\) se escriben en minúsculas ya que se tratan como valores fijos.

Se puede hacer esto debido al hecho que, en la regresión, solamente la media de \(Y\) condicionada en \(x\) es la que interesa.

Se quiere estimar la función de regresión \(r(x) = E (Y \shortmid X = x)\). Existen muchos estimadores de regresión no paramétricos, en la mayoría de estos se estima \(r(x)\) tomando algún tipo de promedio ponderado de los \(Y_i\), dado mayor peso a los puntos cerca de \(x\).

Una versión popular es la del estimador kernel de Nadaraya – Watson. La cual se define como:

\[\widehat{r}(x)=\sum_{i=1}^{n}w_i(x)Y_i\]

En donde \(K\) es un kernel y los pesos \(w_i(x)\) están dados por:

\[w_i(x)=\frac{K(\frac{x-x_i}{h})}{\sum_{j=1}^{n}K(\frac{x-x_j}{h})}\]

La forma de este estimador proviene de estimar primero la densidad conjunta \(f(x, y)\) utilizando la estimación de la densidad kernel:

\[r(x)=\mathbb{E} (Y|X=x)=\int yf(y|x)dy=\frac{\int yf(x,y)dy}{\int f(x,y)dy}\]

1.3.1 Teorema 1


Suponga que \(V(\epsilon i) = \sigma^2\). El riesgo del estimador kernel de Nadaraya-Watson es:

\[R(\widehat{r}_n;r)\approx \frac{h^4}{4} (\int x^2K^2(x)dx)^4\int(r''(x)+2r'(x)\frac{f'(x)}{f(x)})^2dx+\int\frac{\sigma ^2 \int K^2(x)dx}{nhf(x)}dx \]

El ancho de banda óptimo disminuye a una tasa \(n ^{-1/5}\) y con esta elección el riesgo disminuye a una tasa \(n^{-4/5}\). En la práctica, para elegir el ancho de banda \(h\) minimizamos la puntuación de validación cruzada.

\[\widehat{J}(h)=\sum_{i=1}^{n}(Y_i-\widehat{r}_{-1}(x_i))^2\] Donde \(\widehat{r}_{-1}\) es el estimador que obtenemos al omitir la i-ésima variable. Afortunadamente, existe una fórmula abreviada para calcular \(\widehat{J}\)

1.3.2 Teorema 2


\(\widehat{J}\) puede ser escrito como:

\[\widehat{J}=\sum_{i=1}^{n}(Y_i-\widehat{r}(x_i))^2\frac{1}{(1-\frac{K(0)}{\sum_{j=1}^{n}K(\frac{x_i-x_j}{h})})^2}\]

2. Cálculo

Estimación de Densidad de Kernel (Kernel Density Estimator KDE) es uno de los métodos más famosos para la estimación de densidad. KDE amplia la idea del histograma, cada observación aumenta la densidad de probabilidad de la zona donde se encuentra, pero lo hace de forma que las contribuciones se agrupen creando una curva continua y suave.

2.1 Selección del ancho de banda

El ancho de banda es cruacial a la hora de estimar una función de densidad mediante el método KDE. Si su valor es muy bajo, se genera sobreajuste y la función resultante estará demasiado influenciada por el “ruido” de los datos. Si su valor es muy elevado, la función resultante no será capaz de aprender la distribución subyacente Aunque no hay forma de saber de antemano cuál es el valor óptimo de ancho de banda, existen varias estrategias para identificar buenos valores.

Ancho de banda estrecho

Ancho de banda suave

2.2 Reglas empíricas

  • Scott’s rule: \(h \approx 1.06 \cdot \widehat\sigma n^{\frac{-1}{5}}\)
  • Silverman´s rule: \(h =0.9 \cdot min( \widehat\sigma,IQR/1.35)n^\frac {-1}{5}\)_

Regla de Silverman

Validación cruzada

Requiere de mayor tiempo de computación pero es útil para cualquier tipo de distribución. Validación cruzada es la opción disponible en Scikit learn, utilizando como métrica de comparación en log-likelihood.

Es importante descartar que, cuando se dispone de pocos datos, conviene utilizar one leave out crossvalidation.

2.3 Tipo de Kernel

El kernel es la función que determina cómo se distribuye la influencia de cada observación, por lo tanto, puede tener un impacto notable en la estimación de la función de densidad resultante. Aunque en la gran mayoría de casos se emplea un kernel gaussiano (distribución normal), existen otras posibilidades. Las implementadas en scikit learn son:

  • Gaussian: asigna los pesos siguiendo la distribución normal con una desviación estándar equivalente al ancho de banda.

\(K(x;h) \propto exp \left(- \frac{x^2}{2h^2}\right)\)

  • Triangular: \(K(u)=(1-|u|)\)

  • Box: \(K(u)=\frac{1}{2}\)

  • Triweight: \(K(u)=\frac{35}{32}(1-u^2)^3\)

Algunos tipos de Kernel

(kernels <- eval(formals(density.default)$kernel))
## [1] "gaussian"     "epanechnikov" "rectangular"  "triangular"   "biweight"    
## [6] "cosine"       "optcosine"
plot (density(0, bw = 1), xlab = "",
      main = "Otros tipos de Kernel")
for(i in 2:length(kernels))
  lines(density(0, bw = 1, kernel =  kernels[i]), col = i)
legend(1.5,.4, legend = kernels, col = seq(kernels),
       lty = 1, cex = .5, y.intersp = 1)

2.4 Ejercicio

Para el siguiente ejercicio manual, se tienen dos observaciones (1 y 3); con un ancho de banda de 1.

Datos

\(x_1=1\)

\(x_2=3\)

\(n=2\)

\(h=1\)

\[\widehat f(x)=\frac{1}{n}\sum_{i=1}^{n}K_{h}(x-x_{i})=\frac{1}{nh}\sum_{i=1}^{n}K\frac{(x-x_{i})}{h}\]

hist(c(1,3), breaks = 1, col = "#D1EEEE", main = "Ejemplo")

Kernel Triangular proviniente de la Distribución Triangular, \(K(u)=(1-|u|)\)

Evaluando:

\[K_1=\frac{1}{nh}\frac{(1-|u|)}{1}\in (0;2)\]

\[u=x-x_1\] \[K_1=\frac{1}{(2)(1)}\frac{(1-|x-x_1|)}{1}\in (0;2)\]

Se procede a integrar considerando \(x_1=1\)

\(\int_0^2 (1-|x-1|)dx=\int_0^1 (1-(1-x))dx+ \int_1^2 (1-(x-1))dx=\)

\(\frac{x^2}{2}|_{0}^{1}+(2x-\frac{x^2}{2})|_{1}^{2}\)

\(\frac{1}{2}-(0)+(4-2)-(2-\frac{1}{2})=\frac{1}{2}+2-\frac{3}{2}=1\)

Ahora con \(K_2\)

\[K_2=\frac{1}{(2)(1)}\frac{(1-|x-3|)}{1}\in (2;4)\]

\(\int_2^4 (1-|x-3|)dx=\int_2^3 (1-(1-3))dx+ \int_3^4 (1-(x-3))dx= 1\)

Entonces:

\[\int \widehat f(x)dx=\frac{1}{n}\int K_1dx+\frac{1}{n}\int K_2dx\] Con \(n=2\):

\(\widehat f(x)=\frac{1}{2}*1+\frac{1}{2}*1=1\)

Función de densidad estimada:

\[ \widehat f(x)=\left\{ \begin{array}{ll} 0; & \text{si }x < 0 \\ \frac{1}{2}K_1*1-|x-1| & \text{si }x\in (0;2)\\ \frac{1}{2}K_2*1-|x-3| & \text{si }x\in (2;4)\\ 0; & \text{si }x > 4 \end{array} \right. \]

\[ \widehat f(x)=\left\{ \begin{array}{ll} \frac{1}{2}*0; & \text{si }x < 0 \\ \frac{1}{2}*x; & \text{si }x\in (0;1)\\ \frac{1}{2}*(2-x) & \text{si }x\in (1;2)\\ \frac{1}{2}*(x-2) & \text{si }x\in (2;3)\\ \frac{1}{2}*(4-x) & \text{si }x\in (3;4)\\ \frac{1}{2}*(0) ; & \text{si }x > 4 \end{array} \right. \]

kf <- density(c(1,3), kernel = "triangular")

x11()
par(mfrow=c(1,2))
plot (density(c(1,3)), col="red",main="F.D Original")
plot(kf,main="K. Triangular")

3. Ejemplos

3.1 Ejemplo de la estimación Kernel

Para el ejemplo se utilizó la base pretederminada de mtcars.

library(ggplot2)
library(tidyverse)
library(dbplyr)
base <- mtcars
ggplot(data = base,
       mapping = aes(x = base$mpg)) +
  geom_histogram(bins = 9,col='black', fill='#00FFFF', alpha=0.9) +
  labs(title = '1. Histograma ',
       x = 'Kilometraje',
       y = 'Frecuencia',
       caption = 'Rstudio')

La Función de densidad de una variable aleatoria, generalmente desconocida. Para obtener una aproximación general se realiza en primera instancia un histograma de frecuencias.

El histograma es una representación gráfica de una variable en forma de barras, donde la superficie de cada barra es proporcional a la frecuencia de los valores representados, esto con la finalidad de dar una imagen del comportamiento de la variable.

ggplot(base, aes(x = base$mpg)) + 
       geom_density(color = "black", 
                    fill = "#7AC5CD",  
                    alpha = 0.9)+
  labs(title = '2. Densidad',
       x = 'Kilometraje',
       y = 'Densidad',
       caption = 'Rstudio')

La primera aproximación con respecto a la Función de Densidad, es que no sigue una distribución normal.

La curva es suave, sin embargo, tiene dos bultos.

Descripción

La función genérica kdensity {stats} calcula las estimaciones de la densidad del núcleo. Su método por defecto lo hace con el kernel y el ancho de banda dados para observaciones univariantes.

Argumentos

x: Vector numérico que contiene los datos.

bw: El ancho de banda de suavizado que se utilizará. Los kernels se escalan de tal manera que esta es la desviación estándar del kernel de suavizado.

bw: también puede ser una cadena de caracteres que da una regla para elegir el ancho de banda. Véase bw.nrd. El valor por defecto, “nrd0”, se ha mantenido por razones históricas y de compatibilidad, más que como una recomendación general.

kernel: La función del núcleo. Puede elegirse de la lista de núcleos incorporados o hacerse a medida.

El argumento kernel de la función density usa por defecto un kernel gaussiano (kernel = “gaussian”), pero existen otros tipos de kernel disponibles, como “rectangular”, “triangular”, “epanechnikov”, “biweight”, “cosine” y “optcosine”. La seleción dependerá de tus datos, pero el valor por defecto es el más recomendado.

start: Inicio paramétrico. Puede elegirse de la lista de inicios paramétricos incorporados o hacerse a medida. Ver parametric_starts() para más detalles.

na.rm:
Lógico; si es TRUE, los NAs serán eliminados de x.

plot_start:

lógico; si es VERDADERO, traza el inicio paramétrico en lugar de la estimación de la densidad del núcleo.

library("kdensity")
k1 <- kdensity(base$mpg, start="gaussian", kernel = "gaussian", bw= "nrd0")
k1
## 
## Call:
## kdensity(x = base$mpg, bw = "nrd0", kernel = "gaussian", start = "gaussian")
## 
## Data:      base$mpg (32 obs.)
## Bandwidth: 2.477 ('nrd0')
## Support:   (-Inf, Inf)
## Kernel:    gaussian
## Start:     gaussian

Comparación entre un Kernel Gaussiano - Kernel párametrico

plot(k1, main = "Comparación entre un Kernel Gaussiano - Kernel párametrico", cex=0.8, col="#68228B")
lines(k1, plot_start = TRUE, col="red")
rug(base$mpg)
legend("topright", legend = c("Gaussiano", "E. Parámetrico"),
       lwd = 3, col = c("#68228B", "red"))

x11()
par(mfrow=c(2,3))
plot (density(c(1,3)), col="#68228B",main="F.D Original")
plot (density(base$mpg,kernel="epanechnikov"), col="red", main = "E.epanechnikov")
plot (density(base$mpg,kernel="triangular",bw=0.5), col="#6E8B3D", main = "E. triangular")
plot (density(base$mpg,kernel="cosine"), col="#CD5555", main = "E. Coseno")
plot (density(base$mpg,kernel="rectangular"), col="#CD2990", main = "E. rectangular")
plot (density(base$mpg,kernel="gaussian"), col="#CD2990", main = "E. Gaussiano")

3.2 Ejemplo de regresión no parámetrica

En Rstudio para la estimación de la regresión kernel de Nadaraya-Watson se utiliza el comando ksmooth {stats} - Kernel Regression Smoother.

Las funciones de núcleo usadas en la estimación no paramétrica de la regresión son las mismas que en la densidad.

En el caso particular del ajuste de un polinomio de grado cero, se obtiene el estimador de Nadaraya Watson, o estimador núcleo de la regresión.

Argumentos

x: valores x de entrada. Se admiten vectores largos.

y: entrada de valores y. Se admiten vectores largos.

kernel: el núcleo a utilizar. Puede ser abreviado.

bandwidth: el ancho de banda.

Utilizando la base de datos “cars” de R, que contine las variables “dist” (distancia de parada) y “speed” (velocidad), vamos a realizar la representación gráfica de la regresión kernel realizada con el estimador de Nadaraya–Watson con diferentes parámetros de suavizado.

data(cars)
correlacion <- cor(cars$speed, cars$dist)
correlacion
## [1] 0.8068949

Existe dependencia lineal positiva.

Mientras mayor sea la velocidad mas tardan en detenerse.

R1 <- ksmooth(cars$speed, cars$dist, kernel =  "normal", bandwidth = 5)

head(R1$x, 5)
## [1] 4.000000 4.212121 4.424242 4.636364 4.848485
head(R1$y, 5)
## [1] 7.939636 8.269956 8.639827 9.054622 9.502761
plot(cars$speed, cars$dist)
lines(ksmooth(cars$speed, cars$dist, kernel =  "normal", bandwidth = 2), col = "red")
lines(ksmooth(cars$speed, cars$dist, "normal", bandwidth = 5), col ="blue")

El estimador del parámetro de suavizado h tiene una importancia crucial en el aspecto y propiedades del estimador de función de regresión. Valores pequeños de h dan mayor flexibilidad al estimador y le permiten acercarse a todos los datos observados, pero originan altos errores de predicción (sobre-estimación), valores mas altos de h ofrecerán un menor grado de ajustes a los datos pero predicican mejor, pero si h es demasiado elevado tendremos una falta de ajuste a los datos (sub estimación)

est <- approx(cars, xout = cars$speed)$cars$dist # predictivos
resid <- cars$dist - est # residuales
r.squared <- 1 - sum(resid^2)/sum((cars$dist - mean(cars$dist))^2)
r.squared
## [1] 1

El modelo explica que a mayor velocidad la distancia en que se detiene el auto es mayor, tomando en cuenta que la base de datos es de 1920.

Bibliografía

Chen, Y.-C. (2018). Introduction to Nonparametric Statistics.

Monfor, P. (2010). Estimación núcleo de distribuciones: existencia y propiedades asintoticas del ancho de banda optimo . Badajoz.

Nadayara, & Watson. (s.f.). Estimación núcleo y métodos de posicionamiento de los países. En Nadayara, & Watson.

Parra, F. (s.f.). Estimadores nucleos y polinomios locales.

Wasserman, L. (2013). All of statistics.

Anders Munk-Nielsen. (2016, 11 noviembre). Nonparametric Kernel regression [Vídeo]. YouTube. https://www.youtube.com/watch?v=ncF7ArjJFqM&ab_channel=AndersMunk-Nielsen

Sebastián Gómez Clases. (2021, 10 junio). Kernel, Estimación de función de densidad no paramétrica, II [Vídeo]. YouTube. https://www.youtube.com/watch?v=LyvsQf5_ywU&t=1681s&ab_channel=Sebasti%C3%A1nG%C3%B3mezClases

Webel od. (2018, 24 septiembre). Intro to Kernel Density Estimation [Vídeo]. YouTube. https://www.youtube.com/watch?v=x5zLaWT5KPs&ab_channel=webelod

https://www.icane.es/c/document_library/get_file?uuid=c2e9fff1-72d5-42ad-b391-bacb3ebe9dbe&groupId=10138