density_2nd_deriv.R

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

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

d <- density(ran, bw=5) # density uses fast fourier transform; super-fast


# brute-force kernel density estimate (slow! say 2000 ms vs 5 ms)
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)
}  
  
k <- mydensity(ran, d$x, bw=5)

# superficially, they look the same
plot(d$x, d$y, type="l")
lines(k$x, k$y, col="red", lty=2)

# but there are big differences
plot(d$x, d$y - k$y, type="l")

# 1st and 2nd derivatives using the density result
library(ggplot2)
first_derivative <- with(d, differentiate(x, y))
(p1 <- ggplot(first_derivative, aes(x, dyByDx)) + geom_line())

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

# 1st and 2nd derivatives using mydensity
first_derivative_b <- with(k, differentiate(x, y))
(p1b <- ggplot(first_derivative_b, aes(x, dyByDx)) + geom_line())

second_derivative_b <- with(first_derivative_b, differentiate(x, dyByDx))
(p2b <- p1 %+% second_derivative_b + ylab("d2yByDx2"))


# plot the two together
par(las=1, mar=c(5.1, 6.1, 0.6, 0.6))
plot(second_derivative, type="l", col="blue", ylab="")
lines(second_derivative_b, col="red", lwd=2)
title(ylab="estimated 2nd derivative", line=4.5)

Code related to my answer to this stackexchange question.