[R] density ellipses?

[Thaddeus Tarpey]

Dear Kaspar,

Below is some code which will plot an ellipse.  This code was written by
Bernard Flury and Marco Bee originally in S+ for Flury's book "A First
Course in Multivariate Statistics" (Springer 1997).  You would have to
plug in the mean and covariance matrix for m and A respectively for a
density ellipse.

#
#  Plot points satisfying this equation of an ellipse:
#                   -1             2
#           (x-m)' A   (x-m)  <-  c  .
#
#*****************************  CUT HERE  ********************************


# First define all parameters 

A <- matrix(c(1,1,1,2), ncol = 2)		 # define matrix A
m <- c(3, 4)					 # define vector m
const <- 1			 # define constant 
k <- 1000		  # define number of points on the ellipse

# procedure ELLIPS
# procedure ELLIPS

ellips <- function(A, m, const, k)
{
# input:  A	positive definite symmetric matrix of dimension 2 by 2
#         m	column vector of length 2, center of ellipse
#         const	positive constant
#         k	number of points on ellipse (must be at least 2)
# output: x 	a (k by 2) matrix whose rows are the coordinates
#		of k points on the ellipse (y-m)'*A^{-1}*(y-m) = c^2

r <- A[1, 2]/sqrt(A[1, 1] * A[2, 2])
Q <- matrix(0, 2, 2)			  # construct matrix Q for
Q[1, 1] <- sqrt(A[1, 1] %*% (1+r)/2)	# transformation of circle
Q[1, 2] <- -sqrt(A[1, 1] %*% (1-r)/2)		      # to ellipse
Q[2, 1] <- sqrt(A[2, 2] %*% (1+r)/2)
Q[1, 1] <- sqrt(A[2, 2] %*% (1-r)/2)
alpha <- seq(0, by = (2 * pi)/k, length = k)	   # define angles
Z <- cbind(cos(alpha), sin(alpha)) 	   # points on unit circle
X <- t(m + const * Q %*% t(Z))  # coordinates of points on ellipse
X
}					 # end of procedure ellips

# call procedure ellips

X <- ellips(A, m, const, k)		   

# graph the results

win.graph()
plot(X[,1], X[,2])

#***********************************************************

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
> 
> -- 
> 
> Kaspar Pflugshaupt
> Geobotanisches Institut
> Zuerichbergstr. 38
> CH-8044 Zuerich
> 
> Tel. ++41 1 632 43 19
> Fax  ++41 1 632 12 15
> 
> mailto:pflugshaupt at geobot.umnw.ethz.ch
> privat:pflugshaupt at mails.ch
> http://www.geobot.umnw.ethz.ch
> 
> -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
> 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
> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
> 


*****************************************************************************
Thaddeus Tarpey                                   Phone: (937) 775-2861
Wright State University                             Fax: (937) 775-2081
Department of Mathematics and Statistics
Dayton,  Ohio  45435
Home-page: http://www.math.wright.edu/People/Thad_Tarpey/thad.htm
*****************************************************************************

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._