After 3 weeks and close to 2 days of constant thinking about density functions, I have finally cracked the density function in base R.My head was spinning in all directions and I could not get my hands on cracking it.

Finally I did what any failed person would do to learn something. Start from scratch

Firstly I wrote my own R code to do kernel density estimation

> gauss.wt <- function(x) {
+     1/(sqrt(2 * pi)) * exp(-x^2/2)
+ }
> set.seed(1977)
> n <- 512
> obs.c <- sample(c(1, 2), n, T)
> mus <- c(3, 7)
> sigs <- c(1, 1)
> x <- sapply(obs.c, function(x) rnorm(1, mus[x], sigs[x]))
> n <- length(x)
> h <- bw.nrd0(x)
> prob <- sapply(x, function(z) (1/(n * h)) * sum(gauss.wt((z -
+     x)/h)))
> z <- cbind(x, prob)
> z <- z[order(x), ]
> par(mfrow = c(1, 1))
> plot(prob ~ x, type = "l", data = z, xlim = range(x), ylim = c(0,
+     0.2), col = "red", lwd = 3)

Cracked-002.jpg

Ok, it looks like KDE

Now I wanted to compare with base R output

> res <- density(x)
> prob1 <- approx(res$x, res$y, z[, 1])

So, I was able to get the density values at all the relevant points.

> par(mfrow = c(1, 1))
> plot(prob ~ x, type = "l", data = z, xlim = range(x), ylim = c(0,
+     0.2), col = "grey", lwd = 8)
> points(z[, 1], prob1$y, col = "red", lty = "dashed", cex = 0.2,
+     pch = 19)

Cracked-004.jpg

This was the first mini achievement. I could manually replicate the density function

Now the next task was to customize density function of R and code my own. I download the R source code, read through the C file that was being used in the density function and replicated the C code in R.

> density.rad <- function(x, bw = "nrd0", n = 512, cut = 3) {
+     x <- as.vector(x)
+     weights <- rep.int(1/n, n)
+     n <- max(n, 512)
+     if (n > 512)
+         n <- 2^ceiling(log2(n))
+     bw <- bw.nrd0(x)
+     from <- min(x) - cut * bw
+     to <- max(x) + cut * bw
+     lo <- from - 4 * bw
+     up <- to + 4 * bw
+     BinDist <- function(x, weights, lo, up, n) {
+         y <- numeric(2 * n)
+         xdelta <- (up - lo)/n
+         for (i in seq_along(x)) {
+             xpos <- (x[i] - lo)/xdelta
+             ix <- floor(xpos)
+             fx <- xpos - ix
+             wi <- weights[i]
+             y[ix] <- y[ix] + (1 - fx) * wi
+             y[ix + 1] <- y[ix + 1] + (fx) * wi
+         }
+         return(y)
+     }
+     y <- BinDist(x, weights, lo, up, n)
+     kords <- seq.int(0, 2 * (up - lo), length.out = 2L * n)
+     kords[(n + 2):(2 * n)] <- -kords[n:2]
+     kords <- dnorm(kords, sd = bw)
+     kords <- fft(fft(y) * Conj(fft(kords)), inverse = TRUE)
+     kords <- pmax.int(0, Re(kords)[1L:n]/length(y))
+     xords <- seq.int(lo, up, length.out = n)
+     x <- seq.int(from, to, length.out = n)
+     list(x = x, y = approx(xords, kords, x)$y)
+ }
> res <- density.rad(x)
> prob1 <- approx(res$x, res$y, z[, 1])

So, I was able to get the density values at all the relevant points.

> par(mfrow = c(1, 1))
> plot(prob ~ x, type = "l", data = z, xlim = range(x), ylim = c(0,
+     0.2), col = "grey", lwd = 8)
> points(z[, 1], prob1$y, col = "red", lty = "dashed", cex = 0.2,
+     pch = 19)

Cracked-007.jpg

So, finally I have achieved a mini milestone, i.e. understanding the density function at a certain depth.

I am still unclear on the convolution front. Let me take some time and understand that aspect. But for now , let me blog this. After a very long time , I feel that I have cracked something nice. I had been going through a lot of troubles to understand this.