7
$\begingroup$

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

the first derivative looks smooth, as expected

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

the second derivative is unexpectedly noisy

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 second derivative created directly from the probability density function is smooth

$\endgroup$
3
  • $\begingroup$ The noise effect doesn't happen if you start with a synthetic probability density function created with dnorm. $\endgroup$ May 7, 2014 at 11:42
  • 4
    $\begingroup$ 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. $\endgroup$
    – whuber
    May 7, 2014 at 15:00
  • $\begingroup$ @whuber I clearly recall such a thread, but I can't find it either. $\endgroup$
    – Glen_b
    May 7, 2014 at 22:46

2 Answers 2

6
$\begingroup$

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. enter image description here

$\endgroup$
1
$\begingroup$

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. enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.