Calculating the area of a confidence ellipse in a certain region of space

Question

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")

enter image description here

The general answer is to do some annoying integrals (Mathematica/Maple/Sage/maybe Yacas could help); alternatively you could do Monte Carlo integration ... — Ben Bolker, May 14 '13 at 22:04
Probably, but it should be even simpler to sample points inside the blue rectangular region (i.e. by breaking it into three sub-areas of known relative area). — Ben Bolker, May 14 '13 at 22:15
Not true, because some of those points will fall outside the ellipse. — , May 14 '13 at 23:02
I suppose if you rotate the whole figure so that the larger axis of the ellipse is vertical, you can integrate (since both the upper half of the ellipse and the curve delineating the blue area will be functions)...? — Frank, May 15 '13 at 00:11

Josh O'Brien · Accepted Answer · 2013-05-15T01:59:54.517

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

Calculating the area of a confidence ellipse in a certain region of space

1 Answers1

Linked