# Why does a kernel density estimate of a normally distributed vector have a non-smooth second derivative?

I had some normally distributed data:

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


I took a kernel density estimate with a fairly high bandwidth:

kernel_density_of_x <- density(x, bw = 5)


Then I differentiated it:

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))


This looked just as expected:

library(ggplot2)
(p1 <- ggplot(first_derivative, aes(x, dyByDx)) + geom_line()) When I differentiated again, I expected another smooth curve, but I saw an odd cyclical effect.

second_derivative <- with(first_derivative, differentiate(x, dyByDx))
(p2 <- p1 %+% second_derivative + ylab("d2yByDx2")) I tried a few different options for the kernel argument, but the noisiness persisted.
Dropping the bandwidth down to, for example, 0.5 gave a lower frequency noise that dominated the plot (making it nonsense).

Dropping the number of sampling points down from n = 512 to n = 32 stopped the issue, but that causes other problems.

Why does this effect occur? Is it an artifact of the density function, or have I done something silly?

We can redraw the plot using the probability density function of the normal distribution that x was generated from to see the shape I expected:

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() • The noise effect doesn't happen if you start with a synthetic probability density function created with dnorm. – Richie Cotton May 7 '14 at 11:42
• It's an artifact due to approximations made by the underlying R code; the discrete second derivative operator greatly magnifies the small errors. This has been discussed in another thread which I cannot at the moment find ... When you think about it, smoothness of the derivative (of a kernel density estimate as computed in R) isn't even defined in this context, because everything in sight is a discrete set of points, anyway. – whuber May 7 '14 at 15:00
• @whuber I clearly recall such a thread, but I can't find it either. – Glen_b -Reinstate Monica May 7 '14 at 22:46

As whuber commented, what you're seeing is due to the approximation (via a fast fourier transform) that density uses.

If you calculate the kernel density estimate by brute force and compare it to the estimate that density gives, you'll see this sort of cyclic pattern in the differences.

The 2nd differences in the result from density result in an exaggeration of the effect. The second differences from the brute-force kernel density estimate don't show that pattern, but the calculation is like 3 orders of magnitude slower.

I put some code in this gist. My brute-force kernel density estimate is the following:

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)
}


Here's a picture of the second differences, blue from density and red from brute-force approach. In your codes in the gist you used a hard-entered bandwidth (bw=5) which is very different from a optimal one returned for example as:

library(ks)
h <- hpi(x)


which is actually about h=0.85.

It is interesting to execute your codes where bw is set to h. Then the final figure is as below. Here the differences are not so drastic. 