[R] density ellipses?

About this list Date view Thread view Subject view Author view Other groups

Subject: [R] density ellipses?
From: Hans Ehrbar (ehrbar@econ.utah.edu)
Date: Thu 23 Mar 2000 - 06:11:13 EST 

Message-Id: <200003222011.NAA14208@marx.econ.utah.edu>

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@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@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@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

About this list Date view Thread view Subject view Author view Other groups

This archive was generated by hypermail 2b25 : Thu 23 Mar 2000 - 20:01:16 EST