3

因此，我从 R 中的多元正态分布中采样，并试图弄清楚如何使用包装车中的 ellipse() 函数计算其 95% 置信椭圆。

这是我正在运行的一些代码：

mu = c(0,0)
sigma = matrix(c(20,0,0,45),nrow=2)

z = rmvnorm(10000,mu,sqrt(sigma))

par(mfrow=c(1,2))
plot(z)
ellipse(mu,sqrt(sigma*qchisq(.05,2)),radius=1)
dataEllipse(z,levels=.95)

所以基本上我希望 ellipse 命令复制 dataEllipse 命令。如果有人有任何建议，将不胜感激！

编辑：使用 Dwins 代码并将其合并到我自己的代码中：

library(car)
library(mvtnorm)

mu = c(0,0)
sigma = matrix(c(20,0,0,45),nrow=2)

z = rmvnorm(10000,mu,sqrt(sigma))

dataEllipse(z,levels=.95)
car::ellipse(mu, sigma*qchisq(.05,2), col="blue", 
                 radius=sqrt(2 * qf(.975, 2, 9998)) )

在此处输入图像描述

所以你可以看到，椭圆仍然不一样......

4

2 回答 2

2

我猜（尽管你不应该让我这样做）rmvnorm来自“mixtools”，它是在“car”之后加载的。我认为不需要 sqrt() 函数，因为 to 的参数ellipse应该是协方差矩阵。目前它正在绘制但您看不到它，因为您没有将它涂成红色（或任何东西）。此外，'mixtools' 和 'car' 都有ellipse功能，所以如果你想要 car-version（它确实有一个与 mixtools 版本不同的 radius 参数），那么你需要使用双冒号约定来调用它：

library(car); library(mixtools)
car::ellipse(mu, sigma*qchisq(.05,2), col="red", 
                 radius=sqrt(2 * qf(.975, 2, 9998)) )
于 2013-05-21T00:52:17.817 回答
2

由于这篇文章仍在获得意见，我将提供实际答案。此代码片段的最后三行car::dataEllipse完全复制：

library(car)
library(mvtnorm)

mu = c(0,0)
sigma = matrix(c(20,0,0,45),nrow=2)
z = rmvnorm(10000,mu,sigma)

dataEllipse(z,levels=.95)

center <- apply(z, 2, mean)
cov_mat <- cov(z)
ellipse(center, cov_mat, col="red", radius=sqrt(2 * qf(.95, 2, 9999)))

请注意，两者都car::dataEllipse默默地car::ellipse返回点的坐标，因此可以确认这些点确实相等。

于 2020-07-29T20:31:03.770 回答