[R] is there a way to visualize 3D normal distributions?

[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 
-------------------------------- 

> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch 
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Duncan Murdoch
> Sent: Thursday, February 02, 2006 6:22 PM
> To: Michael
> Cc: r-help at stat.math.ethz.ch; Ben Bolker
> Subject: Re: [R] is there a way to visualize 3D normal distributions?
> 
> On 2/2/2006 5:01 PM, Michael wrote:
> > shape3d only gives rigid sphere... not the free form 
> ellipsoid that I 
> > want...
> 
> ? 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
> 
> 
> > 
> > On 2/2/06, Ben Bolker <bolker at ufl.edu> wrote:
> >> Duncan Murdoch <murdoch <at> stats.uwo.ca> writes:
> >>
> >>> On 2/2/2006 3:39 AM, Michael wrote:
> >>>> Hi all,
> >>>>
> >>>> How do I visualize a contour of a tri-variate normal 
> distribution?
> >>>>
> >>>> I just like to see the ellipsoid very much. I hope there 
> is a easy 
> >>>> way
> >> or
> >>>> existing method in R.
> >>> The misc3d package includes a function for 3d contour plots; that 
> >>> should do what you want.
> >>>
> >>   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
> >>
> >> ______________________________________________
> >> R-help at stat.math.ethz.ch mailing list 
> >> https://stat.ethz.ch/mailman/listinfo/r-help
> >> PLEASE do read the posting guide!
> >> http://www.R-project.org/posting-guide.html
> >>
> > 
> > 	[[alternative HTML version deleted]]
> > 
> > ______________________________________________
> > R-help at stat.math.ethz.ch mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide! 
> > http://www.R-project.org/posting-guide.html
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html