在尝试将某些代码从Matlab移植到R时，我遇到了一个问题。该代码的要旨是生成2D内核密度估计，然后使用该估计进行一些简单的计算。在Matlab中，使用函数ksdensity2d.m完成KDE计算。在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不？)
最佳答案
在镜像的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()的帮助文件，这建议查看bandwidth的帮助文件。说的是:
但为什么？
density()表示
width参数是为了与S(R的前身)兼容而存在。 source中的ojit_a中的注释为:## 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.
默认值为高斯。当未指定
density()参数而bw为时，将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
但是，因为
width显然是为了与S兼容而编写的(并且我想它最初是为FOR S编写的，因为它在MASS中)，所以一切最终都除以四。翻到MASS的相关部分(第126页)后，似乎他们可能已经选择了四个来在数据的平滑度和保真度之间取得平衡。总之，我的猜测是
kde2d()除以四以与MASS的其余部分(以及其他最初为S编写的东西)保持一致，并且您的处理方式看起来还不错。关于r - 使用kde2d(R)和ksdensity2d(Matlab)生成的2D KDE的差异，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/30626400/