Reputation: 327
I had posted a similar question earlier. I was trying to determine whether a point lies within an ellipse. Basically I generate some bivariate normal data and create an ellipse. Heres the code I use
library(MASS)
set.seed(1234)
x1<-NULL
x2<-NULL
k<-1
Sigma2 <- matrix(c(.72,.57,.57,.46),2,2)
Sigma2
rho <- Sigma2[1,2]/sqrt(Sigma2[1,1]*Sigma2[2,2])
eta<-replicate(300,mvrnorm(k, mu=c(-2.503,-1.632), Sigma2))
p1<-exp(eta)/(1+exp(eta))
n<-60
x1<-replicate(300,rbinom(k,n,p1[,1]))
x2<-replicate(300,rbinom(k,n,p1[,2]))
rate1<-x1/60
rate2<-x2/60
library(car)
dataEllipse(rate1,rate2,levels=c(0.05, 0.95))
I need to find out whether the pair (p1[,1],p1[,2]) lies within the area of the ellipse above.
Upvotes: 3
Views: 1340
Reputation: 951
dataEllipse
returns the ellipses as polygons, so you could use the point.in.polygon
function from the sp
library to check whether the points are inside the ellipse:
ell = dataEllipse(rate1, rate2, levels=c(0.05, 0.95))
point.in.polygon(rate1, rate2, ell$`0.95`[,1], ell$`0.95`[,2])
When I run the following code...
library(MASS)
set.seed(1234)
x1<-NULL
x2<-NULL
k<-1
Sigma2 <- matrix(c(.72,.57,.57,.46),2,2)
Sigma2
rho <- Sigma2[1,2]/sqrt(Sigma2[1,1]*Sigma2[2,2])
eta<-replicate(300,mvrnorm(k, mu=c(-2.503,-1.632), Sigma2))
p1<-exp(eta)/(1+exp(eta))
n<-60
x1<-replicate(300,rbinom(k,n,p1[,1]))
x2<-replicate(300,rbinom(k,n,p1[,2]))
rate1<-x1/60
rate2<-x2/60
library(car)
ell = dataEllipse(rate1, rate2, levels=c(0.05, 0.95))
library(sp)
point.in.polygon(rate1, rate2, ell$`0.95`[,1], ell$`0.95`[,2])
... I get the following output
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[56] 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[111] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[166] 1 1 1 1 1 1 1 0 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1
[221] 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[276] 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Upvotes: 7