如何将一个核密度估计除以另一个？

Question

我正在尝试将一个KDE除以另一个KDE，以产生一条连续的线，其中任意点X处的Y值等于两个初始KDE在该点X处的值的比率。

下面是一个可复现的示例，按照所需的格式。我想要将d2除以d1，生成一个第三条线，其X值从4到8，Y值范围从零（对于所有Y组件为零的X值）到大约2（当d2的Y组件大约是d1的两倍时，如在X值约为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")

如何完成这个任务？

- overcup

---

我认为plot(d1.d$x, d1.d$y/d2.d$y, type = "l")会产生你想要的线条，但我无法弄清楚如何将density对象与同一图形一起绘制。 - undefined

Answer 1

在两个密度的联合范围中定义from=和to=，可以对y值进行算术运算。

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

数据：

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

- jay.sf

---

@jay.sf，目前来看，比率术语仍然是颠倒的。你是怎么得出“fac.”的值为1/20的呢？虽然它似乎有效，但我无法弄清楚你是如何确定使用这个特定的值的。 - undefined

@overcup 我使用了ratio*fac以及axTicks(2)/fac，这样应该能保持一致。1/20只是看起来不错 :) - undefined

@jay.sf，啊，所以fac主要/完全是为了美观？我以为你使用了一个fac的值，以便缩放比例线，使其下方的面积等于1。 - undefined

要得到一个正确的PDF，你还需要除以积分。不确定这是否是OP想要的，但是你可以通过在定义的区域上进行integrate(approxfun(...), ...)来计算密度的积分，然后将密度除以该数值。 - undefined

@jay.sf，对不起，我完全没有注意到你在图表的右侧使用了一个单独的比例尺。你使用1/20比例尺现在有意义了。 - undefined

Answer 2

你可以使用approxfun()将density()的结果转换为函数。你只需要对其中一个密度估计进行这样的操作，然后可以直接使用另一个密度估计的x和y值。

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)
d2.f <- approxfun(d2.d$x, d2.d$y)

plot(d1.d, 
     main = "KDE of d1 (black) and d2 (blue)",
     ylim = c(0,3))
lines(d2.d, col = "blue")
lines(d1.d$x, d2.f(d1.d$x)/d1.d$y, col = "red")
legend("topleft", legend = c("d1", "d2", "d2/d1"), 
       col = c("black", "blue", "red"), lty=1)

^{2023-11-13创建，使用reprex v2.0.2生成}

- user2554330