# [R] rgl: ellipse3d with axes

Michael Friendly friendly at yorku.ca
Wed Sep 24 18:32:54 CEST 2008

```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, scale, scale)
if (!missing(centre))
if (length(centre) != 3) scale <- rep(centre, length.out=3)
result <- translate3d(result, centre, centre, centre)
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,

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

Duncan Murdoch wrote:
> 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)
> 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