2013-05-14 2 views
4

Мне было интересно, если у кого-то была идея о том, как рассчитать синюю заштрихованную область внутри моего эллипса доверия. Любые предложения или места для поиска очень ценятся. Кроме того, я надеюсь найти общую формулу, так как эллипс не обязательно должен лежать в этой области в приложении (т. Е. Эллипс мог бы быть больше). Вот мой код для моей картины, если это поможет:Вычисление области доверительного эллипса в определенной области пространства

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

+1

Общий ответ, чтобы сделать некоторые раздражающие интегралы (Mathematica/Maple/Sage/возможно Yacas может помочь); альтернативно, вы можете сделать интеграцию в Монте-Карло ... –

+0

Существуют ли простые способы выборки точек внутри эллипса? –

+0

Возможно, но еще проще пробовать точки внутри синей прямоугольной области (т. Е. Разбивая ее на три подзоны известной относительной области). –

ответ

4

Если вы можете преобразовать эллипс и синюю область для SpatialPolygons объектов, то это подпруга, используя функцию из rgeos пакета , чтобы рассчитать их площадь, площадь их пересечения и все другие виды интересных величин.

К сожалению, получение объектов в надлежащей форме требует немного тяжелого подъема:

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