6

I have this set of x and y coordinates:

x<-c(1.798805,2.402390,2.000000,3.000000,1.000000)
y<-c(0.3130147,0.4739707,0.2000000,0.8000000,0.1000000)
as.matrix(cbind(x,y))->d

and I want to calculate the ellipsoid that contains this set of points, I use the function ellipsoidhull() in the package "cluster", and I get:

> ellipsoidhull(d)
'ellipsoid' in 2 dimensions:`
 center = ( 2.00108 0.36696 ); squared ave.radius d^2 =  2`
and shape matrix =
x 0.66590 0.233106
y 0.23311 0.095482
  hence, area  =  0.60406

However it's not obvious to me how I can get from these results, the lengths of the semi-major axes of this ellipse.

Any idea?

Thank you very much in advance.

Tina.

user18441
  • 643
  • 1
  • 7
  • 15
  • 1
    Might be worth reading thru this question http://stackoverflow.com/questions/3417028/ellipse-around-the-data-in-matlab and the links there for some info on relating the eigenvalues of the covariance matrix to the axes of the ellipse. – Carl Witthoft Aug 16 '13 at 17:28
  • user18441 I treat the question as a geometrical problem, so maybe it is worth that you read the link showed in the above comment and use better tools ( statistical tools) to deal with it. – agstudy Aug 16 '13 at 17:39
  • 1
    Even better, the wikipedia page "ellipsoid" defines the shape matrix in terms of the axes (via the eigenvalues), so you should be able to calculate the radii explicitly. – Carl Witthoft Aug 16 '13 at 17:42

2 Answers2

5

You can do this:

exy <- predict(ellipsoidhull(d)) ## the ellipsoid boundary
me <- colMeans((exy))            ## center of the ellipse

Then you compute the minimum and maximum distance to get respectively minor and major axis: 

dist2center <- sqrt(rowSums((t(t(exy)-me))^2))
max(dist2center)     ## major axis
[1] 1.264351
> min(dist2center)   ## minor axis
[1] 0.1537401

EDIT plot the ellipse with the axis:

plot(exy,type='l',asp=1)
points(d,col='blue')
points(me,col='red')
lines(rbind(me,exy[dist2center == min(dist2center),]))
lines(exy[dist2center == max(dist2center),])

enter image description here

Newbie_R
  • 655
  • 7
  • 22
agstudy
  • 119,832
  • 17
  • 199
  • 261
  • Thank you agstudy! to plot them I guess I need to trace a line between the center of the ellipse and the point the furthest away and the closest from it. Any idea of how this could be easily be done? – user18441 Aug 16 '13 at 17:26
  • Heh. An entirely geometry-free solution. Is perhaps `me <- colMeans((exy))` a better center? The interior points in the sample shouldn't contribute, I think. – IRTFM Aug 16 '13 at 17:45
  • @DWin right. I edit my answer. But, For some reasons I get 2 points for the max but only one point for the min. – agstudy Aug 16 '13 at 17:54
  • Maybe you need to look at both min and "next-min" , i.e. the penultimate-minimum if there is such a word? – IRTFM Aug 16 '13 at 18:10
  • one last thing, using the code as it is now, I don't seem to get the semi-axes well plotted. The long semi axis appears to occupy the total length of the long axis and the minor axis fails to be at 90º with the long axis. What can be wrong? – user18441 Aug 16 '13 at 18:18
  • @MikeWise What do you mean? What did you try that don't work? – agstudy Mar 15 '17 at 18:55
  • `x <- predict(ellipsoidhull(matrix(rnorm(120),ncol=2)))` works, but `ncol=3` emits the error "not yet implemented for p>=3" – Mike Wise Mar 15 '17 at 19:02
  • Or is there a way around this? – Mike Wise Mar 15 '17 at 19:02
5

The square of the semi-axes are the eigenvalues of the shape matrix, times the average squared radius.

x <- c(1.798805,2.402390,2.000000,3.000000,1.000000)
y <- c(0.3130147,0.4739707,0.2000000,0.8000000,0.1000000)
d <- cbind( x, y )
library(cluster)
r <- ellipsoidhull(d)
plot( x, y, asp=1, xlim=c(0,4) )
lines( predict(r) )
e <- sqrt(eigen(r$cov)$values)
a <- sqrt(r$d2) * e[1]  # semi-major axis
b <- sqrt(r$d2) * e[2]  # semi-minor axis
theta <- seq(0, 2*pi, length=200)
lines( r$loc[1] + a * cos(theta), r$loc[2] + a * sin(theta) )
lines( r$loc[1] + b * cos(theta), r$loc[2] + b * sin(theta) )
Vincent Zoonekynd
  • 31,893
  • 5
  • 69
  • 78
  • Thank you Vincent. Sorry my ignorance. "a" and "b" the lengths of the axes of the ellipse and therefore to get the semi-axes they need to be divided by 2. Is this correct? – user18441 Aug 16 '13 at 18:33
  • Yeah yeah, steal my commments-answer :-). Very nice details in this code. – Carl Witthoft Aug 16 '13 at 19:02
  • @user18441: the comments in my code were incorrect, these are already the semi-axes -- no need to divide by 2. – Vincent Zoonekynd Aug 16 '13 at 19:21
  • A screen shot would be nice. I spent 2-3 minutes trying to see how it corresponded to the plot below before I realized it was a different answer. – Mike Wise Mar 15 '17 at 10:26