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
参数在足够小的 dx
和 dy
步上计算 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
我正在使用包 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
参数在足够小的 dx
和 dy
步上计算 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