KDE – Método De Plug-In Directo
Un resumen rápido: en el Manual de optimización de KDE, derivamos la fórmula para el ancho de banda óptimo y el AMISE mínimo, de la siguiente manera:
$${{h}_{\text{opt}}}=\text{ }\sqrt[5]{\frac{R(K)}{\mu _{2}^{2}(K)R({{f}^{(2)}})n}}$$
$$\text{AMISE}({{h}_{\text{opt}}})=\frac{5}{4}\sqrt[5]{\frac{R({{f}^{(2)}})R{{(K)}^{4}}\mu _{2}^{2}(K)}{{{n}^{4}}}}$$
Donde:
- $R({{f}^{(2)}})$ se define así:
$$R({{f}^{(2)}})=\int\limits_{-\infty }^{\infty }{[{{f}^{(2)}}(x)}{{]}^{2}}dx$$ - $R(K)$ y $\mu _{2}^{2}(K)$ son cantidades constantes conocidas determinadas por la selección de la función kernel (p. ej., gaussiana).
El principal problema al usar la fórmula anterior en la práctica es que no sabemos $R({{f}^{(2)}})$: La integral de la segunda derivada al cuadrado de la función de densidad de probabilidad subyacente $f(x)$, la cual estamos tratando de estimar.
¿Cómo podemos superar este problema? Podemos estimar $R({{f}^{(2)}})$ utilizando el propio KDE.
Antes de continuar, debemos presentar un buen truco matemático: suponiendo que $\forall r>0$, ${{f}^{(r)}}(\pm \infty )\to 0$, entonces se cumple la siguiente relación:
$$\int{{{[{{f}^{(s)}}(x)]}^{2}}}dx={{(-1)}^{s}}\int{{{f}^{(2s)}}(x)f(x)dx}$$
Entonces, si vamos a aplicarlo a $R({{f}^{(2)}})$:
$$R({{f}^{(2)}})=\int{{{[{{f}^{(2)}}(x)]}^{2}}}dx=\int{{{f}^{(4)}}(x)f(x)dx}$$
Recordemos que $f(x)$ es la función de densidad de probabilidad, por lo que podemos expresar $R({{f}^{(2)}})$ como a continuación:
$$\text{R}({{f}^{(2)}})=\text{E}[{{f}^{(4)}}(x)]={{\psi }_{4}}$$
En efecto, convertimos $R({{f}^{(2)}})$ en la expectativa de la 4-ésima derivada de $f(x)$, que se puede estimar así ${{\hat{\psi }}_{4}}$ de forma no paramétrica utilizando los datos de la muestra.
$${{\hat{\psi }}_{4}}=\frac{1}{n}\sum\limits_{i=1}^{n}{{{{\hat{f}}}^{(4)}}({{x}_{i}})}$$
Donde ${{\hat{f}}^{(4)}}(.)$ es un estimador basado en datos de la cuarta derivada de $f(x)$.
${{\psi }_{r}}$ Estimación
Usando un estimador KDE con función kernel $L(.)$ y el ancho de banda $g$, el $\hat{f}(x;g)$ se define así:
$$\hat{f}(x;g)=\frac{1}{ng}\sum\limits_{i=1}^{n}{L(\frac{x-{{X}_{i}}}{g})}$$
$${{\hat{f}}^{(1)}}(x;g)=\frac{1}{n{{g}^{2}}}\sum\limits_{i=1}^{n}{{{L}^{(1)}}(\frac{x-{{X}_{i}}}{g})}$$
$$…$$
$${{\hat{f}}^{(r)}}(x;g)=\frac{1}{n{{g}^{r+1}}}\sum\limits_{i=1}^{n}{{{L}^{(r)}}(\frac{x-{{X}_{i}}}{g})}$$
Conectando ${{\hat{f}}^{(r)}}(x;g)$, el estimador de la derivada r-ésima esperada ${{\hat{\psi }}_{r}}$ se expresa así:
$${{\hat{\psi }}_{r}}=\frac{1}{{{n}^{2}}{{g}^{r+1}}}\sum\limits_{i=1}^{n}{\sum\limits_{j=1}^{n}{{{L}^{(r)}}(\frac{{{X}_{i}}-{{X}_{j}}}{g})}}$$
La siguiente pregunta es qué es $g$? Es la función de kernel $L(.)$ (y el ancho de banda $g$). ¿Es el mismo que usamos para el KDE final? No necesariamente, ya que nuestro objetivo es minimizar el error en el ${{\hat{\psi }}_{r}}$ estimado.
Bajo ciertos supuestos de regularidad (Wand y Jones (1995)), el sesgo asintótico y la varianza de ${{\hat{\psi }}_{r}}$ se obtienen, por lo que podemos calcular el error cuadrático medio asintótico (AMSE):
$$\text{AMSE}[{{\hat{\psi }}_{r}}(g)]=\left[ \frac{{{L}^{(r)}}(0)}{n{{g}^{r+1}}}+\frac{{{\mu }_{2}}(L){{\psi }_{r+2}}\times {{g}^{2}}}{4} \right]+\frac{2R({{L}^{(r)}}){{\psi }_{o}}}{{{n}^{2}}{{g}^{2r+1}}}+\frac{4}{n}\left[ \int{[{{f}^{(r)}}}(x){{]}^{2}}f(x)dx-\psi _{r}^{2} \right]$$
Y el ancho de banda óptimo de AMSE es:
$${{g}_{\text{AMSE}}}=\sqrt[r+k+1]{-\frac{k!{{L}^{(r)}}(0)}{{{\mu }_{k}}(L){{{\hat{\psi }}}_{r+k}}\times n}}$$
Donde:
- $k$ es el número de etapas, siempre que ${{\mu }_{k}}(L)>0$.
- Para $k=0$ (Método Silverman), comenzamos con una distribución conocida (por ejemplo, Gaussiana), calculamos ${{\hat{\psi }}_{4}}$ (analíticamente), y calculamos el ancho de banda óptimo.
- Para $k=2$ (Método de Plug-in directo), empezamos con distribución conocida (ej., Gaussiana). Entonces:
- Calcule ${{\hat{\psi }}_{8}}$ (analíticamente).
- Etapa 1:
- Calcule el valor de ancho de banda óptimo para estimar ${{\hat{\psi }}_{6}}$
- Calcule ${{\hat{\psi }}_{6}}$.
- Etapa 2:
- Calcule el valor de ancho de banda óptimo para estimar ${{\hat{\psi }}_{8}}$
- Calcule ${{\hat{\psi }}_{4}}$.
- Calcule el ancho de banda óptimo de KDE.
- Para $k=4$ (método Direct Plug-in), comenzamos con una distribución conocida (por ejemplo, Gaussiana). Después:
- Calcule ${{\hat{\psi }}_{12}}$ (analíticamente).
- Etapa 1:
- Calcule el valor de ancho de banda óptimo para estimar ${{\hat{\psi }}_{10}}$.
- Calcule $ {{\hat{\psi }}_{10}}$.
- Etapa 2:
- Calcule el valor de ancho de banda óptimo para estimar ${{\hat{\psi }}_{8}}$.
- Calcule ${{\hat{\psi }}_{8}}$.
- Etapa 3:
- Calcule el valor de ancho de banda óptimo para estimar ${{\hat{\psi }}_{6}}$.
- Calcule ${{\hat{\psi }}_{6}}$.
- Etapa 4:
- Calcule el valor de ancho de banda óptimo para estimar ${{\hat{\psi }}_{4}}$.
- Calcule ${{\hat{\psi }}_{4}}$.
- Calcule el ancho de banda óptimo de KDE.
Por lo general, dos etapas ($k=2$) se consideran una buena compensación entre el sesgo (mitigados a medida que $k$ aumenta), y varianza (aumenta con $k$).
Método de Plug-in Directo (Sheather & Jones)
Este es el método propuesto por Sheather y Jones (1991), donde consideran $L=K$ y $k=2$, produciendo lo que llamamos Direct Plug-In (DPI). El algoritmo es:
- Usando los datos de la muestra, calcule $\hat{\sigma }=\min (s,{{\hat{\sigma }}_{\text{IQR}}})$.
- Suponga una distribución subyacente gaussiana (ej., $f(x)=\phi (x)=\frac{1}{\sqrt{2\pi }\sigma }\exp \left( -\frac{{{(x-\mu )}^{2}}}{2{{\sigma }^{2}}} \right)$), luego calcule (analíticamente) ${{\hat{\psi }}_{8}}$.
$${{\hat{\psi }}_{8}}=\int\limits_{-\infty }^{\infty }{{{\phi }^{(8)}}(x)\phi (x)dx}=\frac{105}{32\sqrt{\pi }{{{\hat{\sigma }}}^{9}}}$$ - Calcule el ancho de banda óptimo para ${{\hat{\psi }}_{6}}$, ${{g}_{1}}$:
$${{g}_{1}}=\sqrt[9]{-\frac{2{{K}^{(6)}}(0)}{{{\mu }_{2}}(K){{{\hat{\psi }}}_{8}}\times n}}$$ - Estime ${{\hat{\psi }}_{6}}$:
$${{\hat{\psi }}_{6}}=\frac{1}{{{n}^{2}}g_{1}^{7}}\sum\limits_{i=1}^{n}{\sum\limits_{j=1}^{n}{{{K}^{(6)}}(\frac{{{X}_{i}}-{{X}_{j}}}{{{g}_{1}}})}}$$ - Calcule el ancho de banda óptimo para ${{\hat{\psi }}_{4}}$, ${{g}_{2}}$:
$${{g}_{2}}=\sqrt[7]{-\frac{2{{K}^{(4)}}(0)}{{{\mu }_{2}}(K){{{\hat{\psi }}}_{6}}\times n}}$$ - Estime ${{\hat{\psi }}_{4}}$:
$${{\hat{\psi }}_{4}}=\frac{1}{{{n}^{2}}g_{2}^{5}}\sum\limits_{i=1}^{n}{\sum\limits_{j=1}^{n}{{{K}^{(4)}}(\frac{{{X}_{i}}-{{X}_{j}}}{{{g}_{2}}})}}$$ - Ahora, usando ${{\hat{\psi }}_{4}}$ como un estimado para $R({{f}^{(2)}})$:
$${{h}_{\text{DPI}}}=\text{ }\sqrt[5]{\frac{R(K)}{\mu _{2}^{2}(K){{{\hat{\psi }}}_{4}}\times n}}$$
Conclusión
En este artículo, asumimos que la función kernel $k$ no solo es simétrica sino que también tiene cuatro (4) derivadas continuas. La suposición excluye el uso de muchas funciones kernel (p. ej., uniforme, triangular, bipeso, tripeso) pero, afortunadamente, la Gaussiana cumple las condiciones y, a menudo, se usa con el método DPI.
El método de complemento directo de Sheather y Jones es popular en la práctica para un amplio conjunto de casos, y ofrece un buen rendimiento para densidades uniformes, al menos en simulación.
Referencias
- Silverman, B.W. (1986). Density Estimation for Statistics and Data Analysis. Chapman & Hall/CRC London.
- W. Zucchini, Applied smoothing techniques, Part 1 Kernel Density Estimation., 2003.
- Byeong U. Park and J. S. Marron. Comparison of Data-Driven Bandwidth Selectors. Journal of the American Statistical Association Vol. 85, No. 409 (Mar., 1990), pp. 66-72 (7 pages).
- S.J. Sheather and M.C. Jones. A reliable data-based bandwidth selection method for kernel density estimation. J. Royal Statist. Soc. B, 53:683-690, 1991.
