[R] density ellipses?
Thaddeus Tarpey
thaddeus.tarpey at wright.edu
Wed Mar 22 19:23:10 CET 2000
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list