La diferencia entre estadística paramétrica y no paramétrica está basada en el conocimiento o desconocimiento de la distribución de probabilidad de la variable que se pretende estudiar.
En otras palabras la estadística paramétrica utiliza cálculos y procedimientos asumiendo que conoce como se distribuye una variable aleatoria \(X\) en particular para determinado análisis. Por el otro lado, la estadística no paramétrica utiliza métodos para conocer como se distribuye un fenómeno en particular para después utilizar metodologías de la estadística paramétrica.
Ambas son ramas de la estadística inferencial clásica.
Tomado: Economipedia (2023)
Estas secciones son tomadas de Wasserman (2006) .
Definición 2.1 Muestra Aleatoria Simple
Una muestra aleatoria simple (m.a.s) es un conjunto de variables aleatorias, independientes e idénticamente distribuidas, obtenidas a partir de la variable aleatoria \(X\) y que se distribuyen igual que la misma.
Definición 2.2 Función de Distribución Empírica
Sea \(X_1,X_2,...,X_n \sim F\), una muestra aleatoria simple (m.a.s) con \(F\) función de distribución real. La Función de Distribución Empírica \(\hat F_n\) es la Función de distribución Acumulada (CDF) tal que pondera cada punto \(X_i\) con \(\frac{1}{n}\) (peso), es decir:
\[\hat F_n(x)=\frac{\sum_{i=1}^n I_{(X_i \leq x)}}{n} \tag{2.1}\]
donde:
\[ I_{(X_i \leq x)}=\begin{cases} 1 & X_i \leq x\\ 0 & X_i > x \end{cases} \]
Ejemplo 2.1
Tomemos una muestra de tamaño 1000, con \(X \sim N(0,1)\)
library(functional)
set.seed(2024)
BD=rnorm(1000,0,1)
EmpiricalCDF<-function(x,sample){
sum(sample <= x, na.rm =TRUE )/length(sample)
}
aux<-function(x){
return(lapply(x,Curry(EmpiricalCDF,sample=BD)))
}
curve(aux,from=-3,to=3,ylab="CDF",col=2,main="CDF Empírica vs Teórica")
curve(pnorm,from=-3,to=3,add=TRUE, col=4)
legend(x="topleft",legend=c("empírica","téorica"),fill=c(2,4))
Para estimar \(P(-1 \leq X \leq 1) \approx \hat F_n(1)-\hat F_n(-1) \approx 0.692\)
(EmpiricalCDF(1,BD)-EmpiricalCDF(-1,BD))[1] 0.692Observe que la aproximación depende de la cantidad de los datos de la muestra para valores con pocas observaciones podría ocurrir que ciertos intervalos de probabilidad sean cero. Siguiendo el ejemplo:
Si estimamos la probabilidad que la variable aleatoria (v.a.) \(X\) esté entre 2.75 y 2.80, el resultado sería 0, según nuestra aproximación, pero el dato teórico sería de 0.0004246329.
Otro caso extremo sería tener pocas observaciones, supongamos que nuestra muestra son solo 100 números, veamos como aumenta la diferencia entre la curva teórica y la curva empírica.
library(functional)
set.seed(2024)
BD2=rnorm(100,0,1)
EmpiricalCDF<-function(x,sample){
sum(sample <= x, na.rm =TRUE )/length(sample)
}
aux<-function(x){
return(lapply(x,Curry(EmpiricalCDF,sample=BD2)))
}
curve(aux,from=-3,to=3,ylab="CDF",col=2,main="CDF Empírica (100) vs Teórica")
curve(pnorm,from=-3,to=3,add=TRUE, col=4)
legend(x="topleft",legend=c("empírica","téorica"),fill=c(2,4))
Si estimamos la probabilidad que la v.a. \(X\) esté entre -1.40 y -1.25 el valor teórico sería de 0.02489311
cat("Valor Teórico: \n" )Valor Teórico: (pnorm(-1.25)-pnorm(-1.40))[1] 0.02489311cat("Valor muestra de 1000: \n" )Valor muestra de 1000: (EmpiricalCDF(-1.25,BD)-EmpiricalCDF(-1.40,BD))[1] 0.02cat("Valor muestra de 100: \n" )Valor muestra de 100: (EmpiricalCDF(-1.25,BD2)-EmpiricalCDF(-1.40,BD2))[1] 0Teorema 2.1 Teorema Esperanza, Varianza y Convergencia
Para cualquier valor fijo de \(x\):
\[E[\hat F_n(x)]=F(x) \tag{2.2}\] \[Var[\hat F_n(x)]=\frac{F(x)(1-F(x))}{n} \tag{2.3}\] \[\hat F_n(x)\xrightarrow{p}F_n(x) \tag{2.4}\]
Teorema 2.2 Teorema Glivenko-Cantelli
Sea \(X_1,X_2,...,X_n \sim F\). Entonces cuando \(n\to \infty\):
\[||\hat F_n(x)-F(x)||_\infty=\sup_x |\hat F_n(x)-F(x)|\xrightarrow{p}0 \tag{2.5}\]
es decir \(\hat F_n(x)\) converge uniformemente a \(F(x)\).
Para Teorema 2.2 se puede usar la ley débil de los Grandes Números para demostrar la convergencia, pero también se puede usar la Ley Fuerte, en ese caso se puede probar que la convergencia es casi segura. Es decir que para cada \(x_i\) existe un conjunto de medida cero \(A_x\) tal que:
\[\lim_{n\to \infty}\hat F_n(x) \xrightarrow{c.s.} F(x)\]
Dada la convergencia en probabilidad y la convergencia uniforme podemos enuciar la desigualdad DKW para estimar intervalo de confianza para \(F\).
Teorema 2.3 Teorema Desigualdad Dvoretzky-Kiefer-Wolfowitz (DKW)
Sea \(X_1,X_2,...,X_n \sim F\) Entonces para cualquier \(\epsilon >0\):
\[P \left( \sup_x |\hat F_n(x)-F(x)|>\epsilon \right) \leq 2e^{-2n\epsilon^2} \tag{2.6}\]
Un intervalo de confianza no paramétrico de \(F\) con un nivel de significancia de \(\alpha\), esta definido por:
\[L(x)=max(\hat F_n(x)-\epsilon_n,0)\] \[U(x)=min(\hat F_n(x)+\epsilon_n,1)\]
con \(\epsilon_n=\sqrt{\frac{1}{2n}\ln(\frac{2}{\alpha})}\)
Por la desigualdad DKW, Teorema 2.3, se tiene que:
\[P \left( L(x) \leq F(x) \leq U(X) \right) \geq 1-\alpha\]
Ejemplo 2.2 Continuación Ejemplo 2.1
Continuando con el ejemplo desarrollado en el documento, para el caso de \(n=1000\), para determinar un intervalo de confianza del 95% (significancia = 5%), nuestro error \(\epsilon_n\) = \(\sqrt{\frac{1}{2(1000)}\ln(\frac{2}{0.05}))} = 0.0018\) y para el caso de \(n=100\) y el mismo nivel de significancia sería de \(0.0184\)
library(functional)
set.seed(2024)
BD2=rnorm(100,0,1)
EmpiricalCDF<-function(x,sample){
sum(sample <= x, na.rm =TRUE )/length(sample)
}
ep<-1/200*log(2/0.05)
aux<-function(x){
return(lapply(x,Curry(EmpiricalCDF,sample=BD2)))
}
Lx<-function(x){
max(aux(x)[[1]]-ep,0)
}
Ux<-function(x){
min(aux(x)[[1]]+ep,1)
}
aLx<-function(x){lapply(x,Lx)}
aUx<-function(x){lapply(x,Ux)}
curve(aux,from=-3,to=3,ylab="CDF",col=2,main="CDF Empírica (100) vs Teórica")
curve(pnorm,from=-3,to=3,add=TRUE, col=4)
curve(aLx,from=-3,to=3,add=TRUE, col=9)
curve(aUx,from=-3,to=3,add=TRUE, col=9)
legend(x="topleft",legend=c("empírica","téorica"),fill=c(2,4))
Definición 2.3 Funcional estadístico
Un funcional estadístico \(T(F)\) es cualquier función de \(F\).
Ejemplo 2.3 Ejemplos Funcionales Estadísticos
Sea una variable aleatoria \(X\) con cdf \(F\), algunos ejemplos de funcionales estadísticos son:
-La media: \(\mu=\int xdF(x)\).
-La varianza: \(\sigma^2=\int (x-\mu)^2dF(x)\)
-mediada: \(m=F^{-1}(1/2)\).
Definición 2.4 Principio plug-in
La estimación plug-in de un parámetro \(\theta=T(F)\), donde \(T\) es un funcional estadístico de \(F\) se define como:
\[\hat \theta = T(\hat F_n)\] el estimador plug-in es un buen estimador cuando sólo conocemos una muestra \(X\), sin embargo este estimador no nos indica nada sobre al precisión de la estimación.
Definición 2.5 Funcional lineal
Si \(T(F)=\int r(x)dF(x)\) para alguna función \(r(x)\) entonces \(T\) se le llama funcional lineal.
Esto dado que si F y G son funciones de distribución satisface:
\[T(aF+G)=aT(F)+T(G)\]
Dado que \(\hat F_n(x)\) es discreta, con un peso de \(1/n\) para cada \(X_i\) de la muestra, tenemos:
Si \(T(F)=\int r(x)dF(x)\) el estimador plug-in del funcional lineal es:
\[T(\hat F_n(x))=\int r(x)d\hat F_n(x)=\frac{1}{n}\sum_{i=1}^nr(X_i) \tag{2.7}\]
Para muchos cálculos requerimos encontrar un estimador para el error estándar ( desviación típica o desviación estándar) \(se\) del estimador plug-in \(T(\hat F_n(x))\), en algunos casos no es fácil de obtener esta medición, por lo que se requiere métodos como bootstrap que veremos más adelante para estimarlo, por el momento asumiremos que conocemos el estimador \(\hat{se}\).
Intervalo de Confianza con base Normal
Si \(T(\hat F_n(x)) \sim N(T(F),\hat{se}^2 )\), el intervalo de confianza con base normal con grado de significancia \(\alpha\) es el intervalo:
\[T(\hat F_n(x)) \pm z_{\alpha/2}\hat{se} \tag{2.8}\]
Ejemplo 2.4 media
Si tomamos \(r(x)=x\) la función identidad obtenemos que \(T(F)=\int r(x)dF(x)=\int xdF(x)= \int xf(x)dx=\mu_x\) en caso continuo o \(\sum xp(X=x)=\mu_x\) en caso discreto. El estimador plug-in de \(\mu\) es
\[\hat\mu=\frac{1}{n}\sum_{i=1}^nr(X_i)=\frac{1}{n}\sum_{i=1}^n(X_i)=\bar X\] En caso que conozcamos la varianza \(\sigma^2\) sabemos que \(se=\sigma/\sqrt{n}\), si no lo conocemos podríamos estimar \(\sigma\) por medio de \(\hat \sigma\), entonces \(\hat{se}=\hat \sigma/\sqrt{n}\).
Ejemplo 2.5 Varianza
Sea \(\sigma^2=T(F)=E[x^2]-(E[x])^2\), estimador plug-in es:
\[\hat \sigma^2=\int x^2 d\hat F_n(x)-\left( \int x d\hat F_n(x) \right)^2\] \[=\frac{1}{n}\sum_{i=1}^n(X_i)^2-\left(\frac{1}{n}\sum_{i=1}^n(X_i)\right)^2\] \[=\frac{1}{n}\sum_{i=1}^n(X_i-\bar X)^2=S^2\]
Ejemplo 2.6 Ejemplo Asimetría (Skewness)
\[\kappa=\frac{E[(x-\mu)^3]}{\sigma^3}\] \[\hat \kappa = \frac{\frac{1}{n}\sum_{i}(X_i-\hat \mu)^3}{\hat \sigma^3}\]
Realizar Ejercicio capitulo 7 Wasserman (2006)
La estimación de las curvas de densidad no paramétricas, utiliza técnicas similares para problemas de regresión no paramétricas que se conocen como curvas de estimación o suavizado de curvas (smoothing). Esta parte es tomada del Wasserman (2006), Aguilar (s.f.) y Ojeda (2014)
Definición 2.6 Definición Error Cuadrático integrado (ISE)
Sea \(g\) una función de densidad desconocida y \(\hat g_n\) un estimador de \(g\). Observe que \(\hat g_n(x)\) varía según los datos de la muestra.
Una función de pérdida de información es el error cuadrático integrado (ISE) y se define como:
\[L(g,\hat g_n)=\int(g(u)-\hat g_n(u))^2du\]
Definición Error Cuadrático Medio integrado (MISE)
Se define como una medida de riesgo al error cuadrático medio integrado (MISE) como:
\[R(g,\hat g_n)=E[L(g,\hat g_n)]\]
Lema 2.1 Riesgo Total Sea \(b(x)\) el sesgo y \(v(x)\) la varianza de \(\hat g_n(x)\) respectivamente para un \(x\) fijo. Entonces:
\[R(g,\hat g_n)=\int b^2(x)dx+\int v(x)dx\]
el Lema 2.1 nos indica que existe una relación entre el sesgo de nuestro estimador y su varianza, situación similar con el Error cuadrático Medio, que se conoce como compensación sesgo-varianza o bais-variance tradeoff. Que el Riesgo o Error Total lo podemos separar como la suma de dos errores el del sesgo y el de la varianza.
Si buscamos una \(\hat g_n(x)\) con poca variabilidad, el cuadrado del sesgo aumenta, esta condición se conoce como subajuste (underfitting), esto ocurre cuando un modelo es muy sencillo y no captura toda la información de la muestra, la curva tiende a ser muy suavizada
Si buscamos una \(\hat g_n(x)\) con un sesgo mínimo, la varianza tiene a aumentar, esta condición se conoce como sobreajuste (overfitting), esto ocurre cuando un modelo captura toda la información de la muestra, donde prácticamente la curva \(\hat g_n(x)\) coincide prácticamente con todos los todos los puntos de la muestra, la curva no es suave.
Un ajuste balanceado o equilibrado, es cuando el Riesgo se minimiza.
A nivel del Riesgo o Error total, se puede graficar el tradeoff de un modelo como el siguiente:
Problema
Dada una m.a.s. de \(N\) observaciones \(X=x_1,x_2,...,x_N\), se quiere estimar una función de densidad \(g\) desconocida.
Una primera aproximación para estimar la función \(g\) es utilizando un histograma.
Realizaremos un ejemplo con los datos de una distribución de asaltos en Estados Unidos, como primer acercamiento realizaremos un histograma.
df<-USArrests
summary(df) Murder Assault UrbanPop Rape
Min. : 0.800 Min. : 45.0 Min. :32.00 Min. : 7.30
1st Qu.: 4.075 1st Qu.:109.0 1st Qu.:54.50 1st Qu.:15.07
Median : 7.250 Median :159.0 Median :66.00 Median :20.10
Mean : 7.788 Mean :170.8 Mean :65.54 Mean :21.23
3rd Qu.:11.250 3rd Qu.:249.0 3rd Qu.:77.75 3rd Qu.:26.18
Max. :17.400 Max. :337.0 Max. :91.00 Max. :46.00 cat("histograma \n")histograma hist(df$Assault, main="Histograma Asaltos EEUU")
En general para la muestra \(X=x_1,x_2,...,x_N\) del problema, se debe hallar un intervalo de soporte \(S=[a,b]\), tal que, \(a \leq x_i \leq b \; \forall x_i \in X\).
El intervalo \(S\) se divide en \(K\) partes iguales (particiones o bins) de longitud \(\frac{b-a}{K}\) tal que:
\[ I_k=\left[ a+(k-1)\left( \frac{b-a}{K}\right), a+k\left( \frac{b-a}{K}\right) \right] \; k=1,...K\]
Sea \(N_k\) el número de observaciones de la muestra \(X\) en el intervalo \(I_k\), es decir:
\[ N_k=\sum_{i=1}^N I_{\left( a+(k-1)\left( \frac{b-a}{K}\right), a+k\left( \frac{b-a}{K}\right) \right)} \; k=1,...,K\] Si tomamos la proporción de \(N_k\) con respecto al total de observaciones \(N\) esta es un estimador de que \(x_i\) pertenezca a \(I_k\), es decir:
\[E\left[\frac{N_k}{N}\right]=P\left( a+(k-1)\left( \frac{b-a}{K}\right)\leq x_i \leq a+k\left( \frac{b-a}{K}\right) \right)\] Como estamos con el caso de los histogramas, para cada \(I_k\), todos los puntos son iguales, un estimador para la densidad \(g(x)\) sería.
\[\hat g(x)=\frac{N_k/N}{(b-a)/K}\]
cat("histograma como función de densidad \n")histograma como función de densidad hist(df$Assault, main="Histograma Asaltos EEUU",freq=FALSE)
Estimemos el sesgo de \(\hat g(x)\).
\[E[\hat g(x)]=\frac{1}{(b-a)/K}E\left[\frac{N_k}{N}\right]=\frac{1}{(b-a)/K}\int_{a+(k-1)\left( \frac{b-a}{K}\right)}^{a+k\left( \frac{b-a}{K}\right)}g(x)dx\]
Por el teorema de valor medio de integrales, existe un \(\tilde a\) tal que:
\[g(\tilde a)=\frac{1}{(b-a)/K}\int_{a+(k-1)\left( \frac{b-a}{K}\right)}^{a+k\left( \frac{b-a}{K}\right)}g(x)dx\] Por lo que para un determinado valor c, su sesgo \(b(c)\) sería \(g(c)-g(\tilde a)\).
Observe que cuando \(K\) aumenta disminuye la longitud del intervalo y esto reduce el sesgo \(b(c)\), es decir:
\[\lim_{K\to\infty}g(\tilde a)=g(c).\]
Se podría pensar que entre más grande sea \(K\), es decir que sea más pequeño \(h=\frac{b-a}{K}\) es mejor, sin embargo, esta situación produce un sobreajuste al modelo aumentado la variabilidad del modelo.
La probabilidad que \(X_i\) pertenezca a \(I_k\) esta dada por \(\hat p=g(\tilde a)\frac{b-a}{K}\), en otras palabras, la probabilidad sigue una distribución binomial con parámetro \(\hat p\), por lo que:
\[var[\hat p]=\hat p (1-\hat p)\frac{1}{N}=g(\tilde a)\frac{b-a}{K} \left( 1-g(\tilde a)\frac{b-a}{K}\right)\frac{1}{N}\]
\[\Rightarrow var[g(\tilde a)]=var\left[ \frac{K}{b-a}\hat p\right]=\left(\frac{K}{b-a} \right)^2var(\hat p)=g(\tilde a)\left( \frac{K}{b-a}-g(\tilde a)\right)\frac{1}{N}\]
Observe que si \(K\) aumenta, la \(var[g(\tilde a)]\) aumenta, y si \(N\) aumenta la varianza disminuye.
Un problema con el histograma para estimar la densidad de una variable es que los puntos adyacentes puede ser muy distintos si están de un lado u otro de la frontera de un intervalo. En otras palabras, hay cambios discontinuos en el estimador de la densidad justo en la frontera de cada bin (partición). Esto genera mayor sesgo cerca de la frontera del intervalo que en el centro.
Si asumimos para un punto \(c\) que es el centro del intervalo \(I_k\) de longitud \(h\) como referencia obtendríamos como función de estimación para \(g(x)\) a:
\[\hat g (c) = \sum_{i=1}^N I_{\left( c-h \leq x_i \leq c+h \right)}\frac{1}{2hN}\] Esta metodología se conoce como estimador de densidad por kernel uniforme.
cat("Kernel uniforme \n")Kernel uniforme D=density(df$Assault,kernel="rectangular")
hist(df$Assault, main="DistribucionAsaltos-KernelUniforme", freq=FALSE)
lines(D,lwd=2,lty=1,col ="blue")
Este estimador de una densidad, es usual para problemas de 2 o más dimensiones, en esta metodología para el valor \(c\), se toma la observación más cercana dada una distancia específica \(d\) y posteriormente se escogen los \(K\) vecinos más cercanos según la distancia.
Sea \(d_K(x)=\underset{t}{\mathrm{argmin}} I_{(d(x_i,x)\leq t)}\geq K\) dado un t. y
\[\hat g (c)=\frac{K-1}{2d_K(x)N}\]
donde las distancias más usuales en esta metodología es la euclidiana y la de manhattan ( distancia \(L_1\)).
Leer: Zhao Puning Lai Lifeng “Analysis of KNN Density Estimation” ( Zhao & Lai (2020) ).
Rosenblatt en 1956 en su publicación “Remarks on some nonparametric estimatees of a density function” ( Rosenblatt (1956) ) propone una metodología para aproximar la densidad de una muestra más suave y con convergencia más rápida que la de los histogramas.
Una función kernel (\(K(x)\)) cumple con la siguientes características:
Además de estas características vamos a solicitar que \(K(x)\) sea una función de densidad, y que cumpla que:
Sea \(K(x)\) una función kernel, definida de la siguiente manera:
\[K(x)= \begin{cases} f(x) & x \in [-a,a] \\ 0 & \text{en otro caso} \end{cases}\]
Sea \(h \in \mathbb R^{+}\), un kernel parametrizado en \(h\) es:
\[K_h(x)=\begin{cases} \frac{1}{h}K\left( \frac{x}{h} \right) & x \in [-ah,ah] \\ 0 & \text{en otro caso} \end{cases}\]
Es claro que \(K_h(x)\) sigue cumpliendo las propiedades de un Kernel, es decir que \(K_h(x)\) es:
Para determinar que si es una función que se puede utilizar para la metodología de estimación de densidades por kernel, basta con determinar si:
Prueba
\[\int K_h(x)dx=\int_{-ah}^{ah} \frac{1}{h}K\left( \frac{x}{h} \right)dx\] tome \(u=x/h \Rightarrow du = \frac{1}{h}dx\)
\[\int K_h(x)dx=\int_{-a}^{a} K\left( u \right)du=1\]
\[E_{kh}[x]=\int xK_h(x)dx=\int_{-ah}^{ah} \frac{x}{h}K\left( \frac{x}{h} \right)dx=h\int_{-a}^{a} uK\left( u \right)du=0\]
\[var_{kh}[x]=\int x^2K_h(x)dx=\int_{-ah}^{ah} \frac{x^2}{h}K\left( \frac{x}{h} \right)dx=h^2\int_{-a}^{a} u^2K\left( u \right)du < \infty\]
La parametrización de un kernel como el indicado, sigue siendo un kernel. El paramétro \(h\) posee un rol importante para la estimación de densidades:
Si \(h>1\), la amplitud de función \(K\) se aumenta, pero el factor \(1/h\) reduce la amplitud de \(K\) para mantener el área igual a 1.
Si \(h<1\), se reduce la amplitud de \(K\), pero el factor \(1/h\) amplía el rango para que se mantenga el área igual a 1.
Podemos generalizar el concepto de kernels, cambiando la propiedad de simetría respecto a eje \(y\), a un ser simétrica con respecto a la recta \(y=x_i\) con \(x_i \in mathbb R\). El kernel parametrizado en \(h\) y centrado en \(x_i\) es
\[K_h(x)=\begin{cases} \frac{1}{h}K\left( \frac{x-x_i}{h} \right) & x \in [x_i-ah,x_i+ah] \\ 0 & \text{en otro caso} \end{cases}\]
Se puede verificar que esta traslación cumple con:
Ejercicio para el lector.
Para la muestra \(X=x_1,x_2,...,x_N\), se debe hallar un intervalo de soporte \(S=[a,b]\), tal que, \(a \leq x_i \leq b \quad \forall x_i \in X\).
Si utilizamos una función kernel \((k(x))\), de peso para cada punto de observación, con una longitud \(h\) podemos estimar una función de densidad \(f(x)\) en el punto \(c\) como:
\[\hat f(c)=\frac{1}{N}\sum_{i=1}^NK(x_i-c)=\frac{1}{hN}\sum_{i=1}^NK\left(\frac{x_i-c}{h} \right) \tag{2.9}\]
Análizaremos el sesgo de \(\hat f(c)\):
\[E[\hat f(c)]=E\left[\frac{1}{hN}\sum_{i=1}^NK\left(\frac{x_i-c}{h} \right)\right]=\frac{1}{h}E\left[K\left(\frac{x_i-c}{h} \right)\right]=\frac{1}{h}\int K\left(\frac{x-c}{h} \right)f(x)dx\] si tomamos \(y=\frac{x-c}{h}\), entonces:
\[E[\hat f(c)]=\int K(y)f(c+hy)dy\]
si \(h \to 0 \Rightarrow E[\hat f(c)]=f(c)\).
Dado que \(h\) es pequeño, es lo que buscamos, podemos realizar un desarrollo de Taylor alrededor de punto \(c\) ( asumiendo que \(f''(x)<\infty\)), es decir:
\[f(c+hy)=f(c)-hyf'(c)+\frac{1}{2}h^2y^2f''(c)+o(h^2)\]
Por lo que
\[E[\hat f(c)]=\int K(y)f(c+hy)dy=\] \[f(c)\int k(y)dy+hf'(c)\int yk(y)dy+\frac{1}{2}f''(c)\int y^2k(y)dy+o(h^2)\] \[\Rightarrow b(\hat f(c))=\frac{1}{2}h^2f''(c)var_k[y]+o(h^2)\]
Similarmente si estimamos la varianza de \(\hat f(c)\), tenemos que:
\[var[\hat f(c)]=var\left[\frac{1}{hN}\sum_{i=1}^NK\left(\frac{x_i-c}{h} \right)\right]\leq\] \[\frac{1}{Nh^2}var\left[K\left(\frac{x_i-c}{h} \right)\right]\leq\frac{1}{Nh^2}\int K\left(\frac{x-c}{h} \right)^2f(x)dx\]
Si tomamos \(y=\frac{x-c}{h}\) y realizando un desarrollo de taylor de orden 1, obtenemos que:
\[var[\hat f(c)]=\frac{1}{hN}f(c)\sigma_k+o\left(\frac{1}{hN} \right)\] donde \(\sigma_k= \int K^2(y)dy\)
la varianza disminuye cuando \(h \to 0\) y \(N \to \infty\).
nota: Recordemos que el Error cuadrático medio (MSE o ECM) se define como
\(MSE(\hat \theta)=E[(\hat \theta-\theta)^2]=var(\hat \theta)+b^2(\hat \theta)\)
Para el caso de \(\hat f(c):\)
\(MSE=\frac{1}{4}(f''(c)var_k[y])^2+o(h^4)+\frac{1}{hN}f(c)\sigma_k+o\left(\frac{1}{hN} \right)=O(h^4)+O\left(\frac{1}{hN} \right)\)
El término \(\frac{1}{4}(f''(c)var_k[y])^2+\frac{1}{hN}f(c)\sigma_k\) se le llama error cuadrático medio asintótico (AMSE)
El término \(h\) se conoce como ventana o ancho de banda o Bandwidth y como vimos anteriormente con el caso del Histograma, la selección de esta variable es importante para no caer en un sobreajuste o en un subajuste que suavice adecuadamente las curvas de densidad.
Si minimizamos el \(h\) del AMSE de \(\hat f(c)\) obtenemos:
\[h_{opt}=\left( \frac{4f(c)\sigma_k}{N(f''(c)var_k[y])^2}\right)^{1/5}=C*N^{-1/5}\] Es decir el \(MSE_{opt}= O(N^{-1/5})\).
En otras palabras, el error cuadrático medio disminuye cuando N aumenta pero a una velocidad muy despacio a un factor de 1/5.
También se puede optimizar el bandwidth, usando el MISE.
Los principales métodos empíricos para obtener el \(h\) optimo son:
Scott’s rule \(h \approx 1.06 \hat S N^{(-1/5)}\) donde \(\hat S\) es la cuasidesviación estandar.
Silverman’s rule \(h \approx 0.9 min( \hat S, \frac{IQR}{1.35}) N^{(-1/5)}\) donde \(IQR\) es el rango intercuartil ( diferencia entre el tercer y primer cuartil)
Para el caso de Silverman’s rule se asume que el núcleo (kernel) es gaussiano.
Se puede realizar metodología de validación cruzada (cross-validation) para obtener el optimo, estas metodología puede ser sesgada o insesgadas.
Improved Sheather & Jones: Es un algoritmo recursivo que minimiza el ASIME ( error cuadrático medio asintótico integrado) que depende del cuadrado de \(f''(x)\), es recomendado para datos no normales o multimodales.
Existen varios modelos de kernels utilizados para estimar la densidad de ciertos datos. En el siguiente gráfico se muestra un ejemplo de varios tipos de kernels
(kernels <- eval(formals(density.default)$kernel))[1] "gaussian" "epanechnikov" "rectangular" "triangular" "biweight"
[6] "cosine" "optcosine" plot(density(0, from = -1.2, to = 1.2, width = 2, kernel = "gaussian"),
type = "l", ylim = c(0,1), xlab = "",
main = "tipos de Kernels de la funcion Density en R", lwd=2)
for(i in 2:length(kernels))
lines(density(0, width = 2, kernel = kernels[i]), col = i, lwd=2)
legend("topright", legend = kernels, col = seq(kernels), lty = 1,lwd = 2,cex = 0.5)
Kernel Rectagular o Uniforme
Es un rectángulo que centrado sobre cada punto. Al interactuar con los kernels de los otros puntos, el efecto de la suma es un cambio abrupto.
\[K(x)=0.5, \; x \in[-1,1]\]
D=density(df$Assault,kernel="rectangular")
hist(df$Assault, main="DistribucionAsaltos-Kernel Uniforme", freq=FALSE)
lines(D,lwd=2,lty=1,col ="blue")
Kernel Triángular
Es un triángular centrado sobre cada punto. Al interactuar con los otros kernels el efecto es lineal, pero más suave que los rectángulos.
\[K(x)=1-|x|, \; x \in[-1,1]\]
D=density(df$Assault,kernel="triangular")
hist(df$Assault, main="DistribucionAsaltos-Kernel Triangular", freq=FALSE)
lines(D,lwd=2,lty=1,col ="blue")
Kernel Epanechnikov
Es uno de los kernels más estudiados. Es un segmento del arco de la parábola centrada en cada observación
\[K(x)=\frac{3}{4}(1-x^2), \; x \in[-1,1]\]
D=density(df$Assault,kernel="epanechnikov")
hist(df$Assault, main="DistribucionAsaltos-Kernel Epanechnikov", freq=FALSE)
lines(D,lwd=2,lty=1,col ="blue")
Kernel Normal o Gaussiano
Para este Kernel, se define sobre todo \(\mathbb R\) por lo que cada kernel influye en todos los otros kernels, genera una curva continua y suave
\[K(x)=\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}, \; x \in[-\infty,\infty]\]
D=density(df$Assault,kernel="gaussian")
hist(df$Assault, main="DistribucionAsaltos-Kernel Normal-Gaussiano", freq=FALSE)
lines(D,lwd=2,lty=1,col ="blue")
Kernel Biweight o Cuártico
\[K(x)=\frac{15}{16}(1-x^2)^2, \; x \in[-1,1]\]
D=density(df$Assault,kernel="biweight")
hist(df$Assault, main="DistribucionAsaltos-Kernel Biweight", freq=FALSE)
lines(D,lwd=2,lty=1,col ="blue")
Kernel Arco coseno
\[K(x)=\frac{\pi}{4}cos\left(\frac{\pi}{4}x\right), \; x \in[-1,1]\]
D=density(df$Assault,kernel="cosine")
hist(df$Assault, main="DistribucionAsaltos-Kernel Coseno", freq=FALSE)
lines(D,lwd=2,lty=1,col ="blue")
En el siguiente gráfico se puede comparar cada uno de los estimadores de la densidad de los asaltos en Estado Unidos.
D1=density(df$Assault,kernel="gaussian") #esta es por defecto
D2=density(df$Assault,kernel="cosine")
D3=density(df$Assault,kernel="biweight")
D4=density(df$Assault,kernel="epanechnikov")
D5=density(df$Assault,kernel="rectangular")
D6=density(df$Assault,kernel="triangular")
hist(df$Assault, main="DistribucionAsaltos-comparativo", freq=FALSE)
lines(D1,lwd=2,lty=1,col="blue")
lines(D2,lwd=2,lty=1,col="red")
lines(D3,lwd=2,lty=1,col="green")
lines(D4,lwd=2,lty=1,col="gray")
lines(D5,lwd=2,lty=1,col="orange")
lines(D6,lwd=2,lty=1,col="gold")
Para la construcción de funciones kernels, se debe considerar ciertas características geométricas y las propiedades de este tipo de funciones.
Kernel Cúbico Natural
En este ejemplo consideremos una función cúbica, definida entre \([-1,1]\) y 0 en otro caso, manteniendo la continuidad hasta la primera derivada con extremos libres sin curvaturas (\(K''(1)=K(-1)=0\)).
Estimaremos el polinomio cúbico \(K(x)=ax^3+bx^2+cx+d\) del tramo derecho (\(x \in [0,1]\)). Tenemos la siguientes derivadas:
Tenemos las condiciones geométricas es decir:
Por otro lado, \(K(0)=d\) la altura de la densidad del Kernel. Esto implica que \(K(x)=\frac{d}{2}x^3-\frac{3d}{2}x^2+d\).
Para que \(K(x)\) sea una función de densidad se debe cumplir que:
\(\int_{0}^{1}K(x)dx=\frac{1}{2} \Rightarrow d=\frac{4}{5}\)
Entonces:
\[K(x)=\frac{2}{5}(x^3-3x^2+2), \;x \in [0,1]\] \[K(-x)=\frac{2}{5}(-x^3-3x^2+2), \;x \in [-1,0]\] Entonces el Kernel Cúbico Natural seria
\[K(x)=\frac{2}{5}(|x|^3-3x^2+2), \;x \in [-1,1]\]
Kcubico <- function(x) {
if(abs(x)<=1){t=(2/5)*(abs(x)^3-3*x^2+2)} else {
t=0}
return(t)
}
t<-seq(-1.5,1.5,0.01)
vK<-Vectorize(Kcubico)
plot(x=t,y=vK(t),type = "l", main = "Kernel Cubico Natural")
Kernel Cúbico Plano
Supondremos en este caso que en vez de una curvatura nula, en los extremos, es que \(K''(0)=0\).
Estimaremos el polinomio cúbico \(K(x)=ax^3+bx^2+cx+d\) del tramo derecho (\(x \in [0,1]\)). Tenemos la siguientes derivadas:
Tenemos las condiciones geométricas es decir:
Por otro lado, \(K(0)=d\) la altura de la densidad del Kernel. Esto implica que \(K(x)=-dx^3+d\).
\(\int_{0}^{1}K(x)dx=\frac{1}{2} \Rightarrow d=\frac{2}{3}\)
Entonces el Kernel Cúbico Plano seria
\[K(x)=\frac{2}{3}(1-|x|^3), \;x \in [-1,1]\]
Kcubico <- function(x) {
if(abs(x)<=1){t=(2/3)*(1-abs(x)^3)} else {
t=0}
return(t)
}
t<-seq(-1.5,1.5,0.01)
vK<-Vectorize(Kcubico)
plot(x=t,y=vK(t),type = "l", main = "Kernel Cubico Plano")
Kernel Cúbico Sujeto
Este kernel asume que existe continuidad en los extremos con la primera derivada, es decir que es una conexión suave. (\(K'(1)=0\))
Ejercicio Verificar que el Kernel Sujeto es de la forma:
\[K(x)=2|x|^3-3x^2+1, \;x \in [-1,1]\]
Kcubico <- function(x) {
if(abs(x)<=1){t=2*abs(x)^3-3*x^2+1} else {
t=0}
return(t)
}
t<-seq(-1.5,1.5,0.01)
vK<-Vectorize(Kcubico)
plot(x=t,y=vK(t),type = "l", main = "Kernel Cubico Sujeto")
Ejercicios
Calcule los valores de a, b, c para que la función \(K(x)=bcos(ax)+c\),sea un kernel de densidad en el intervalo [-1,1] asumiendo que la segunda derivada en los extremos es nula (sin curvatura) y con \(b\neq 0\) y \(a \in ]0,\pi]\)
Calcule los valores de a, b, c para que la función \(K(x)=ax^2+bx+c\),sea un kernel de densidad en el intervalo [-1,1].
Calcule los valores de a, b, c para que la función \(K(x)=bcos(ax)+c\),sea un kernel de densidad en el intervalo [-1,1] asumiendo la unión con los extremos es suave y con \(b\neq 0\) y \(a \in ]0,\pi]\)
Sea \(x \in [x_1-h,x_N+h]\) una variable aleatoria con distribución \(\hat f(x)\). Demuestre que \(E[x]=\bar X\), para una m.a.s \(X=x_1,x_2,...,x_N\)
Ejercicio 1 página 325 del libro.(Capitulo 20) Wasserman (2006)
Dada una muestra aleatoria simple de la una variable aleatoria X, con los siguientes valores:
\[1 ,2.5 ,2.5 ,3,3, 3 ,4\] Usando un Kernel Triangular con bandwith de 1.5. Estime el valor de la función de densidad para los puntos 1.75, 2, 3.5