计算特定空间区域中置信椭圆的面积

问题描述

我想知道是否有人对如何计算置信椭圆内的蓝色阴影区域有所了解.任何建议或寻找的地方，我们将不胜感激.另外，我希望找到一个通用公式，因为在应用中椭圆不必一定位于该区域内(即椭圆可能更大).这是我的图片代码(如果有帮助的话):

I was wondering if someone had an idea on how to calculate the blue shaded area inside my confidence ellipse. Any suggestions or places to look are greatly appreciated. Also, I am hoping to find a general formula since the ellipse does not necessarily have to lie in that region in application (i.e., the ellipse could have been bigger). Here is my code for my picture if it helps:

library(car)

x = c(7,4,1)
y = c(4,6,7)

plot(x,y,xlim=c(0,10),ylim=c(0,10))
rect(x,y,max(x)+100000,max(y)+100000,col="lightblue",border=NA)
points(x,y,col="red",pch=19)

ellipse(center=c(3.5,5),shape=matrix(c(1,.5,.5,2),nrow=2),radius=3,col="green")

推荐答案

如果您可以将椭圆和蓝色区域转换为SpatialPolygons对象，那么使用 rgeos 中的函数就可以了打包，以计算它们的面积，相交的面积以及所有其他有趣的数量.

If you can convert the ellipse and the blue area to SpatialPolygons objects, then it's a cinch, using functions from the rgeos package, to calculate their area, the area of their intersection, and all other sorts of interesting quantities.

不幸的是，以适当的形式获取对象需要一些繁重的工作:

Unfortunately, getting the objects in the proper form requires a bit of heavy lifting:

library(car)
library(sp)
library(rgeos)

## Function for creating a SpatialPolygons object from data.frame of coords
xy2SP <- function(xy, ID=NULL) {
    if(is.null(ID)) ID <- sample(1e12, size=1)
    SpatialPolygons(list(Polygons(list(Polygon(xy)), ID=ID)),
                    proj4string=CRS("+proj=merc"))
}

## Ellipse coordinates
plot.new()  # Needed by ellipse()
ell <- ellipse(center=c(3.5,5),shape=matrix(c(1,.5,.5,2),nrow=2),radius=3)
dev.off()   # Cleaning up after plot.new()

## Three rectangles' coordinates in a length-3 list
x <- c(7,4,1)
y <- c(4,6,7)
mx <- max(x) + 1e6
my <- max(y) + 1e6
rr <- lapply(1:3, function(i) {
    data.frame(x = c(x[i], x[i], mx, mx, x[i]),
               y = c(y[i], my, my, y[i], y[i]))
})


## Make two SpatialPolygons objects from ellipse and merged rectangles
ell <- xy2SP(ell)
rrr <- gUnionCascaded(do.call(rbind, lapply(rr, xy2SP)))

## Find area of their intersection
gArea(gIntersection(ell, rrr))
# [1] 10.36296

这篇关于计算特定空间区域中置信椭圆的面积的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！