[R] density ellipses?

[Hans Ehrbar]

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