¿Por qué un núcleo de estimación de la densidad de una distribución normal vector tiene un no-suave segunda derivada?

Question

He tenido algunos distribución normal de los datos:

mu <- 3
sigma <- 5
x <- rnorm(1e5, mu, sigma)

Me tomó un núcleo de estimación de la densidad bastante con un alto ancho de banda:

kernel_density_of_x <- density(x, bw = 5)

Entonces me diferenciadas:

differentiate <- function(x, y)
{
  diffOfX <- diff(x)
  data.frame(
    x      = x[-length(x)] + (diffOfX / 2), 
    dyByDx = diff(y) / diffOfX
  )
}

first_derivative <- with(kernel_density_of_x, differentiate(x, y))

Esto parecía tan esperado:

library(ggplot2)
(p1 <- ggplot(first_derivative, aes(x, dyByDx)) + geom_line())

Cuando me diferenciadas de nuevo, me esperaba otra curva suave, pero vi a un extraño efecto cíclico.

second_derivative <- with(first_derivative, differentiate(x, dyByDx))
(p2 <- p1 %+% second_derivative + ylab("d2yByDx2"))

He intentado un par de opciones diferentes para la kernel argumento, pero el ruido persiste.
Dejar caer el ancho de banda, por ejemplo, 0.5 dio una menor frecuencia de ruido que dominó la trama (lo que es absurdo).

Bajando el número de puntos de muestreo hacia abajo desde n = 512 a n = 32 parado el tema, pero que la causa de otros problemas.

¿Por qué este efecto produce? Es un artefacto de la density función, o he hecho algo tonto?

---

Podemos dibujar el gráfico mediante la utilización de la función de densidad de probabilidad de la distribución normal que x fue generado a partir de ver la forma que yo esperaba:

xx <- seq.int(-20, 20, 0.1)
pdf_of_xx <- dnorm(xx, mu, sigma)
first_derivative_of_xx <- differentiate(xx, pdf_of_xx)
second_derivative_of_xx<- with(first_derivative_of_xx, differentiate(x, dyByDx))
ggplot(second_derivative_of_xx, aes(x, dyByDx)) + geom_line()

Answer 1

Como whuber comentó, lo que estamos viendo es debido a la aproximación (a través de una transformada rápida de fourier) que density utiliza.

Si se calcula el núcleo de la estimación de la densidad por la fuerza bruta y la comparamos con la estimación de que el density da, vas a ver este tipo de patrón cíclico en las diferencias.

La 2ª diferencias en el resultado de density resultado en una exageración de los efectos. La segunda de las diferencias de la fuerza bruta estimación de densidad de kernel no muestran que el patrón, pero el cálculo es como de 3 órdenes de magnitud más lento.

Pongo el código en este gist. Mi fuerza bruta estimación de densidad de kernel es la siguiente:

mydensity <-
function(dat, x, bw=5) # dat=the data; x=points to calculate density estimate
{
  y <- vapply(x, function(a) mean(dnorm(a, dat, bw)), 1)
  data.frame(x=x, y=y)
}

Aquí una foto de la segunda de las diferencias, el azul de la density y el rojo de la fuerza bruta enfoque.

Answer 2

En sus códigos en el gist se utiliza un duro introducido por el ancho de banda (bw=5) que es muy diferente de un óptimo devuelto por ejemplo como:

library(ks)
h <- hpi(x)

que en realidad es acerca de la h=0.85.

Es interesante para la ejecución de los códigos de donde bw se establece en h. A continuación, la cifra final es la siguiente. Aquí las diferencias no son tan drásticos.