[R] density ellipses?

[Ben Bolker]

  Here are two ellipse() functions I've written in the past, probably very
inefficient but they might give you something to start with.  One takes
the center, semiminor and major axes, and rotation, the other takes a
variance-covariance matrix (it was designed for confidence ellipses using
the likelihood ratio test).
   Hope they give you a start.

ellipse _ function(x=1,y=1,r1=1,r2=1,ang=0,narc=200,...)
{
  rmat <- c(cos(ang),-sin(ang),sin(ang),cos(ang))
  calcpt _ function(x,y,r1,r2,ang1) {
    r _ rmat %*% c(r1*cos(ang1),r2*sin(ang1))
    c(x+r[1],y+r[2])
  }
  arcs _ matrix(sapply(seq(0,2*pi,len=narc),
                       function(q)calcpt(x,y,r1,r2,q)),byrow=TRUE,ncol=2)
  lines(arcs,...)
}

vcellipse <- function(ctr,varcov,correlation=NULL,
                      critval=qchisq(0.95,2),nx=200,
                    fill=FALSE,...) {
  # draw the ellipse defined by the center ctr, and either a
  # variance vector c(sigma^2(x),sigma^2(y)) and a correlation,
  # OR a variance-covariance matrix, plus an optional critical value.
  # nx is the number of separate arcs used to draw the ellipse
  # plot solution for y of equation: h11*x^2+2*h12*x*y+h22*y^2=critval
  #
  # reconstruct hessian matrix
  if (!is.null(correlation) && (correlation>1 || correlation<(-1)))
    stop(paste("Correlation=",correlation," is out of bounds: parameter mixup?"))
  if (!is.null(correlation) && is.vector(varcov)) {
    v <- matrix(nrow=2,ncol=2)
    v[1,1] <- varcov[1]
    v[2,2] <- varcov[2]
    v[1,2] <- correlation*sqrt(v[1,1]*v[2,2])
    v[2,1] <- v[1,2]
  } else if (is.matrix(varcov)) {
    v <- varcov
  }
  h <- solve(v)
  eps <- 1e-7
  ##  xrad  <- sqrt(solve(hess)[1,1]*critval)-eps
  xrad <- sqrt(-h[2,2]*critval/(h[1,2]^2-h[1,1]*h[2,2]))-eps  # x range
  ##  segments(ctr[1],ctr[2],ctr[1]+xrad,ctr[2],...)
  ## det >= 0
  ## (h12^2-h11*h22)*x1^2 > -h22*critval
  ## x1^2 > -h22*critval/(h12^2-h11*h22)
  ## x1   > sqrt(-h22*critval/(h12^2-h11*h22)),"\n")
  ## same as inverting the matrix!
  x1 <- seq(-xrad,xrad,length=round(nx/2))
  det <- sqrt((h[1,2]^2-h[1,1]*h[2,2])*x1^2+h[2,2]*critval)/h[2,2]
  # calculate "top" and "bottom" halves (plus/minus solutions to quadratic)
  y1 <- -h[1,2]/h[2,2]*x1 + det 
  x2 <- rev(x1)
  y2 <- -h[1,2]/h[2,2]*x2 - det
  # translate ellipse to center
  x <- c(x1,x2)+ctr[1]
  y <- c(y1,y2)+ctr[2]
  x <- c(x,x[1])   # complete ellipse
  y <- c(y,y[1])
  if (fill)
    polygon(x,y,...)
  else
    lines(x,y,...)
}

On Wed, 22 Mar 2000, Kaspar Pflugshaupt wrote:

> Hello,
> 
> has anybody written a function to plot density ellipses (95%, 99% or
> anything) in a scatterplot? I found nothing in any package, nor in the list
> archives. 
> 
> There does seem to be a contributed package "ellipse" for S-Plus (on
> S-Archive), but it does a lot more than what I would need. Still, if anybody
> ported it to R, I'd be grateful for a link. I'm a bit afraid to try the port
> myself, since the routines do some magic with postscript preambles that I
> don't understand.
> 
> Or is there a simple way to implement it myself? The axis directions could
> be extracted from princomp, but how would one plot the ellipses?
> 
> Thanks for any tips
> 
> Kaspar
> 
> 

-- 
Ben Bolker                                  bolker at zoo.ufl.edu
Zoology Department, University of Florida   http://www.zoo.ufl.edu/bolker
318 Carr Hall/Box 118525                    tel: (352) 392-5697
Gainesville, FL 32611-8525                  fax: (352) 392-3704

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._