CODEWITHSUNDEEP
rellipse
user2005253
user2005253

Reputation:

Getting the ellipse function to match the dataEllipse function in R

So I am sampling from a multivariate normal distribution in R, and am trying to figure out how to calculate its 95% confidence ellipse using the ellipse() function in the package car.

Here is some code I am running:

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)

So basically I want the ellipse command to replicate the dataEllipse command. If anyone has any suggestions that would be greatly appreciated!

Edit: Using Dwins code and combining it within my own:

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

enter image description here

So as you can see, the ellipses are still not the same...

Upvotes: 3

Views: 5045

Answers (2)

user2175871
user2175871

Reputation: 121

Since this post is still getting views, I'll provide the actual answer. The last three lines of this code snippet replicates car::dataEllipse exactly:

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

Note that both car::dataEllipse and car::ellipse return the coordinates of the points silently, so one can confirm that the points are indeed equal.

Upvotes: 2

IRTFM
IRTFM

Reputation: 263332

I'm guessing (although you should not have made me do so) that rmvnorm is from 'mixtools' and it was loaded after 'car'. I do not think the sqrt() function is needed since the argument to ellipse is supposed to be a covariance matrix. Also at the moment it is plotting but you cannot see it, because you didn't color it red (or anything). Furthermore both 'mixtools' and 'car' have ellipse functions, so if you want the car-version ( which does have a radius argument unlike the mixtools version) then you need to call it with the double colon convention:

library(car); library(mixtools)
car::ellipse(mu, sigma*qchisq(.05,2), col="red", 
                 radius=sqrt(2 * qf(.975, 2, 9998)) )

Upvotes: 2

Related Questions