在尝试将一些代码从 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?)
最佳答案
在镜像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 的帮助文件.那就是:
...which are all scaled to the width argument of density and so give answers four times as large.
但为什么呢?
density()表示 width
参数的存在是为了与 S(R 的前身)兼容。 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()
除以 4 以与 MASS 的其余部分(以及最初为 S 编写的其他内容)保持一致,并且您正在进行的方式一切看起来都很好。
关于r - 使用 kde2d (R) 和 ksdensity2d (Matlab) 生成的 2D KDE 的差异,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/30626400/