如何使用 'ks' 和 'rgl'从 R 中构建的 3D 核密度图中提取值

Question

我一直在使用'ks'包和'rgl'包来生成3D核密度估计和这些的3D图。第一部分的效果很好(下面是一个简单的例子)。我无法弄清楚的是，如果有可能提取用于构建内核的给定xyz位置的内核值。换句话说，提取3D图中点的值，类似于"光栅"包中用于2D表面的提取命令。有没有人有类似的经验可以给我指出正确的方向?谢谢。dj

library("rgl")
library("ks")
# call the plug-in bandwidth estimator
H.pi <- Hpi(b,binned=TRUE) ## b is a matrix of x,y,z points
# calculate the kernel densities
fhat2 <- kde(b, H=H.pi)
#plot the 50% and 95% kernels in gray and blue
plot(fhat2,cont=c(50,95),colors=c("gray","blue"),drawpoints=TRUE
    ,xlab="", ylab="", zlab="",size=2, ptcol="white", add=FALSE, box=TRUE, axes=TRUE) 


#Structure of fhat2. Original df consists of ~6000 points.  
List of 9
 $ x          : num [1:6173, 1:3] -497654 -497654 -497929 -498205 -498205 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:6173] "50" "57" "70" "73" ...
  .. ..$ : chr [1:3] "x" "max_dep" "y"
$ eval.points:List of 3
  ..$ : num [1:51] -550880 -546806 -542733 -538659 -534586 ...
  ..$ : num [1:51] -7.9 -4.91 -1.93 1.06 4.05 ...
  ..$ : num [1:51] -376920 -374221 -371522 -368823 -366124 ...
$ estimate   : num [1:51, 1:51, 1:51] 0 0 0 0 0 ...
$ H          : num [1:3, 1:3] 3.93e+07 -2.97e+03 8.95e+06 -2.97e+03 2.63e+01 ...
$ gridtype   : chr [1:3] "linear" "linear" "linear"
$ gridded    : logi TRUE
$ binned     : logi FALSE
$ names      : chr [1:3] "x" "max_dep" "y"
$ w          : num [1:6173] 1 1 1 1 1 1 1 1 1 1 ...
 - attr(*, "class")= chr "kde"

Answer 1

试试这个

## from ?plot.kde
library(ks)
library(MASS)
 data(iris)
 ## trivariate example
 fhat <- kde(x=iris[,1:3])
## this indicates the orientation
image(fhat$eval.points[[1]], fhat$eval.points[[2]], apply(fhat$estimate, 1:2, sum))
points(fhat$x[,1:2])
library(raster)
## convert to RasterBrick from raw array
## with correct orientation relative to that of ?base::image
b <- brick(fhat$estimate[,ncol(fhat$estimate):1,], 
    xmn = min(fhat$eval.points[[1]]), xmx = max(fhat$eval.points[[1]]), ymn = min(fhat$eval.points[[2]]), ymx = max(fhat$eval.points[[2]]), 
    transpose = TRUE)
## check orientation
plot(calc(b, sum))
points(fhat$x[,1:2])

现在我们很高兴，因为栅格功率很好。

plot(b)
## note this is a matrix with nrows = nrow(fhat$x), ncols = nlayers(b)
extract(b, fhat$x[,1:2])

Answer 2

答案也可能在eval.points中。看起来您可以在这里输入您自己的点，因此您可以潜在地输入用于构建kde的点或一组全新的点。