3

我之前发布了一个类似的问题。我试图确定一个点是否在椭圆内。基本上我会生成一些二元法线数据并创建一个椭圆。这是我使用的代码

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

我需要找出这对 (p1[,1],p1[,2]) 是否位于上面的椭圆区域内。

4

2 回答 2

7

dataEllipse将椭圆返回为多边形,因此您可以使用库中的point.in.polygon函数sp来检查点是否在椭圆内:

ell = dataEllipse(rate1, rate2, levels=c(0.05, 0.95)) 
point.in.polygon(rate1, rate2, ell$`0.95`[,1], ell$`0.95`[,2])

当我运行以下代码时...

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

...我得到以下输出

  [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
于 2013-08-09T19:51:12.570 回答
1

只需找到 C 并减去半径

在此处输入图像描述 在此处输入图像描述

于 2013-08-09T20:01:01.073 回答