Purpose
To Draw the Location Dispersion Ellipsoid

> library(ellipse)
> library(mnormt)
> n <- 10000
> sample.mean <- c(0, 0)
> sample.cov <- matrix(data = c(1, 0.6, 0.6, 1), nrow = 2, ncol = 2)
> x1 <- rmnorm(n = n, mean = sample.mean, varcov = sample.cov)
> sample.cov <- matrix(data = c(1, -0.6, -0.6, 1), nrow = 2, ncol = 2)
> x2 <- rmnorm(n = n, mean = sample.mean, varcov = sample.cov)
> sample.cov <- matrix(data = c(1, 0.9, 0.9, 1), nrow = 2, ncol = 2)
> x3 <- rmnorm(n = n, mean = sample.mean, varcov = sample.cov)
> sample.cov <- matrix(data = c(1, -0.9, -0.9, 1), nrow = 2, ncol = 2)
> x4 <- rmnorm(n = n, mean = sample.mean, varcov = sample.cov)
> par(mfrow = c(2, 2))
> plot(x1, pch = 19, col = "blue", xlab = "X", ylab = "Y", main = expression(paste(rho,
+     " = ", 0.6), ))
> plot(x2, pch = 19, col = "blue", xlab = "X", ylab = "Y", main = expression(paste(rho,
+     " = ", -0.6), ))
> plot(x3, pch = 19, col = "blue", xlab = "X", ylab = "Y", main = expression(paste(rho,
+     " = ", 0.9), ))
> plot(x4, pch = 19, col = "blue", xlab = "X", ylab = "Y", main = expression(paste(rho,
+     " = ", -0.9), ))

C2_1-001.jpg

Draw a Location Dispersion Ellipsoid

> sample.cov <- matrix(data = c(1, 0.9, 0.9, 1), nrow = 2, ncol = 2)
> x <- rmnorm(n = n, mean = sample.mean, varcov = sample.cov)
> scaled <- 2
> x1 <- sample.mean[1] - sqrt(sample.cov[1, 1]) * scaled
> x2 <- sample.mean[1] + sqrt(sample.cov[1, 1]) * scaled
> y1 <- sample.mean[2] + sqrt(sample.cov[2, 2]) * scaled
> y2 <- sample.mean[2] - sqrt(sample.cov[2, 2]) * scaled
> par(mfrow = c(1, 1))
> plot(x, pch = 19, col = "blue", xlab = "X", ylab = "Y", main = expression(paste(rho,
+     " = ", 0.9), ), xlim = c(-4, 4), ylim = c(-4, 4))
> vecs <- eigen(sample.cov)$vectors
> evals <- eigen(sample.cov)$values
> length1 <- sqrt(evals[1])
> length2 <- sqrt(evals[2])
> segments(sample.mean[1], sample.mean[2], sample.mean[1] + vecs[1,
+     1] * length1 * scaled, sample.mean[2] + vecs[2, 1] * length1 *
+     scaled, lwd = 3, col = "green")
> segments(sample.mean[1], sample.mean[2], sample.mean[1] + vecs[1,
+     2] * length2 * scaled, sample.mean[2] + vecs[2, 2] * length2 *
+     scaled, lwd = 3, col = "green")
> par(new = T)
> plot(ellipse(0.9), type = "l", col = "yellow", lwd = 4, xlim = c(-4,
+     4), ylim = c(-4, 4))

C2_1-002.jpg

Drawing Ellipse

> angles <- seq(from = 0, to = 2 * pi, by = pi/500)
> ang.points <- matrix(data = NA, nrow = length(angles), ncol = 2)
> for (i in seq_along(angles)) {
+     ang.points[i, 1] <- cos(angles[i]) * sqrt(evals[1])
+     ang.points[i, 2] <- sin(angles[i]) * sqrt(evals[2])
+ }
> plot(ang.points %*% t(vecs), xlab = "", ylab = "", xlim = c(-2,
+     2), ylim = c(-2, 2))

C2_1-003.jpg

Draw a Location Dispersion Ellipsoid with out using ellipse

> sample.cov <- matrix(data = c(1, 0.9, 0.9, 1), nrow = 2, ncol = 2)
> x <- rmnorm(n = n, mean = sample.mean, varcov = sample.cov)
> scaled <- 2
> x1 <- sample.mean[1] - sqrt(sample.cov[1, 1]) * scaled
> x2 <- sample.mean[1] + sqrt(sample.cov[1, 1]) * scaled
> y1 <- sample.mean[2] + sqrt(sample.cov[2, 2]) * scaled
> y2 <- sample.mean[2] - sqrt(sample.cov[2, 2]) * scaled
> par(mfrow = c(1, 1))
> plot(x, pch = 19, col = "blue", xlab = "X", ylab = "Y", main = expression(paste(rho,
+     " = ", 0.9), ), xlim = c(-4, 4), ylim = c(-4, 4))
> vecs <- eigen(sample.cov)$vectors
> evals <- eigen(sample.cov)$values
> length1 <- sqrt(evals[1])
> length2 <- sqrt(evals[2])
> segments(sample.mean[1], sample.mean[2], sample.mean[1] + vecs[1,
+     1] * length1 * scaled, sample.mean[2] + vecs[2, 1] * length1 *
+     scaled, lwd = 3, col = "green")
> segments(sample.mean[1], sample.mean[2], sample.mean[1] + vecs[1,
+     2] * length2 * scaled, sample.mean[2] + vecs[2, 2] * length2 *
+     scaled, lwd = 3, col = "green")
> par(new = T)
> angles <- seq(from = 0, to = 2 * pi, by = pi/500)
> ang.points <- matrix(data = NA, nrow = length(angles), ncol = 2)
> for (i in seq_along(angles)) {
+     ang.points[i, 1] <- cos(angles[i]) * sqrt(evals[1]) * scaled
+     ang.points[i, 2] <- sin(angles[i]) * sqrt(evals[2]) * scaled
+ }
> plot(ang.points %*% t(vecs), xlab = "", ylab = "", xlim = c(-4,
+     4), ylim = c(-4, 4), col = "yellow1")

C2_1-004.jpg