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
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)
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
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
