rgl: ellipse3d with axes

Thanks Duncan (& others)

Here is a function that does what I want in this case, and tries to do 
it to work generally
with ellipse3d.  (Note that I reverse the order of centre and scale 
'cause I was bitten
by trying ellipse3d.axes(cov, mu))

# draw axes in the data ellipse computed by ellipse3d
ellipse3d.axes <-
function (x, centre = c(0, 0, 0), scale = c(1, 1, 1), level = 0.95,
    t = sqrt(qchisq(level, 3)), which = 1:3, ...)
{
    stopifnot(is.matrix(x)) # should test for square, symmetric
    cov <- x[which, which]
    eigen <- eigen(cov)
    # coordinate axes, (-1, 1), in pairs
    axes <- matrix(
      c(0, 0, -1,   0, 0, 1,
        0, -1, 0,   0, 1, 0,
       -1, 0, 0,    1, 0, 0),  6, 3, byrow=TRUE)

    # transform to PC axes
    axes <- axes %*% sqrt(diag(eigen$values)) %*% t(eigen$vectors)
    result <- scale3d(axes, t, t, t)
    if (!missing(scale))
        if (length(scale) != 3) scale <- rep(scale, length.out=3)
        result <- scale3d(result, scale[1], scale[2], scale[3])
    if (!missing(centre))
        if (length(centre) != 3) scale <- rep(centre, length.out=3)
        result <- translate3d(result, centre[1], centre[2], centre[3])
    segments3d(result, ...)
    invisible(result)
}

Test case:

library(rgl)
data(iris)
iris3 <- iris[,1:3]
cov <- cov(iris3)
mu <- mean(iris3)
col <-c("blue", "green", "red")[iris$Species]
plot3d(iris3, type="s", size=0.5, col=col, cex=2, box=FALSE)
plot3d( ellipse3d(cov, centre=mu, level=0.68), col="gray", alpha=0.2,  
add = TRUE)

axes <- ellipse3d.axes(cov, centre=mu)

One thing I can't explain, compared to your example is why the my axes 
extend outside the ellipse,
whereas yours didn't.

One final remark-  I knew that axes %*% chol(cov) did not give the 
orthogonal PC axes I wanted,
but at least it gave me something on the right scale and location. But 
these axes also turn out to be
useful for visualizing multivariate scatter and statistical concepts.  
chol() gives the factorization of
cov that corresponds to the Gram-Schmidt orthogonalization of a data 
matrix -- orthogonal axes
in the order x1,  x2|x1,  x3|x1, x2, ..., and vector length and 
orientation in this coordinate system
correspond to Type I SS in linear models.
Thus, I could see generalizing my ellipse3d.axes function further to 
allow a type=c("pca", "chol")
argument.

-Michael

That's easy, but it doesn't give you the principal axes of the 
ellipse.  Just use

axes %*% chol(cov)

If you start with a unit sphere, this will give you points on its 
surface, but not the ones you want.  For those you need the SVD or 
eigenvectors.  This looks like it does what you want:

axes <- matrix(
    c(0, 0, 0, # added origin
       0, 0, -1,   0, 0, 1,
       0, -1, 0,   0, 1,  0,
       -1, 0, 0,   1, 0, 0),  7, 3, byrow=TRUE)
axes <- axes[c(1,2,1,3,1,4,1,5,1,6,1,7),]  # add the origin before each

cov <- cov(trees)
eigen <- eigen(cov)
shade3d(ellipse3d(cov, t=1, alpha=0.2, col='red'))
segments3d(axes %*% sqrt(diag(eigen$values)) %*% t(eigen$vectors))

Duncan Murdoch

Michael Friendly     Email: friendly AT yorku DOT ca 
Professor, Psychology Dept.
York University      Voice: 416 736-5115 x66249 Fax: 416 736-5814
4700 Keele Street    http://www.math.yorku.ca/SCS/friendly.html
Toronto, ONT  M3J 1P3 CANADA