2
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – caracal Mar 25 '14 at 8:42
  • 1
    $\begingroup$ The following threads (most of which have good answers) all seem to contain relevant information: stats.stackexchange.com/questions/89512, stats.stackexchange.com/questions/67422, stats.stackexchange.com/questions/69350, and stats.stackexchange.com/questions/29860. $\endgroup$ – whuber Mar 25 '14 at 16:10
  • $\begingroup$ 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. $\endgroup$ – anthr Mar 25 '14 at 19:58
3
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – anthr Mar 25 '14 at 20:02
  • $\begingroup$ 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! $\endgroup$ – Nick Stauner Mar 25 '14 at 20:05
  • $\begingroup$ In case it is useful to you - I found some examples in Mathematica here. $\endgroup$ – anthr Mar 26 '14 at 3:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.