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

【问题标题】：Difference in 2D KDE produced using kde2d (R) and ksdensity2d (Matlab)使用 kde2d (R) 和 ksdensity2d (Matlab) 生成的 2D KDE 的差异
【发布时间】：2015-08-18 00:50:48
【问题描述】：

在尝试将一些代码从 Matlab 移植到 R 时，我遇到了问题。代码的要点是生成一个 2D 核密度估计，然后使用该估计进行一些简单的计算。在 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？）

【问题讨论】：

标签： r matlab kernel-density

【解决方案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() 的帮助文件，它建议查看bandwidth 的帮助文件。那就是：

...它们都缩放到密度的宽度参数，因此给出答案是原来的四倍。

但是为什么呢？

density() 表示存在 width 参数是为了与 S（R 的前身）兼容。 source for density() 中的 cmets 读取：

## 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 书的相关部分后（大约 p.126），他们似乎选择了四个来在数据的平滑度和保真度之间取得平衡。

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

【讨论】：