Search code examples
roptimizationareaellipse

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


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


Solution

  • 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