R:双变量(或多变量)核密度的概率/数值积分

R: probability / numerical integral of bivariate (or multivariate) kernel density

我正在使用包 ks 进行核密度估计。这是一个简单的例子:

n <- 70
x <- rnorm(n)

library(ks)
f_kde <- kde(x) 

我实际上对我的输入数据各自的超出概率感兴趣,可以很容易地由具有 f_kde:

的 ks 返回
p_kde <- pkde(x, f_kde)

这是在 ks 中使用辛普森规则通过数值积分完成的。不幸的是,他们只在一维案例中实施了这一点。在双变量情况下,ks 中没有实现任何返回概率的方法:

y <- rnorm(n)
f_kde <- kde(data.frame(x,y))
# does not work, but it's what I am looking for:
p_kde <- pkde(data.frane(x,y), f_kde) 

我无法在 Whosebug 中找到任何包或帮助搜索来解决 R 中的这个问题(存在 Python 的一些建议,但我想将其保留在 R 中)。任何代码行或包推荐表示赞赏。尽管我最感兴趣的是双变量案例,但也欢迎任何有关多变量案例的想法。

kde允许multidimensional kernel estimate,所以我们可以用kde来计算pkde
为此,我们使用 eval.points 参数在足够小的 dxdy 步上计算 kde :这为我们提供了 dx*dy 上的局部密度估计 正方形。
我们验证乘以正方形表面的估计总和几乎等于 1:

library(ks)
set.seed(1)
n <- 10000
x <- rnorm(n)
y <- rnorm(n)
xy <- cbind(x,y)

xmin <- -10
xmax <- 10
dx <- .1

ymin <- -10
ymax <- 10
dy <- .1

pts.x <- seq(xmin, xmax, dx)
pts.y <- seq(ymin, ymax, dy)
pts <- as.data.frame(expand.grid(x = pts.x, y = pts.y))
f_kde <- kde(xy,eval.points=pts)

pts$est <- f_kde$estimate

sum(pts$est)*dx*dy
[1] 0.9998778

您现在可以在 pts 数据框中查询所选区域的累积概率:

library(data.table)
setDT(pts)
# cumulative density
pts[x < 1 & y < 2 , .(pkde=sum(est)*dx*dy)]
        pkde
1: 0.7951228

# average density around a point
tolerance <-.1
pts[pmin(abs(x-1))<tolerance & pmin(abs(y-2))<tolerance, .(kde = mean(est))]
          kde
1: 0.01465478