[R] Observation in a confidence ellipse

Rolf Turner rolf.turner at xtra.co.nz
Sun May 29 01:28:42 CEST 2011


On 29/05/11 05:45, Jessica Minkue wrote:
> Hello everyone
> I really need some help here. I made a confidence ellipse using the function ellipse from the package ellipse:
>
> ellipse(SD, centre=colMeans(pcsref),t=sqrt((p * (n-1)/(n-p))*qf(0.99, p,n-p))
>
>
> Now, I want to write a function whom return TRUE or FALSE if a given observation is in the confidence ellipse. But I have no clue how to do it
> Can anyone help me?

First of all, you are probably way off base talking about 
***confidence*** ellipses here.
If you are testing whether observations are inside the ellipses, then 
you are most
likely interested in ***prediction*** ellipses.

It is vital that you understand the difference.

But to answer your question:  It would be easy enough to do it from scratch.
Let your ellipse be defined by

     (x - mu)' M(x-mu) = c

where mu is the ``centre'', x is a 2-vector (x_1,x_2)', M is a positive 
definite matrix,
c is a positive constant, and " ' " denotes ``transpose''.  Then a point 
x = (x_1,x_2)' is
inside the ellipse if and only if (x - mu)' M(x-mu) <= c.

Coding this up is an easy exercise.  If you can't do it, you probably 
shouldn't be messing
with this stuff in the first place.

However if you want to use a sledge-hammer to crack a peanut,
install the "spatstat" package and then do:

require(spatstat)
W <- owin(poly=<your ellipse>)
inside.owin(x1,x2,W) # Where x1 and x2 are the x and y coordinates of 
the points you wish to test.

     cheers,

         Rolf Turner



More information about the R-help mailing list