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

From: John Fox <jfox_at_mcmaster.ca>
Date: Fri 03 Feb 2006 - 12:56:22 EST


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@stat.math.ethz.ch
> [mailto:r-help-bounces@stat.math.ethz.ch] On Behalf Of Duncan Murdoch
> Sent: Thursday, February 02, 2006 6:22 PM
> To: Michael
> Cc: r-help@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@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@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@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@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@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 Received on Fri Feb 03 13:04:33 2006

This archive was generated by hypermail 2.1.8 : Fri 03 Feb 2006 - 16:19:27 EST