Suponer que una distribución pertenece a una familia paramétrica tiene el inconveniente de que se impone de partida un conjunto de propiedades. Por ejemplo, si suponemos que los datos son normales estamos imponiendo que la distribución es simétrica y unimodal, y que aproximadamente un 5 % de los datos dista de la media más de dos desviaciones típicas. Por el contrario, las técnicas estadísticas no paramétricas reducen al máximo las propiedades que se asumen sobre la distribución de los datos.
El objetivo de este capítulo es estimar una función de densidad \(f\) sin hacer hipótesis previas sobre ella a partir de una muestra de v.a.i.i.d. \(X_1,\ldots,X_n\) de una población con función de densidad \(f\). Tal vez los estimadores más famosos para este propósito son los estimadores de núcleo que se trataron en el primer capítulo como herramienta exploratoria. Recordamos de nuevo la definición: \[\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 asumimos que es una función de densidad, \(K\geq 0\) y \(\int K = 1\)).
Se recomienda recordar lo que se explicó en la sección 1.5.4 antes de leer este capítulo. Aquí profundizamos en las propiedades estadísticas de estos estimadores. Además del interés que puedan tener por sí mismos, proporcionan un ejemplo relativamente simple de algunas nociones importantes y recurrentes en estadística y aprendizaje automático:
Sabemos que un estimador de núcleo es una función de densidad (sección 1.5.4) y por lo tanto define una distribución de probabilidad. En esta sección estudiamos con más detalle cuál es esta distribución: intuitivamente, es una versión suavizada (una convolución) de la distribución de probabilidad discreta correspondiente a la función de distribución empírica (sección 2.5).
Dadas las observaciones \(X_1,\ldots,X_n\) consideramos la correspondiente función de distribución empírica, \(F_n\), que es la que da probabilidad \(1/n\) a cada observación \(X_i\), y definimos la variable \(U = Y + Z\), donde
Veamos cuál es la función de distribución de \(U\), que denotamos por \(\hat{F}\) (nótese que depende de los datos a través de la distribución empírica). Por la fórmula de la probabilidad total, \[\hat{F}(x) = \mbox{P}(Y+Z\leq x) = \frac{1}{n} \sum_{i=1}^n \mbox{P}(Z\leq x-Y|Y=X_i) = \frac{1}{n} \sum_{i=1}^n \mbox{P}(Z\leq x-X_i).\] En la última igualdad se usa la independencia entre \(Y\) y \(Z\). Usando que la densidad de \(Z\) es \(K_h\), tenemos \[\frac{1}{n} \sum_{i=1}^n \mbox{P}(Z\leq x-X_i) = \frac{1}{n} \sum_{i=1}^n \int_{-\infty}^{x-X_i} K_h(z) dz = \frac{1}{nh} \sum_{i=1}^n \int_{-\infty}^{x-X_i} K(z/h) dz.\] Ahora, si \(K\) es continua, por el teorema fundamental del cálculo se cumple: \[\hat{F}'(x) = \frac{1}{nh} \sum_{i=1}^nK\left(\frac{x-X_i}{h}\right) = \hat{f}(x).\]
En resumen, la receta para generar realizaciones de una variable con densidad \(\hat{f}\) es:
rnorm() de R si el núcleo es gaussiano.Consideramos un ejemplo con los datos del fichero faithful de R. Se trata de la duración de las erupciones y del tiempo hasta la siguiente erupción del geyser Old Faithful del parque de Yellowstone. Se representa la nube de puntos de las dos variables y el estimador que hemos descrito (fijando el valor \(h = 0.15\) y el núcleo gaussiano, que es la opción por defecto):
graf1 <- ggplot(faithful) +
geom_point(aes(x = eruptions, y = waiting))
graf2 <- ggplot(faithful) +
geom_density(aes(x = eruptions),
bw = 0.15,
fill = 'olivedrab4',
size = 1.2)
graf1 / graf2 # requiere el paquete patchwork
kernel = c("gaussian", "epanechnikov", "rectangular", "triangular", "biweight", "cosine", "optcosine"). Úsalos para ver cómo cambia el aspecto del estimador al cambiar el núcleo.rnucleo <- function(muestra, h){
# genera n observaciones de una distribución correspondiente a un estimador
# del núcleo (calculado con 'muestra') con núcleo gaussiano y parámetro de suavizado h
n <- length(muestra)
y = sample(muestra, n, rep = TRUE) + rnorm(n, sd = h)
return(y)
}Los ejercicios anteriores muestran que la elección del parámetro de suavizado es crucial para obtener un buen estimador, mientras que la elección de núcleo no es tan importante. En la sección siguiente vamos a estudiar qué condiciones debe cumplir \(h\) para que el estimador sea consistente.
Como criterios para evaluar la calidad de un estimador, se suelen utilizar el sesgo y la varianza. Si lo que queremos estimar es \(f(x)\), para un valor fijo de \(x\):
Para que \(\hat f(x)\) sea (débilmente) consistente basta demostrar que tanto el sesgo como la varianza de \(\hat{f}(x)\) convergen a cero. La interpretación como convolución de la sección anterior es la causa de que haya un sesgo en la estimación. Cuanto mayor es \(h\) más se separa \(F_n \ast K_h\) de \(F_n\) y podemos esperar que el sesgo sea una función creciente de \(h\). Intuitivamente, si queremos que el sesgo converja a cero para tamaños muestrales grandes debemos seleccionar \(h=h_n\to 0\) si \(n\to \infty\). Además, para que la varianza tienda a cero necesitamos que \(nh_n\) vaya a infinito. Es decir, el papel que hace \(n\) en los estimadores paramétricos habituales, aquí lo hace \(nh_n\) (lo que refleja que al suavizar perdemos una parte de la información muestral). El siguiente resultado formaliza estas intuiciones:
Prueba. Basta demostrar que tanto el sesgo como la varianza de \(\hat{f}(x)\) convergen a cero. Demostramos en primer lugar que el sesgo tiende a cero, esto es, \(\lim_{n\to\infty}\mbox{E}[\hat{f}(x)] = f(x).\)
Haciendo un cambio de variable muy sencillo tenemos que \[\mbox{E}[\hat{f}(x)] = \mbox{E}\left[\frac{1}{h}K\left(\frac{x-X}{h}\right)\right] = \frac{1}{h}\int K\left(\frac{x-u}{h}\right)f(u) du = \int K(w)f(x-hw)dw. \tag{7.1}\] Dado que \(f\) es continua en \(x\) y \(h\to 0\), si pudiéramos intercambiar el límite con la integral tendríamos \[\lim_{n\to\infty}\mbox{E}[\hat{f}(x)] = \int K(w) f(x) dw = f(x).\] El intercambio se puede hacer por el teorema de la convergencia dominada, ya que \(K(w)f(x-hw)\leq M K(w)\), que es una función integrable.
A continuación, vamos a demostrar que \(\lim_{n\to\infty} nh\mbox{Var}[\hat{f}(x)] <\infty\), y por lo tanto \(\lim_{n\to\infty} \mbox{Var}[\hat{f}(x)] = 0\) (puesto que estamos suponiendo que \(nh\to\infty\)).
Obsérvese que \[nh\mbox{Var}[\hat{f}(x)] = h\left(\mbox{E}\left[\frac{1}{h^2}K^2\left(\frac{x-X}{h}\right)\right]- \left(\mbox{E}\left[\frac{1}{h}K\left(\frac{x-X}{h}\right)\right]\right)^2 \right).\]
Basta demostrar que el primer término converge a un límite finito. Para el segundo término podemos aplicar la primera parte de la demostración.
Procedemos de forma análoga,
\[\mbox{E}\left[\frac{1}{h}K^2\left(\frac{x-X}{h}\right)\right] = \int K^2(w)f(x-hw)dw \to f(x)\int K^2(w)dw<\infty.\] El teorema de la convergencia dominada se aplica para justificar el intercambio de la integral con el límite. Si combinamos las dos partes de la demostración se tiene \[nh\mbox{Var}[\hat{f}(x)]\to f(x)\int K^2(w)dw<\infty,\ \ \ \ \ \mbox{si}\ h\to 0,\ nh\to\infty. \tag{7.2}\]
Junto con la media y la mediana, la moda es otra medida habitual de tendencia central. En estadística elemental se suele definir la moda como el valor que más se repite. Sin embargo, esta definición es insatisfactoria para distribuciones continuas para las que no hay apenas repeticiones en los datos. Algún tipo de suavizado es necesario para definir la moda muestral de variables continuas.
Supongamos que para la densidad \(f\) hay un único valor \(\theta\) tal que \(f(\theta) = \max_{x\in\mathbb{R}} f(x)\). El parámetro \(\theta\) es la moda poblacional.
Dado un estimador de núcleo \(\hat{f}\) se define la moda muestral como el valor \(\hat{\theta}\) tal que \(\hat f(\hat \theta) = \max_{x\in\mathbb{R}}\hat f(x)\).
Bajo ciertas condiciones de regularidad este estimador es consistente:
Teorema 7.2 Sea \(\hat f\) el estimador de núcleo obtenido a partir de una muestra de v.a.i.i.d. \(X_1,\ldots,X_n\) de una distribución con función de densidad \(f\). Se supone que \(f\) es uniformemente continua y la moda poblacional es única. Además, se supone que el parámetro de suavizado verifica \(h_n\to 0\) y \(nh_n^2\to \infty\). Entonces, la moda muestral \(\hat\theta\) converge a la poblacional \(\theta\) en probabilidad.
La demostración de este resultado se puede encontrar en el artículo original de Parzen (1962).
El resultado de consistencia que hemos dado anteriormente se refiere a lo que ocurre localmente en un punto \(x\) prefijado, pero nosotros estamos interesados en el comportamiento global de \(\hat f\) como estimador de la función de densidad.
Un criterio local que tiene en cuenta tanto el sesgo como la varianza es el error cuadrático medio \(\mbox{E}[(\hat{f}(x) - f(x))^2]\). Sabemos que: \[\mbox{ECM}(x) = \mbox{Sesgo}^2[\hat{f}(x)] + \mbox{Var}[\hat{f}(x)].\] Como queremos estimar la función de densidad \(f\) globalmente, no sólo para un valor fijo de \(x\), consideramos el error cuadrático medio integrado: \[\mbox{ECMI}(\hat{f}) = \int \mbox{ECM}(x) dx = \int \mbox{Sesgo}^2(x) dx + \int \mbox{Var}[\hat{f}(x)]dx.\] Como vamos a ver que el ECMI depende sobre todo del parámetro de suavizado \(h\), del tamaño muestral \(n\) y de algunas características de suavidad de la función que se desea estimar.
Si aplicamos el teorema de Fubini, \[\mbox{ECMI}(\hat{f}) = \int \mbox{E}[(\hat{f}(x) - f(x))^2] dx=\mbox{E}\left[\int(\hat f(x)-f(x))^2 dx\right]=\mbox{E}\left[\|\hat f-f\|_2^2\right],\] donde \(\|f\|_2 := \big(\int f(x)^2 dx\big)^{1/2}\). El ECMI es el valor esperado de la distancia \(L_2\) entre \(f\) y \(\hat{f}\). Dado que no todas las funciones de densidad verifican \(\int f(x)^2 dx <\infty\), algunos autores argumentan que para valorar estimadores de la densidad es más apropiado usar criterios \(L_1\) de la forma \[\mbox{E}\left[ \int |\hat{f}(x) - f(x)| dx \right],\] en los que se reemplaza la diferencia al cuadrado por el valor absoluto. En análisis funcional se usa la notación \(\|f\|_1 = \int|f(x)|dx\). Con esta notación, el criterio \(L_1\) se escribe \(\mbox{E}(\|\hat{f}-f\|_1)\). La teoría para el criterio \(L_1\) es más complicada debido a la no derivabilidad del valor absoluto, así que aquí nos vamos a centrar en el ECMI.
Es muy importante seleccionar \(h=h_n\) de forma que se equilibren el sesgo y la varianza. La figura 7.1 ilustra dos malas elecciones de \(h\) que llevan a infrasuavizado y sobresuavizado, respectivamente, de manera que el sesgo y la varianza están desequilibrados, junto con otra elección más adecuada. Para cada elección de \(h\) se ha calculado el estimador para un cierto número de muestras independientes. También se ha representado la densidad de la que proceden los datos (normal). En el gráfico superior se ha seleccionado un valor de \(h\) demasiado grande lo que lleva a una varianza pequeña pero a un estimador muy sesgado. Para las muestras generadas los estimadores siempre salen parecidos y siempre infraestiman en la parte central y sobreestiman en las colas. En el gráfico intermedio se ha elegido \(h\) demasiado pequeña con lo que el estimador resulta muy variable aunque el sesgo es pequeño. Los estimadores sales muy distintos para cada muestra pero su promedio más o menos coincide con la función que queremos estimar. En el gráfico inferior se ha aplicado un método automático de selección del parámetro de suavizado (debido a Sheather y Jones) que da un resultado más equilibrado.
En la figura 7.2 se representa cómo depende el ECMI del parámetro de suavizado. El objetivo es determinar aproximadamente el valor \(h^*\) para el que se minimiza el ECMI.
Como el problema es demasiado complicado, en lugar de optimizar el ECMI directamente vamos a minimizar una aproximación. En lo que resta de esta sección se deducen expresiones aproximadas de los términos de sesgo y varianza del ECMI de los estimadores de núcleo, de manera que nos sirvan de guía para la adecuada selección de \(h\).
Vamos a suponer que el núcleo es una función simétrica, con \(\int K(u) du=1\), \(\int uK(u)du = 0\), \(\sigma^2_K = \int u^2K(u)du<\infty\), y \(\|K\|_2^2 = \int K(u)^2 du <\infty\). Por ejemplo, un núcleo gaussiano cumple todas estas propiedades. También suponemos que \(f\) es derivable dos veces con derivada continua. Entonces, el sesgo de \(\hat{f}(x)\) se puede aproximar a partir de la ecuación (7.1). Para ello, consideramos los dos primeros términos del desarrollo de Taylor de \(f\). Si \(h\) es pequeño, \[\mbox{E}[\hat{f}(x)] \approx \int K(w)\left[f(x) -whf'(x) + \frac{w^2h^2}{2}f''(x)\right]dw=f(x) +\frac{h^2}{2}f''(x)\sigma^2_K.\] En la última iguadad estamos usando las hipótesis sobre \(K\).
Como consecuencia de la ecuación anterior, \[\mbox{Sesgo}^2[\hat{f}(x)] \approx \frac{h^4}{4} f''(x)^2 \sigma^4_K,\ \ \mbox{si}\ \ h\ \mbox{pequeño}.\] Respecto a la varianza tenemos, a partir de la ecuación (7.2), \[\mbox{Var}[\hat{f}(x)] \approx \frac{1}{nh} f(x) \|K\|_2^2,\ \ \mbox{si}\ \ h\ \mbox{pequeño},\ nh\ \ \mbox{grande}.\]
Una posible estrategia para seleccionar \(h\) es elegir el valor para el que se minimiza la aproximación al ECMI. Si derivamos respecto de \(h\) en la ecuación (7.3) e igualamos a cero, el valor \(h\) óptimo resultante es \[h^* = \left(\frac{\|K\|_2^2}{\sigma_K^4\|f''\|_2^2}\right)^{1/5} n^{-1/5}, \tag{7.4}\] y sustituyendo este valor en la expresión aproximada del ECMI tenemos \[\mbox{ECMI}^* \approx \frac{5}{4} \sigma_K^{4/5} \|K\|_2^{8/5}\|f''\|_2^{2/5} n^{-4/5}. \tag{7.5}\]
Nos planteamos en esta sección encontrar el núcleo óptimo, es decir, la función \(K\) para la que se minimiza la expresión de \(\mbox{ECMI}^*\) dada por la ecuación (7.5).
La parte que depende del núcleo es \(\phi(K) = \sigma_K^{4/5} \|K\|_2^{8/5}\). Esta función es invariante ante cambios de escala del núcleo (probar como ejercicio que \(\phi(K_h)=\phi(K)\), para \(h>0\), donde \(K_h(\cdot)=K(\cdot/h)/h\)). Por lo tanto, para encontrar el núcleo óptimo se puede suponer sin pérdida de generalidad que \(\sigma_K=1\) y resolver el siguiente problema de cálculo de variaciones: \[\min \|K\|_2^2 \ \ \mbox{s.a.}\ \ K\geq 0,\ \int K(u) du=1,\ \int uK(u)du = 0,\ \sigma^2_K = 1.\] Resulta que la solución de este problema viene dada por el llamado núcleo de Epanechnikov: \[K^*(t) = \frac{3}{4\sqrt{5}}\left(1-\frac{t^2}{5}\right),\ \ \ \mbox{si}\ \ -\sqrt{5}\leq t\leq \sqrt{5}.\] Este resultado, aunque interesante teóricamente, no tiene mucha importancia práctica puesto que en realidad la elección del núcleo no afecta demasiado al ECMI. Si calculamos el cociente \(\phi(K^*)/\phi(K)\) para distintos núcleos, los valores suelen estar por encima de \(0.9\). Por ejemplo, para el núcleo rectangular el cociente es 0.9295 y para el núcleo gaussiano, 0.9512. Esto significa que no hay mucha diferencia entre usar el núcleo óptimo y, por ejemplo, el núcleo gaussiano.
Se han propuesto multitud de métodos para seleccionar \(h\) en la práctica. En esta sección hacemos un resumen de las principales ideas involucradas sin entrar en demasiados detalles técnicos. Básicamente hay dos tipos de procedimientos: los métodos plug-in que resultan de estimar la parte desconocida en la expresión de \(h^*\) (esto es, \(R(f'') \equiv \|f''\|_2^2\), que esencialmente es una medida de curvatura de \(f\)) por una estimación o aproximación basada en alguna hipótesis más o menos adecuada, y los métodos basados en validación cruzada.
El lector interesado en el problema de selección del parámetro de suavizado puede consultar la revisión de Jones, Marron y Sheather (1996).
En este apartado hay distintas estrategias en función de las hipótesis que estemos dispuestos a asumir sobre \(f\):
Si \(f\) es la densidad de una v.a. normal con media \(\mu\) y varianza \(\sigma^2\), entonces se puede demostrar que \[R(f'') = \frac{3}{8\sqrt{\pi}\sigma^5}.\] Reemplazando este valor en la expresión de \(h^*\) y suponiendo que \(K\) corresponde a una densidad normal estándar resulta: \[h^* = \sigma(4/(3n))^{1/5}\approx \sigma (1.0456) n^{-1/5}.\] De aquí, un posible valor para el parámetro de suavizado es \(\hat{\sigma} (1.0456) n^{-1/5}\), donde \(\hat{\sigma}\) es un estimador de \(\sigma\) obtenido a partir de la muestra. La opción propuesta por Silverman es \(\hat{\sigma} = \min\{s,\hat{\sigma}_{ri}\}\), donde \(\hat{\sigma}_{ri}\) es el rango intercuartílico (estandarizado para que converja a \(\sigma\) bajo normalidad): \[\hat{\sigma}_{ri} = \frac{X_{([0.75n])} - X_{([0.25n])}}{\Phi(0.75)-\Phi(0.25)}.\] La función de R que realiza estos cálculos es bw.nrd (una ligera variación es bw.nrd0):
set.seed(100)
n <- 100
x <- rnorm(n)
bw.nrd(x = x)
#> [1] 0.3982919
# El mismo resultado (salvo redondeo) que
ri <- diff(quantile(x, c(0.25, 0.75))) / diff(qnorm(c(0.25, 0.75)))
1.0456 * n^(-1/5) * min(sd(x), ri)
#> [1] 0.390266
# Una ligera variación (se multiplica por 0.9 en lugar de 1.0456)
bw.nrd0(x = x)
#> [1] 0.3381724Otra posibilidad es fijar un parámetro de suavizado piloto \(g\) (usando por ejemplo la regla descrita en el apartado anterior) para obtener un estimador de núcleo preliminar \(\hat{f}_g(x)\) y luego utilizar \(\hat{R}(f'') = \int \hat{f}''_g(x)^2 dx\).
El método de Sheather y Jones es un refinamiento de esta idea. Está implementado en R: bw.SJ. Es uno de los métodos más recomendados:
bw.SJ(x = x)
#> [1] 0.3663063En este caso concreto (en el que los datos simulados realmente son normales) no hay mucha diferencia con los métodos anteriores.
Los métodos de validación cruzada (cross validation) son muy socorridos cuando se trata de fijar el valor de una constante de ajuste. Son métodos en los que la muestra se divide en dos partes y se usa una de ellas para obtener información del procedimiento calculado con la otra. Este esquema se suele repetir muchas veces y se promedian los resultados obtenidos. Veamos cómo se aplica un procedimiento de este tipo en el caso del parámetro de suavizado de un estimador de núcleo.
Recordemos que \[\mbox{ECMI}(\hat{f})= \mbox{E}\left[\int (\hat{f}(x; h) - f(x) )^2 dx\right] = \mbox{E}\left[\int \hat{f}(x; h)^2 dx\right] - 2\mbox{E}\left[\int \hat{f}(x ; h)f(x)dx\right] + \int f(x)^2dx.\] En la ecuación anterior se ha hecho explícita la dependencia de \(h\) para mayor claridad. El último término no depende de \(h\), por lo que el objetivo es minimizar la diferencia de los dos primeros. Este valor se puede estimar mediante \[C(h) = \int \hat{f}(x; h)^2 dx - \frac{2}{n}\sum_{i=1}^n \hat{f}_{(-i)}(X_i; h),\] donde \(\hat{f}_{(-i)}\) denota el estimador de núcleo calculado con todas las observaciones a excepción de \(X_i\).
Para explicar la estimación del segundo término observamos lo siguiente: si tuviéramos una nueva observación \(X\), independiente de las observaciones \(X_1,\ldots,X_n\) usadas para calcular \(\hat{f}\) entonces \[\int \hat{f}(x; h)f(x)dx = \mbox{E}[\hat{f}(X; h)| X_1,\ldots,X_n].\] En la realidad no disponemos de esa observación adicional pero la igualdad anterior sugiere estimar \(\int \hat{f}(x; h)f(x)dx\) mediante \(n^{-1}\sum_{i=1}^n \hat{f}_{(-i)}(X_i; h).\)
Finalmente, para determinar el parámetro de suavizado se resuelve numéricamente el siguiente problema de minimización: \[\hat{h} = \mbox{arg}\min_{h>0} C(h) = \mbox{arg}\min_{h>0} \left[\int \hat{f}(x; h)^2 dx - \frac{2}{n}\sum_{i=1}^n \hat{f}_{(-i)}(X_i; h)\right].\] Desgraciadamente, el problema anterior es difícil ya que la función objetivo puede tener muchos mínimos locales, y el resultado puede depender del método de resolución utilizado.
Este método está implementado en bw.ucv:
bw.ucv(x = x)
#> [1] 0.4224203En el capítulo sobre descripción de datos ya vimos una manera de representar gráficamente los estimadores de núcleo usando ggplot2. En este apartado profundizamos un poco más en la implementación de estos estimadores en R.
El comando que se utiliza para calcular los estimadores de núcleo es density. El primer argumento es la muestra. Del resto de argumentos, los más importantes son bw (el valor de \(h\), por defecto bw.nrd0 explicado más arriba) y kernel (\(K\), por defecto gaussiano). El resultado es una lista cuyos elementos más importantes son x (las coordenadas de los puntos en los que se calcula \(\hat{f}\), una malla en un intervalo que depende del rango muestral) e y (los valores \(\hat{f}(x)\)).
A continuación calculamos (y representamos con comandos de R base) los estimadores de núcleo correpondientes a tres métodos diferentes de selección del parámetro de suavizado:
# Generamos n datos de distribución chi2 con 3 gl
n <- 100
x <- rchisq(n, 3)
# Opciones por defecto
nucleo <- density(x)
plot(nucleo, ylim = c(0,0.25))
curve(dchisq(x,3), from=0, to=15, lty=2, add=TRUE)
# Parámetro de suavizado de Sheather-Jones
nucleo_SJ <- density(x, bw = "SJ")
plot(nucleo_SJ, ylim = c(0,0.25))
curve(dchisq(x,3), from=0, to=15, lty=2, add=TRUE)
# Parámetro de suavizado por validación cruzada
nucleo_vc <- density(x, bw = "ucv")
plot(nucleo_vc, ylim = c(0,0.25))
curve(dchisq(x,3), from=0, to=15, lty=2, add=TRUE)
Si los datos son vectores de dimensión \(d\), el estimador de núcleo multivariante se define: \[\hat{f}(x) = \frac{1}{n|H|^{1/2}}\sum_{i=1}^n \tilde{K}(H^{-1/2}(x-X_i)),\ \ \ x=(x_1,\ldots,x_d)\in\mathbb{R}^d.\] donde \(\tilde{K}\) es ahora un núcleo (una función de densidad) multivariante, \(H=\Sigma^{1/2}\) para una matriz \(\Sigma\) definida positiva \(d\times d\), y \(|H|\) es el determinante de \(H\).
Con la generalidad con la que hemos definido \(\hat{f}\), el usuario tendría que elegir toda una matriz de constantes de ajuste. Por ello, las siguientes simplificaciones son habituales:
Con las simplificaciones anteriores, el estimador se reduce a: \[\hat{f}(x) = \frac{1}{nh^d}\sum_{i=1}^n \prod_{j=1}^d K\left(\frac{x_j - X_{i,j}}{h}\right),\ \ \ x=(x_1,\ldots,x_d)\in\mathbb{R}^d.\]
Para visualizar las curvas de nivel de un estimador de núcleo bidimensional se puede usar la función geom_density_2d en ggplot2. Internamente, este comando usa la función MASS::kde2d para calcular el estimador. Por defecto, se utiliza un producto de núcleos gaussianos con el mismo parámetro de suavizado en cada coordenada, seleccionado mediante bw.nrd. Los detalles de uso se pueden consultar aquí y aquí.
Veamos un ejemplo con los datos del geyser Old faithful:
ggplot(faithful, aes(x = eruptions, y = waiting)) +
geom_point() +
xlim(0.5, 6) + ylim(40, 110) +
geom_density_2d()
A no ser que se disponga de tamaños muestrales enormes, en dimensiones altas es difícil encontrar datos en muchas zonas del espacio muestral. Esto hace que se deterioren las propiedades de los estimadores. Por ejemplo, bajo condiciones de regularidad poco exigentes, puede probarse que para datos de dimensión \(d\) la aproximación asintótica al ECMI óptimo verifica: \[\mbox{ECMI}^* \approx O(n^{-4/(4+d)}).\] El caso \(d=1\) corresponde a la ecuación (7.3) Este resultado sugiere que la convergencia a cero del ECMI se hace más lenta a medida que la dimension de los datos crece.
En el bootstrap no paramétrico tradicional las remuestras se extraen de la función de distribución empírica. En este capítulo hemos estudiado otra manera de estimar la distribución (vía su función de densidad) no paramétricamente. Esto abre la posibilidad de obtener las remuestras a partir de la distribución definida por el estimador de núcleo. A esta versión del bootstrap se le llama bootstrap suavizado. En general, introduce un sesgo a cambio de una reducción de varianza, lo que en algunos problemas puede dar un mejor resultado.
Veamos un ejemplo. En la siguiente simulación se comparan los resultados del bootstrap y del bootstrap suavizado para estimar la desviación típica de la mediana. Se generan 100 muestras de tamaño 21 de una distribución de Cauchy. Recordemos que, en este caso, la desviación típica de la mediana se conoce y se puede aproximar por 0.37 (DasGupta (2008) pag. 473). Se ha señalado este valor en la línea horizontal del gráfico. Para cada muestra se aplica el método bootstrap habitual y el bootstrap suavizado (\(h=0.5\)) para calcular el error típico de la mediana y se representan los resultados mediante diagramas de cajas.
# Funciones ---------------------------------------------------------------
sd_median_bootstrap <- function(x, R = 200){
# Estima el error típico de la mediana mediante bootstrap
# y bootstrap suavizado
# R = number of resamples
n <- length(x)
# Bootstrap
resamples <- matrix(sample(x, n*R, rep = TRUE), ncol = R)
boots_median <- apply(resamples, 2, median)
# Bootstrap suavizado
bandwidth <- 0.5
resamples <- resamples + rnorm(n*R, sd = bandwidth)
boots_sm_median <- apply(resamples, 2, median)
return(c(sd(boots_median), sd(boots_sm_median)))
}
# Simulación --------------------------------------------------------------
set.seed(100)
# Parámetros
R <- 100
n <- 21
gl <- 1
boot_results <- NULL
smoothed_boot_results <- NULL
for (i in 1:R){
sample <- rt(n, gl)
boot_results <- rbind(boot_results, sd_median_bootstrap(sample))
}
metodo <- gl(2, R, labels = c('Boostrap', 'Bootstrap suavizado'))
df <- data.frame(sd_median = c(boot_results[,1], boot_results[,2]),
metodo = metodo)
ggplot(df) +
geom_boxplot(aes(x = metodo, y = sd_median)) +
geom_hline(yintercept = sqrt(0.1367),
size = 1.1, linetype = 2) # DasGupta, pag. 473
true <- sqrt(0.1367)
mean((boot_results[,1] - true)^2)
#> [1] 0.04441418
mean((boot_results[,2] - true)^2)
#> [1] 0.03534701Se observa que el suavizado introduce una ligera tendencia a sobreestimar y a cambio reduce la varianza. El ECM resulta ser inferior, aunque la mejora es modesta. El parámetro de suavizado se ha fijado en \(h=0.5\). Una discusión de cuándo conviene suavizar y cuánto al aplicar bootstrap se puede encontrar en de Angelis y Young (1992).
Ejercicio 7.1 Sea \(X^*\) una v.a. cuya función de densidad es un estimador de núcleo basado en \(X_1,\ldots,X_n\) con parámetro de suavizado \(h\) y núcleo gaussiano.
Calcula la esperanza y la varianza de \(X^*\) (ambas condicionadas a la muestra).
Ejercicio 7.2 Los datos del fichero lipidos.txt corresponden a la concentración de colesterol y triglicéridos (mg/dl) en pacientes evaluados por tener un dolor en el pecho. De estos pacientes, 51 no presentaron evidencia de enfermedad cardiaca mientras que 320 sí la presentaron.
Ejercicio 7.3 Sea \(\hat{f}(x)\) un estimador de núcleo con parámetro de suavizado \(h\) y núcleo simétrico \(K\).
Demuestra que: \[ \int {\hat{f}}^2(x)dx = \frac{1}{n^2h}\sum_{i=1}^n\sum_{j=1}^n K\ast K \left(\frac{X_j-X_i}{h}\right), \] donde \(K\ast K(u)=\int K(u-y)K(y)dy\) es la convolución de \(K\).
Ejercicio 7.4 Sea \(F_n\) la función de distribución empírica correspondiente a una muestra \(X_1,\ldots,X_n\) y sea \(X_1^*,\ldots,X_n^*\) una remuestra extraída de \(F_n\). Se calcula un estimador de núcleo \(\hat{f}_h^*(x)\) a partir de \(X_1^*,\ldots,X_n^*\) y otro \(\hat{f}_h(x)\) a partir de \(X_1,\ldots,X_n\).
Demuestra que \[ \mbox{E}\left[(\hat{f}_h^*(x)-\hat{f}_h(x)) | X_1,\ldots, X_n \right]=0. \] ¿Qué implica este resultado sobre el uso del bootstrap para estimar el sesgo de los estimadores de núcleo?
Ejercicio 7.5 Sean \(X_1,\ldots,X_n\) v.a.i.i.d. de una distribución con función de densidad \(f\) continua. Se considera el estimador de núcleo \(\hat{f}\) con núcleo rectangular \(K(x) = \mathbb{I}_{[-1/2,1/2]}(x)\) y parámetro de suavizado \(h\).
Ejercicio 7.6 Considera una variable aleatoria con distribución beta de parámetros \(\alpha=3\), \(\beta=6\).
Una introducción ya clásica de los estimadores de núcleo de la función de densidad es Silverman (1986). Otros dos libros recomendables sobre el tema son Scott (2015), y Wand y Jones (1994). Chacón y Duong (2018) es un libro de nivel más avanzado y centrado en el caso multivariante.