5

在此处输入图像描述我正在尝试找到一个点位于椭圆内的概率?例如,如果我在 95% 的椭圆体区域中绘制 300 个数据集的双变量数据 (x,y),我如何计算 300 个点中有多少次落在椭圆内?

这是我正在使用的代码

   library(MASS)
   seed<-1234
   x<-NULL
   k<-1
   Sigma2 <- matrix(c(.72,.57,.57,.46),2,2)
   Sigma2
   rho <- Sigma2[1,2]/sqrt(Sigma2[1,1]*Sigma2[2,2])
   rho
   eta1<-replicate(300,mvrnorm(k, mu=c(-1.59,-2.44), Sigma2))

   library(car)
   dataEllipse(eta1[1,],eta1[2,], levels=c(0.05, 0.95))

谢谢你的帮助。

4

1 回答 1

5

我不明白为什么人们会跳上 OP。在上下文中,这显然是一个编程问题:它是关于在给定椭圆内获取数据点的经验频率,而不是理论概率。OP 甚至发布了代码和图表,显示了他们试图获取的内容。

可能是他们不完全理解 95% 椭圆背后的统计理论,但他们没有问这个问题。此外,像这样制作图和计算频率是掌握理论的好方法。

无论如何,这里有一些代码回答了如何计算通过正态分布获得的椭圆内的点的狭义问题(这是基础dataEllipse)。这个想法是通过主成分将您的数据转换为单位圆,然后获取原点一定半径内的点。

within.ellipse <- function(x, y, plot.ellipse=TRUE)
{
    if(missing(y) && is.matrix(x) && ncol(x) == 2)
    {
        y <- x[,2]
        x <- x[,1]
    }

    if(plot.ellipse)
        dataEllipse(x, y, levels=0.95)

    d <- scale(prcomp(cbind(x, y), scale.=TRUE)$x)
    rad <- sqrt(2 * qf(.95, 2, nrow(d) - 1))
    mean(sqrt(d[,1]^2 + d[,2]^2) < rad)
}

也有人评论说,一个 95% 的数据椭圆包含 95% 的定义数据。这当然不是真的,至少对于正规理论椭圆来说是这样。如果您的分布特别差,随着样本量的增加,覆盖频率甚至可能不会收敛到假设水平。考虑一个广义的帕累托分布,例如:

library(evd) # for rgpd

# generalised pareto has no variance for shape > 0.5
z <- sapply(1:1000, function(...) within.ellipse(rgpd(100, shape=5), rgpd(100, shape=5), FALSE))
mean(z)
[[1] 0.97451

z <- sapply(1:1000, function(...) within.ellipse(rgpd(10000, shape=5), rgpd(10000, shape=5), FALSE))
mean(z)
[1] 0.9995808
于 2013-07-13T09:10:12.363 回答