On 24/09/2008 12:32 PM, Michael Friendly wrote:
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.
That's just because you specified level in the ellipse3d call, but not
in the ellipes3d.axes call.
One thing that looks a little strange is that the PC axes don't appear
to be orthogonal: this is because the scaling is not the same on all
coordinates. It might look better doing the first plot as
plot3d(iris3, type="s", size=0.5, col=col, cex=2, box=FALSE, aspect="iso")
Duncan Murdoch
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
______________________________________________
[email protected] 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.
