6.5 Estimación no paramétrica de la densidad en R

Como ya se comentó en la Sección 3.3, en R podemos emplear la función density() del paquete base para obtener una estimación tipo núcleo de la densidad. Los principales parámetros (con los valores por defecto) son los siguientes:

density(x, bw = "nrd0", adjust = 1, kernel = "gaussian", n = 512, from, to)

  • bw: ventana, puede ser un valor numérico o una cadena de texto que la determine (en ese caso llamará internamente a la función bw.xxx() donde xxx se corresponde con la cadena de texto). Las opciones son:

    • "nrd0", "nrd": Reglas del pulgar de Silverman (1986, page 48, eqn (3.31)) y Scott (1992), respectivamente. Como es de esperar que la densidad objetivo no sea tan suave como la normal, estos criterios tenderán a seleccionar ventanas que producen un sobresuavizado de las observaciones.

    • "ucv", "bcv": Métodos de validación cruzada insesgada y sesgada, respectivamente.

    • "sj", "sj-ste", "sj-dpi": Métodos de Sheather y Jones (1991), “solve-the-equation” y “direct plug-in,” respectivamente.

  • adjust: parameto para reescalado de la ventana, las estimaciones se calculan con la ventana adjust*bw.

  • kernel: cadena de texto que determina la función núcleo, las opciones son: "gaussian", "epanechnikov", "rectangular", "triangular", "biweight", "cosine" y "optcosine".

  • n, from, to: permiten establecer la rejilla en la que se obtendrán las estimaciones (si \(n>512\) se emplea fft() por lo que se recomienda establecer n a un múltiplo de 2; por defecto from y to se establecen como cut = 3 veces la ventana desde los extremos de las observaciones).

Utilizaremos como punto de partida el código empleado en la Sección 3.3. Considerando el conjunto de datos precip (que contiene el promedio de precipitación, en pulgadas de lluvia, de 70 ciudades de Estados Unidos).

x <- precip
h <- bw.SJ(x)
npden <- density(x, bw = h)
# npden <- density(x, bw = "SJ")

# plot(npden)
hist(x, freq = FALSE, main = "Kernel density estimation",
     xlab = paste("Bandwidth =", formatC(h)), lty = 2,
     border = "darkgray", xlim = c(0, 80), ylim = c(0, 0.04))
lines(npden, lwd = 2)
rug(x, col = "darkgray")

Alternativamente podríamos emplear implementaciones en otros paquetes de R. Uno de los más empleados es ks (Duong, 2019), que admite estimación incondicional y condicional multidimensional. También se podrían emplear los paquetes KernSmooth (Wand y Ripley, 2019), sm (Bowman y Azzalini, 2019), np (Tristen y Jeffrey, 2019), kedd (Guidoum, 2019), features (Duong y Matt, 2019) y npsp (Fernández-Casal, 2019), entre otros.