使用 kde2d (R) 和 ksdensity2d (Matlab) 生成的 2D KDE 的差异

Question

在尝试将一些代码从 Matlab 移植到 R 时，我遇到了运行问题。代码的要点是生成二维核密度估计，然后使用该估计进行一些简单的计算。在 Matlab 中，KDE 计算是使用函数 ksdensity2d.m 完成的。在 R 中，KDE 计算是使用 MASS 包中的 kde2d 完成的。所以假设我想计算 KDE 并只添加值（这不是我打算做的，但它可以达到这个目的）。在 R 中，这可以通过

    library(MASS)
    set.seed(1009)
    x <- sample(seq(1000, 2000), 100, replace=TRUE)
    y <- sample(seq(-12, 12), 100, replace=TRUE)
    kk <- kde2d(x, y, h=c(30, 1.5), n=100, lims=c(1000, 2000, -12, 12))
    sum(kk$z)

给出的答案是 0.3932732。当在 Matlab 中使用 ksdensity2d 使用相同的精确数据和条件时，答案是 0.3768。通过查看 kde2d 的代码，我注意到带宽除以 4

    kde2d <- function (x, y, h, n = 25, lims = c(range(x), range(y))) 
    {
    nx <- length(x)
    if (length(y) != nx) 
     stop("data vectors must be the same length")
    if (any(!is.finite(x)) || any(!is.finite(y))) 
     stop("missing or infinite values in the data are not allowed")
    if (any(!is.finite(lims))) 
     stop("only finite values are allowed in 'lims'")
    n <- rep(n, length.out = 2L)
    gx <- seq.int(lims[1L], lims[2L], length.out = n[1L])
    gy <- seq.int(lims[3L], lims[4L], length.out = n[2L])
    h <- if (missing(h)) 
    c(bandwidth.nrd(x), bandwidth.nrd(y))
    else rep(h, length.out = 2L)
    if (any(h <= 0)) 
     stop("bandwidths must be strictly positive")
    h <- h/4
    ax <- outer(gx, x, "-")/h[1L]
    ay <- outer(gy, y, "-")/h[2L]
    z <- tcrossprod(matrix(dnorm(ax), , nx), matrix(dnorm(ay), 
     , nx))/(nx * h[1L] * h[2L])
    list(x = gx, y = gy, z = z)
    }

一个简单的检查，看看带宽的差异是否是导致结果差异的原因

    kk <- kde2d(x, y, h=c(30, 1.5)*4, n=100, lims=c(1000, 2000, -12, 12))
    sum(kk$z)

给出 0.3768013（与 Matlab 答案相同）。

那么我的问题是：为什么 kde2d 将带宽除以四？（或者为什么不使用 ksdensity2d？）

Answer 1

在镜像 github source 处，第 31-35 行：

if (any(h <= 0))
    stop("bandwidths must be strictly positive")
h <- h/4                            # for S's bandwidth scale
ax <- outer(gx, x, "-" )/h[1L]
ay <- outer(gy, y, "-" )/h[2L]

和 kde2d(), which suggests looking at the help file for bandwidth 的帮助文件。那就是：

...which are all scaled to the width argument of density and so give answers four times as large.

但是为什么呢？

density() says that the width argument exists for the sake of compatibility with S (the precursor to R). The comments in the source 对于 density() 阅读：

## S has width equal to the length of the support of the kernel
## except for the gaussian where it is 4 * sd.
## R has bw a multiple of the sd.

默认为高斯分布。当 bw 参数未指定且 width 是时，width 被替换为，例如

library(MASS)

set.seed(1)
x <- rnorm(1000, 10, 2)
all.equal(density(x, bw = 1), density(x, width = 4)) # Only the call is different

但是，因为 kde2d() 显然是为了与 S 保持兼容而编写的（我想它最初是为 S 编写的，因为它在 MASS 中），所有内容最终都除以四。翻到 MASS 书的相关部分（大约第 126 页）后，似乎他们可能选择了四个以在数据的平滑度和保真度之间取得平衡。

总而言之，我的猜测是 kde2d() 除以四以与 MASS 的其余部分（以及最初为 S 编写的其他内容）保持一致，并且您处理事情的方式看起来不错.

使用 kde2d (R) 和 ksdensity2d (Matlab) 生成的 2D KDE 的差异

Difference in 2D KDE produced using kde2d (R) and ksdensity2d (Matlab)

matlab

r

kernel-density