3

I'm hoping to be able to take a 3x3 covariance matrix and turn this into an error ellipsoid but so far I haven't been able to achieve this.

I'm very new to R (in fact turned to it to attempt to solve this problem) - and after finding a working example for a 2x2 matrix I've not managed to modify this for 3x3 successfully - can anyone point me in the right direction as to which changes are required please?

The working code for the 2x2 matrix (found here in another question) is below:

ctr    <- c(0, 0)                               
A      <- matrix(c(2.2, 0.4, 0.4, 2.8), nrow=2) 
RR     <- chol(A)                               
angles <- seq(0, 2*pi, length.out=200)          
ell    <- 1 * cbind(cos(angles), sin(angles)) %*% RR  
ellCtr <- sweep(ell, 2, ctr, "+")               

eigVal  <- eigen(A)$values
eigVec  <- eigen(A)$vectors
eigScl  <- eigVec %*% diag(sqrt(eigVal))  
xMat    <- rbind(ctr[1] + eigScl[1, ], ctr[1] - eigScl[1, ])
yMat    <- rbind(ctr[2] + eigScl[2, ], ctr[2] - eigScl[2, ])
ellBase <- cbind(sqrt(eigVal[1])*cos(angles), sqrt(eigVal[2])*sin(angles)) 
ellRot  <- eigVec %*% t(ellBase)                                          
plot((ellRot+ctr)[1, ], (ellRot+ctr)[2, ], asp=1, type="l", lwd=2)
matlines(xMat, yMat, lty=1, lwd=2, col="green")
points(ctr[1], ctr[2], pch=4, col="red", lwd=3)

Many thanks in advance for any help.

anthr
  • 937
  • 1
    You can install packages car and rgl and use function car:::ellipsoid (car doesn't export this function, so you need the :::). It returns an rgl mesh object which can be displayed with shade3d(object, options) and wire3d(object, options). A use case is in car:::scatter3d.default - look for the lines following if (ellipsoid) { and run example(scatter3d) after loading car. – caracal Mar 25 '14 at 08:42
  • 1
    The following threads (most of which have good answers) all seem to contain relevant information: http://stats.stackexchange.com/questions/89512, http://stats.stackexchange.com/questions/67422, http://stats.stackexchange.com/questions/69350, and http://stats.stackexchange.com/questions/29860. – whuber Mar 25 '14 at 16:10
  • Thanks caracal for the suggestions - I'll take a look at the details of this function and see what I can do with it. Thanks too to whuber for the thread suggestions - will read through these. – anthr Mar 25 '14 at 19:58

1 Answers1

3

I'm interested in this too, so I studied a bit and came up with an incomplete answer. This will produce an ellipsoid and mark its radii, but I haven't figured out how to rotate it using the eigenvectors.

Loading the rgl package and running demo(shapes3d) will give you this ellipsoid3d function:

ellipsoid3d <- function(rx=1,ry=1,rz=1,n=30,ctr=c(0,0,0),
                         qmesh=FALSE,
                         trans = par3d("userMatrix"),...) {
   if (missing(trans) && !rgl.cur()) trans <- diag(4)
   degvec <- seq(0,pi,length=n)
   ecoord2 <- function(p) {
     c(rx*cos(p[1])*sin(p[2]),ry*sin(p[1])*sin(p[2]),rz*cos(p[2])) }
   v <- apply(expand.grid(2*degvec,degvec),1,ecoord2)
   if (qmesh) v <- rbind(v,rep(1,ncol(v))) ## homogeneous
   e <- expand.grid(1:(n-1),1:n)
   i1 <- apply(e,1,function(z)z[1]+n*(z[2]-1))
       i2 <- i1+1
   i3 <- (i1+n-1) %% n^2 + 1
   i4 <- (i2+n-1) %% n^2 + 1
   i <- rbind(i1,i2,i4,i3)
   if (!qmesh)
     quads3d(v[1,i],v[2,i],v[3,i],...)
   else return(rotate3d(qmesh3d(v,i,material=list(...)),matrix=trans))
 }

For demonstration, I'll use this covariance matrix: $$\left(\begin{array} {cc}2.2&.4&.2\\.4&2.8&.3\\.2&.3&2.4\end{array}\right)$$

Covariance.Matrix=matrix(c(2.2,.4,.2,.4,2.8,.3,.2,.3,2.4),3)
Eigenvalue.Roots=eigen(Covariance.Matrix)$values^.5->X

require(rgl)
wire3d(ellipsoid3d(X[1],X[2],X[3],qmesh=T),col='lightgreen')
axes3d("x",pos=c(0,0,0),at=c(0,X[1]),labels='',lwd=2,col='maroon')
axes3d("y",pos=c(0,0,0),at=c(0,X[2]),labels='',lwd=2,col='maroon')
axes3d("z",pos=c(0,0,0),at=c(0,X[3]),labels='',lwd=2,col='maroon')

Rotating the ellipsoid3d object is possible by using rotate3d() inside wire3d(), but I'm not sure how to find the angles $\phi,\theta,\psi$ from the eigenvectors. This is my best guess so far (based on this):

Phi=with(eigen(Covariance.Matrix),atan(vectors[3,2]/vectors[3,1]))
Theta=acos(eigen(Covariance.Matrix)$vectors[3,3])
Psi=with(eigen(Covariance.Matrix),atan(vectors[2,3]/-vectors[1,3]))

If these are known, the euler function given here may help produce a matrix for rotate3d's matrix argument...but I'm not sure it works. With sample data, I seem to get a different rotation using John Fox's scatter3d function. If you have a dataset of observation trios, not just a covariance matrix, scatter3d seems to work well.

Nick Stauner
  • 12,342
  • 5
  • 52
  • 110
  • 1
    Many thanks for your comment - even if incomplete this is clearly useful. I was also unsure how to handle the rotations - maybe you've seen it already but these are demonstrated in a 2x2 matrix in the thread I linked too in the original post, although in a different example that I quoted. I will take a look at implementing what you suggest today and seeing if I can resolve this remaining question. Thanks again. – anthr Mar 25 '14 at 20:02
  • Yeah, rotation is much simpler in 2D. In 3D, it seems $\theta$ is an angle of rotation around the new $y$ axis following $\phi$ rotation around the original $z$ axis. And then there's $\psi$, with which I've never dealt before today! – Nick Stauner Mar 25 '14 at 20:05
  • In case it is useful to you - I found some examples in Mathematica here. – anthr Mar 26 '14 at 03:36