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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
On 22-Mar-00 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 > > KasparYou can calculate a 2D kernel density estimate on a regular grid using the bkde2D() function (package Kernsmooth) and then plot highest density regions with contour(). An example giving a 90% region: library(KernSmooth) x <- rnorm(1000) y <- (x+rnorm(1000))/sqrt(2) z <- bkde2D(cbind(x,y), bandwidth=c(0.4,0.4)) zsort <- sort(z$fhat) p <- cumsum(zsort)/sum(zsort) contour(z$x1,z$x2,z$fhat, levels=min(zsort[p>=0.10])) points(x,y) You need to play with the bandwidth. Martyn -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Kaspar Pflugshaupt asked:> 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.The following very elegant function computes confidence ellipses using the F-statistic. It was sent to S-news last year by Roger Koenker (I did some small editing). I looked at the math behind it in my econometrics class, if you go to my class notes http://www.econ.utah.edu/ehrbar/ecmetrcs.pdf (4.6 megabytes) and select Multivariate Normal Density -> Special Case: Bivariate -> Level lines you will find a proof why this works. confelli <- function(b, C, df, level = 0.95, xlab = "", ylab = "", add=T, prec=51) # Plot an ellipse with "covariance matrix" C, center b, and P-content # level according the F(2,df) distribution. # Sent to S-NEWS on May 19, 1999 by Roger Koenker # Department of Economics # University of Illinois # Champaign, IL 61820 # url: http://www.econ.uiuc.edu # email roger at ysidro.econ.uiuc.edu # vox: 217-333-4558 # fax: 217-244-6678. { d <- sqrt(diag(C)) dfvec <- c(2, df) phase <- acos(C[1, 2]/(d[1] * d[2])) angles <- seq( - (PI), PI, len = prec) mult <- sqrt(dfvec[1] * qf(level, dfvec[1], dfvec[2])) xpts <- b[1] + d[1] * mult * cos(angles) ypts <- b[2] + d[2] * mult * cos(angles + phase) if(add) lines(xpts, ypts) else plot(xpts, ypts, type = "l", xlab = xlab, ylab = ylab) } -- Hans G. Ehrbar ehrbar at econ.utah.edu Economics Department, University of Utah (801) 581 7797 (my office) 1645 Campus Center Dr., Rm 308 (801) 581 7481 (econ office) Salt Lake City UT 84112-9300 (801) 585 5649 (FAX) -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._