rkernel-density

How do I divide one kernel density estimate by another?


I am trying to divide one KDE by another KDE to produce a continuous line for which the Y value at any point X is equal to the ratio of the values of the two initial KDEs at that point X.

Below is a reproducible example in the required format. I want to divide d2 by d1, such that a third line is generated with X values from 4 to 8 and Y values ranging from zero (for all X values for which the Y component of d2 is zero) to approximately 2 (where the Y component of d2 is approximately twice that of d1, as occurs at approx. an X value of 7).

set.seed(1)

d1 <- sample(1:10, 30, replace = TRUE)
d2 <- sample(4:8, 30, replace = TRUE)

d1.d <- density(d1, bw = 1)
d2.d <- density(d2, bw = 1)

plot(d1.d, 
     main = "KDE of d1 (black) and d2 (blue)",
     ylim = c(0,0.25))
lines(d2.d, col = "blue")

enter image description here

How can I accomplish this?


Solution

  • Defining from= and to= in both densityes from their joint range enables to do arithmetic operations in the y values.

    > rg <- range(c(d1, d2))
    > 
    > d1.d <- density(d1, bw = 1, from=rg[1], to=rg[2], n=512)
    > d2.d <- density(d2, bw = 1, from=rg[1], to=rg[2], n=512)
    > 
    > ratio <- d2.d$y/d1.d$y
    > 
    > par(mar=c(5, 4, 4, 3))
    > plot(d1.d, main = "KDE of d1 (black) and d2 (blue)", 
    +      ylim=range(c(d1.d$y, d2.d$y)))
    > lines(d2.d, col=4)
    > fac <- 1/20
    > lines(d1.d$x, ratio*fac, type='l', col=2, lty=2)
    > axis(4, axTicks(2), labels=axTicks(2)/fac, col=2, col.axis=2)
    > mtext('ratio', 4, 2, col=2)
    > legend('topleft', legend=c('d1', 'd2', 'ratio'), col=c(1, 4, 2), lty=c(1, 1, 2))
    

    enter image description here

    Data:

    > dput(d1)
    c(9L, 4L, 7L, 1L, 2L, 7L, 2L, 3L, 1L, 5L, 5L, 10L, 6L, 10L, 7L, 
    9L, 5L, 5L, 9L, 9L, 5L, 5L, 2L, 10L, 9L, 1L, 4L, 3L, 6L, 10L)
    > dput(d2)
    c(5L, 7L, 7L, 7L, 5L, 7L, 4L, 4L, 7L, 4L, 5L, 6L, 5L, 5L, 8L, 
    5L, 4L, 6L, 6L, 7L, 6L, 4L, 7L, 8L, 4L, 4L, 7L, 8L, 8L, 7L)