attach(olives)
plot(palmitoleic, oleic, xlim=c(-10,300))
plot.contours2d(palmitoleic, oleic)

# alternatively use library ellipse

library(ellipse)
plot(palmitoleic, oleic, xlim=c(-10,300))
for( l in c(0.5, 0.8, 0.95, 0.99))
   lines(ellipse(cor(cbind(palmitoleic, oleic)), 
                 scale=c(sd(palmitoleic), sd(oleic)),
                 centre=c(mean(palmitoleic), mean(oleic)),
                 level=l))
                 

plot.contours2d <- function(x, y, clevels=c(0.5, 0.8, 0.95, 0.99), label=T, vc=NULL, col="blue") {
   # a crude (half) circle  
   n   <- 100 
   cx  <-seq(-1,1,2/n)[-(trunc(n/2)+1)]
   cy  <-sqrt(1-cx**2)
   
   # replicate circles and add NAs to avoid connection of ellipses
   tx  <-rep(c(cx,  rev(cx), NA), length(clevels))
   ty  <-rep(c(cy, rev(-cy), NA), length(clevels))
   cxy <-cbind(tx, ty)

   rclevels <- rep(rep(clevels, rep(2*n+1, length(clevels))), 2)
         
   # scale the circles to be confidence circles
   cxys <- cxy * sqrt(qchisq(matrix(rclevels, ncol=2, byrow=F),2))
      
   # estimate VC-matrix and generate \Sigma^{-1/2}
   if (is.null(vc))
      vc <- matrix(c(var(x), cov(x, y), cov(x, y), var(y)),ncol=2,byrow=T)
   evc<-eigen(vc)
   vc2<-(evc$vectors)%*%diag(sqrt(evc$values))%*%t(evc$vectors)
  
   # multiply circles to rotate and scale the ellipses 
   scxys <- cxys%*%vc2
   
   # shift means
   scxys[,1] <- mean(x)+scxys[,1]
   scxys[,2] <- mean(y)+scxys[,2]
   
   lines(scxys, col=col)
   if( label ) {
      ind <- (0:(length(clevels)-1))*(2*n+1)+1
      text(scxys[ind,1], scxys[ind,2], rclevels[ind], col="red")
   }
}