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

问题描述

我想知道是否有人知道如何计算置信椭圆内的蓝色阴影区域.任何建议或要看的地方都非常感谢.此外，我希望找到一个通用公式，因为椭圆在应用中不一定位于该区域(即椭圆可能更大).如果有帮助，这是我的图片代码:

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

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