Skip to content

is there a way to visualize 3D normal distributions?

10 messages · Duncan Murdoch, John Fox, Michael +1 more

Duncan Murdoch
#
On 2/2/2006 3:39 AM, Michael wrote:
The misc3d package includes a function for 3d contour plots; that should 
do what you want.

There are some errors in the rendering when plotting multiple contours, 
due to limitations in rgl (the underlying package that does the 
rendering).  The problem is that transparent surfaces need to be 
rendered in order from back to front, and rgl doesn't do that.  But it 
still looks good.

Duncan Murdoch
Ben Bolker
# 
Duncan Murdoch <murdoch <at> stats.uwo.ca> writes:
is contour3d really necessary or could you just plot ellipsoids?
(library(rgl); demo(shapes3d)) -- still a little bit of figuring 
to do, but this should get you most of the way there.

  Ben Bolker
John Fox
# 
Dear Michael,

Some time ago, I posted to r-help a solution to a problem very close to
this, which adapts the scatter3d() function in the Rcmdr package so that it
plots concentration ellipsoids. You'll find the code at
<http://finzi.psych.upenn.edu/R/Rhelp02a/archive/61156.html>.

I hope this helps,
 John

--------------------------------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
--------------------------------
Ben Bolker
#
Michael wrote:
one or the other of us is missing something.

    after running demo(shapes3d)
[to define ellipsoid3d] and rgl.clear(),

s1 <- ellipsoid3d(qmesh = TRUE, trans = diag(4))
shade3d(rotate3d(scale3d(s1, 1, 1, 2),
                  angle=pi/4, x=0, y=1, z=0),
         col = "red")
axis3d()

  produces an ellipsoid with radii (1,1,2) rotated by pi/4 around
the y axis.  I'm too lazy to figure out how to
translate a variance-covariance matrix into
a scaling/rotation matrix (the scale will want
to be something like eigen(v)$values*qnorm(0.975)
for a 95% contour, the rotation matrix is some
function of the eigenvectors), but it seems to
me that this will produce any ("free form"?) ellipsoid
you want.

   cheers
     Ben Bolker
Duncan Murdoch
#
On 2/2/2006 5:01 PM, Michael wrote:
? If you run the demo, you'll see ellipsoids...

You just need to work out the appropriate transform to apply to a sphere 
to get the ellipsoid you want.  I imagine something like this:

sphere <- ellipsoid3d(2,2,2, qmesh=TRUE)
ellipsoid <- translate3d(rotate3d(sphere, matrix=chol(S)), xbar, ybar, zbar)

shade3d(ellipsoid)

is what you want, where S is the covariance matrix, and xbar,ybar,zbar 
have the obvious meaning.  (rotate3d() is used with a matrix that isn't 
a rotation matrix; it may not be obvious that this is allowed, but it is.)

Duncan Murdoch
John Fox
# 
Dear Duncan, Michael, and Ben,

The code to which I referred Michael works by deforming a sphere. Here's the
central function:

ellipsoid <- function(center=c(0, 0, 0), radius=1, shape=diag(3), n=30){
# adapted from the shapes3d demo in the rgl package
  degvec <- seq(0, 2*pi, length=n)
  ecoord2 <- function(p) c(cos(p[1])*sin(p[2]), sin(p[1])*sin(p[2]),
cos(p[2]))
  v <- t(apply(expand.grid(degvec,degvec), 1, ecoord2))
  v <- center + radius * t(v %*% chol(shape))
  v <- rbind(v, rep(1,ncol(v))) 
  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)
  qmesh3d(v, i)
  }

Here are the key lines from scatter3d() for the shape and radius (slightly
edited to remove the context):

	dfn <- 3
	dfd <- length(x) - 1
      radius <- sqrt(dfn * qf(level, dfn, dfd))
      ellips <- ellipsoid(center=c(mean(x), mean(y), mean(z)), 
            shape=cov(cbind(x,y,z)), radius=radius)
      shade3d(ellips, col=surface.col[1], alpha=0.1, lit=FALSE)
      wire3d(ellips, col=surface.col[1], lit=FALSE)

Regards,
 John

--------------------------------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
--------------------------------
Ben Bolker
# 
In response to your last question --

   Duncan already gave you the answer in his last e-mail -- to give
an explicit example, for a particular set of means and
variance-covariance matrix:

library(rgl)
demo(shapes3d)
rgl.clear()
sphere <- ellipsoid3d(2,2,2,qmesh=TRUE)
means <- c(0,0,0)
S <- matrix(c(3,1,-1,1,2,-1,-1,-1,2),nrow=3)
ellipsoid <- translate3d(rotate3d(sphere,matrix=chol(S)*qnorm(0.975)),
                                     means[1],means[2],means[3])
shade3d(ellipsoid)
axis3d()
Duncan Murdoch wrote: