How to plot bivariate normal distribution with expanding ellipses

Question

How to plot bivariate normal distribution with expanding ellipses and add 5%, 25%, 50%, 75% and 95% label in the plot? Thank you!

Answer 1

You can create a contour plot using an R package called mvtnorm.

Let's say you're trying to plot a bivariate normal distribution where mu_x = 1 and mu_y = 1 and variance matrix is c(2,1,1,1). Generate 100 observations for x,y,z. You can create a contour plot for this scenario as such:

library(mvtnorm)
x.points <- seq(-3,3,length.out=100)
y.points <- x.points
z <- matrix(0,nrow=100,ncol=100)
mu <- c(1,1)
sigma <- matrix(c(2,1,1,1),nrow=2)
for (i in 1:100) {
   for (j in 1:100) {
    z[i,j] <- dmvnorm(c(x.points[i],y.points[j]),
    mean=mu,sigma=sigma)
   }
}
contour(x.points,y.points,z)

Answer 2

Here is a solution that computes the contours at the levels that you want

#####
# Compute points and rotation matrix

# input
theta <- c(1, 2)
sigma <- diag(c(3^2, 2^2))
sigma[2, 1] <- sigma[1, 2] <- sqrt(sigma[1, 1] * sigma[2, 2]) * .5 

# we start from points on the unit circle
n_points <- 100
xy <- cbind(sin(seq(0, 2 * pi, length.out = n_points)), 
            cos(seq(0, 2 * pi, length.out = n_points)))

# then we scale the dimensions
ev <- eigen(sigma)
xy[, 1] <- xy[, 1] * 1 
xy[, 2] <- xy[, 2] * sqrt(min(ev$values) / max(ev$values))

# then rotate
phi <- atan(ev$vectors[2, 1] / ev$vectors[1, 1])
R <- matrix(c(cos(phi), sin(phi), -sin(phi), cos(phi)), 2) 
xy <- tcrossprod(R, xy)

# the quantiles you ask for
chi_vals <- qchisq(c(.05, .25, .50, .75, .95), df = 2) * max(ev$values)

#####
# Plot contours
par(mar = c(4.5, 4, .5, .5))
plot(c(-8, 10), c(-4, 8), type = "n", xlab = "x", ylab = "y")
for(r in sqrt(chi_vals))
  lines(r * xy[1, ] + theta[1], r * xy[2, ] + theta[2], lty = 1)

Detailed Explanation

theta <- c(1, 2)
sigma <- diag(c(3^2, 2^2))
sigma[2, 1] <- sigma[1, 2] <- sqrt(sigma[1, 1] * sigma[2, 2]) * .5 

# we start from points on the unit circle
n_points <- 100
xy <- cbind(sin(seq(0, 2 * pi, length.out = n_points)), 
            cos(seq(0, 2 * pi, length.out = n_points)))
par(mar = c(5, 5, 3, 1), mfcol = c(2, 2))
plot(xy[, 1], xy[, 2], xlab = "x", ylab = "y", xlim = c(-8, 10), bty = "l",
     ylim = c(-4, 8), main = "Unit circle", type = "l")
arrows(c(0, 0), c(0, 0), c(0, 1), c(1, 0), length = .05)

# this is very much like PCA. We scale the dimensions such that the first 
# dimension has length one and the others are scaled proportional to the square
# root of their eigenvalue relative to the largest eigenvalue
ev <- eigen(sigma)
scal <- sqrt(min(ev$values) / max(ev$values))
xy[, 2] <- xy[, 2] * scal
plot(xy[, 1], xy[, 2], xlab = "x", ylab = "y", xlim = c(-8, 10), bty = "l",
     ylim = c(-4, 8), main = "Scaled", type = "l")
arrows(c(0, 0), c(0, 0), c(0, 1), c(scal, 0), 
       length = .05)

# then we rotate phi degrees to account for the correlation coefficient. See 
#   https://en.wikipedia.org/wiki/Rotation_matrix
# and notice that we compute the angle of the first eigenvector  
phi <- atan(ev$vectors[2, 1] / ev$vectors[1, 1])
R <- matrix(c(cos(phi), sin(phi), -sin(phi), cos(phi)), 2) 
xy <- tcrossprod(R, xy) # R %*% t(xy)
plot(xy[1, ], xy[2, ], xlab = "x", ylab = "y", xlim = c(-8, 10), bty = "l",
     ylim = c(-4, 8), main = "Rotated", type = "l")
arrs <- tcrossprod(R, matrix(c(0, 1, scal, 0), 2L))
arrows(c(0, 0), c(0, 0), arrs[1, ], arrs[2, ], length = .05)

# the right size of each circle can now be found by taking the wanted 
# quantile from the chi-square distribution with two degrees of freedom 
# multiplied by the largest eigenvalue
plot(c(-8, 10), c(-4, 8), type = "n", xlab = "x", ylab = "y", main = "Final", 
     bty = "l")
chi_vals <- qchisq(c(.05, .25, .50, .75, .95), df = 2) * max(ev$values)
for(r in sqrt(chi_vals))
  lines(r * xy[1, ] + theta[1], r * xy[2, ] + theta[2], lty = 1)