如何将一个核密度估计(KDE)除以另一个？附R语言示例

需求：将d2的核密度估计（KDE）除以d1的KDE，生成一条连续曲线，曲线上任意X点的Y值等于两个KDE在该点的数值之比，X范围限定为4至8，Y值范围从0到约2（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")

实现步骤 #

两个KDE对象的x坐标序列通常不一致，直接做除法会因为维度不匹配报错。核心是先将两个KDE插值到相同的x网格上，再计算比值：

1. 定义统一的x轴范围（4到8），取足够多的点保证曲线连续
2. 对两个KDE的密度值进行插值，得到统一x网格上的对应值
3. 计算比值，处理分母接近0的情况避免NaN/无穷大
4. 绘制比值曲线（可叠加到原图或单独绘制）

完整代码 #

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)

# 1. 定义共同的x网格，4到8之间取1000个点
common_x <- seq(4, 8, length.out = 1000)

# 2. 插值得到统一x上的密度值
d1_interp <- approx(d1.d$x, d1.d$y, xout = common_x)$y
d2_interp <- approx(d2.d$x, d2.d$y, xout = common_x)$y

# 3. 计算比值，分母过小时设为0（避免异常值）
ratio <- ifelse(d1_interp < 1e-6, 0, d2_interp / d1_interp)

# 4. 绘制图像，调整y轴范围以显示比值曲线
plot(d1.d, 
     main = "KDE of d1, d2 and Their Ratio",
     ylim = c(0, 2.2))  # 扩大y轴范围容纳比值
lines(d2.d, col = "blue")
lines(common_x, ratio, col = "red", lwd = 2)

# 添加图例
legend("topright", 
       legend = c("d1 KDE", "d2 KDE", "d2/d1 Ratio"),
       col = c("black", "blue", "red"),
       lwd = c(1, 1, 2))

说明 #

• approx()函数用于线性插值，保证两个KDE在相同x点上有对应密度值
• 加入ifelse(d1_interp < 1e-6, 0, ...)是为了处理d1 KDE接近0的情况，防止出现无穷大或NaN，保证曲线连续
• 调整ylim到c(0,2.2)是因为比值的最大值接近2，需要完整显示曲线

内容的提问来源于stack exchange，提问作者fre1990