我正在尝试从R中的ks包估计kde对象的95％轮廓的面积。
如果使用ks包中的示例数据集，则将按以下方式创建内核对象：
library(ks)
data(unicef)
H.scv <- Hscv(x=unicef)
fhat <- kde(x=unicef, H=H.scv)
我可以使用plot函数轻松绘制25％，50％，75％的轮廓：
plot(fhat)
但是我想估计轮廓内的面积。
我在这里看到了类似的问题，但提出的答案并不能解决问题。
在我的实际应用中，我的数据集是动物的时间序列，我想使用二元正态核来测量该动物的家园范围。我使用ks包是因为它允许使用诸如插件和平滑交叉验证之类的方法来估计内核分发的带宽。
任何帮助将非常感激！
这有两种方法。它们在概念上都相当复杂，但实际上在代码上却非常简单。
fhat <- kde(x=unicef, H=H.scv,compute.cont=TRUE)
contour.95 <- with(fhat,contourLines(x=eval.points[[1]],y=eval.points[[2]],
z=estimate,levels=cont["95%"])[[1]])
library(pracma)
with(contour.95,polyarea(x,y))
# [1] -113.677
library(sp)
library(rgeos)
poly <- with(contour.95,data.frame(x,y))
poly <- rbind(poly,poly[1,]) # polygon needs to be closed...
spPoly <- SpatialPolygons(list(Polygons(list(Polygon(poly)),ID=1)))
gArea(spPoly)
# [1] 113.677
说明
首先，该kde(...)函数返回一个kde对象，该对象是包含9个元素的列表。您可以在文档中阅读有关此内容的信息，也可以str(fhat)在命令行中键入内容；或者，如果使用的是RStudio（强烈建议使用），则可以通过fhat在“环境”选项卡中展开对象来查看此信息。
元素之一是$eval.points，评估内核密度估计的点。默认值为151个等距点。$eval.points本身就是2个向量的列表。因此，fhat$eval.points[[1]]代表“低于5岁以下” fhat$eval.points[[2]]的得分，代表“ Ave life exp”的得分。
另一个元素是$estimate，它具有内核密度的z值，在x和y的每种组合下进行评估。所以$estimate是一个151 X 151矩阵。
如果你打电话kde(...)跟compute.cont=TRUE，你得到的结果的附加元素：$cont，它包含在Z值$estimate对应于每一个百分点，从1％到99％。
因此，您需要提取与95％轮廓相对应的x和y值，然后使用该值计算面积。您将按照以下步骤进行操作：
fhat <- kde(x=unicef, H=H.scv,compute.cont=TRUE)
contour.95 <- with(fhat,contourLines(x=eval.points[[1]],y=eval.points[[2]],
z=estimate,levels=cont["95%"])[[1]])
现在，contour.95x和y值对应于的95％轮廓fhat。有（至少）两种方法可以到达该区域。一个使用pracma包并直接计算。
library(pracma)
with(contour.95,polyarea(x,y))
# [1] -113.677
负值的原因与x和y的顺序有关：polyarea(...)将多边形解释为“孔”，因此它具有负面积。
另一种方法是使用rgeosGIS程序包中的面积计算例程。不幸的是，这要求您首先将坐标转换为“ SpatialPolygon”对象，这有点让人难以忍受。尽管如此，它也很简单。
library(sp)
library(rgeos)
poly <- with(contour.95,data.frame(x,y))
poly <- rbind(poly,poly[1,]) # polygon needs to be closed...
spPoly <- SpatialPolygons(list(Polygons(list(Polygon(poly)),ID=1)))
gArea(spPoly)
# [1] 113.677