R:二元(或多元)核密度的概率/数值积分



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

n <- 70
x <- rnorm(n)
library(ks)
f_kde <- kde(x) 

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

p_kde <- pkde(x, f_kde)

这是在ks中使用Simpson规则进行数值积分来完成的。不幸的是,他们只在1d的情况下实现了这一点。在二元情况下,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) 

我在stackoverflow中找不到任何包或帮助搜索来解决R中的这个问题(有一些关于Python的建议,但我想把它保留在R中(。任何代码行或包装推荐都将不胜感激。尽管我最感兴趣的是双变量情况,但任何关于多变量情况的想法都会受到赞赏。

kde允许多维核估计,因此我们可以使用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

最新更新