## October 4, Ellipses 
rm(list = ls())

Mu <- matrix(c(5,2),2,1)
Sigma <- matrix(c(1, 0.7, 0.7, 1), 2,2, byrow=T)

alpha <- 0.05  # thus 1-alpha probability contour
prob <- 1-alpha
c.sq <- qchisq(prob,2)
eig.val <- eigen(Sigma)$values
eig.vec <- eigen(Sigma)$vectors

p.cen <- t(Mu)
p.1 <- sqrt(c.sq)*sqrt(eig.val[1])*eig.vec[,1]
p.2 <- sqrt(c.sq)*sqrt(eig.val[2])*eig.vec[,2]

Ell <- rbind(p.cen, p.cen+p.1, p.cen+p.2, p.cen+(-p.1), p.cen+(-p.2))

library(mvtnorm)
sam <- rmvnorm(500, Mu, Sigma)
plot(sam[,1], sam[,2], pch=16, xlab="x1", ylab="x2", cex.lab=1.5)
points(Ell[,1], Ell[,2], pch=16, cex=3, col="blue")

library(ellipse)
plot(sam[,1], sam[,2], pch=16, xlab="x1", ylab="x2", cex.lab=1.5)
points(Ell[,1], Ell[,2], pch=16, cex=3, col="blue")
lines(ellipse(Sigma, centre = c(Mu), level=prob), lwd=3)


## Add in another ellipse 
prob.2 <- 0.50
lines(ellipse(Sigma, centre = c(Mu), level=prob.2), lwd=3, col="red")