¿Cómo trazar una elipse a partir de valores propios y vectores propios en R? [cerrado]

15

¿Podría alguien idear el código R para trazar una elipse de los valores propios y los vectores propios de la siguiente matriz

UN=(2.2 2.20.4 0.40.4 0.42.8)
MYaseen208
fuente

Respuestas:

16

Puede extraer los vectores propios y los valores mediante eigen(A). Sin embargo, es más simple usar la descomposición de Cholesky. Tenga en cuenta que al trazar elipses de confianza para datos, los ejes de elipse generalmente se escalan para tener una longitud = raíz cuadrada de los valores propios correspondientes, y esto es lo que da la descomposición de Cholesky.

ctr    <- c(0, 0)                               # data centroid -> colMeans(dataMatrix)
A      <- matrix(c(2.2, 0.4, 0.4, 2.8), nrow=2) # covariance matrix -> cov(dataMatrix)
RR     <- chol(A)                               # Cholesky decomposition
angles <- seq(0, 2*pi, length.out=200)          # angles for ellipse
ell    <- 1 * cbind(cos(angles), sin(angles)) %*% RR  # ellipse scaled with factor 1
ellCtr <- sweep(ell, 2, ctr, "+")               # center ellipse to the data centroid
plot(ellCtr, type="l", lwd=2, asp=1)            # plot ellipse
points(ctr[1], ctr[2], pch=4, lwd=2)            # plot data centroid

library(car)  # verify with car's ellipse() function
ellipse(c(0, 0), shape=A, radius=0.98, col="red", lty=2)

Editar: para trazar también los vectores propios, debe usar el enfoque más complicado. Esto es equivalente a la respuesta de suncoolsu, solo usa notación matricial para acortar el código.

eigVal  <- eigen(A)$values
eigVec  <- eigen(A)$vectors
eigScl  <- eigVec %*% diag(sqrt(eigVal))  # scale eigenvectors to length = square-root
xMat    <- rbind(ctr[1] + eigScl[1, ], ctr[1] - eigScl[1, ])
yMat    <- rbind(ctr[2] + eigScl[2, ], ctr[2] - eigScl[2, ])
ellBase <- cbind(sqrt(eigVal[1])*cos(angles), sqrt(eigVal[2])*sin(angles)) # normal ellipse
ellRot  <- eigVec %*% t(ellBase)                                          # rotated ellipse
plot((ellRot+ctr)[1, ], (ellRot+ctr)[2, ], asp=1, type="l", lwd=2)
matlines(xMat, yMat, lty=1, lwd=2, col="green")
points(ctr[1], ctr[2], pch=4, col="red", lwd=3)

ingrese la descripción de la imagen aquí

lince
fuente
¿Te importaría trazar los valores propios y los vectores propios en esta elipse? Gracias
MYaseen208
@ MYaseen208 Edité mi respuesta para mostrar los vectores propios como los ejes de la elipse. La mitad de la longitud de los ejes es igual a la raíz cuadrada de los vectores propios correspondientes.
caracal
7

Creo que este es el código R que deseas. Tomé prestado el código R de este hilo en la lista de correo r. Básicamente, la idea es: los medios diámetros mayor y menor son los dos valores propios y usted gira la elipse por la cantidad de ángulo entre el primer vector propio y el eje x

mat <- matrix(c(2.2, 0.4, 0.4, 2.8), 2, 2)
eigens <- eigen(mat)
evs <- sqrt(eigens$values)
evecs <- eigens$vectors

a <- evs[1]
b <- evs[2]
x0 <- 0
y0 <- 0
alpha <- atan(evecs[ , 1][2] / evecs[ , 1][1])
theta <- seq(0, 2 * pi, length=(1000))

x <- x0 + a * cos(theta) * cos(alpha) - b * sin(theta) * sin(alpha)
y <- y0 + a * cos(theta) * sin(alpha) + b * sin(theta) * cos(alpha)


png("graph.png")
plot(x, y, type = "l", main = expression("x = a cos " * theta * " + " * x[0] * " and y = b sin " * theta * " + " * y[0]), asp = 1)
arrows(0, 0, a * evecs[ , 1][2], a * evecs[ , 1][2])
arrows(0, 0, b * evecs[ , 2][3], b * evecs[ , 2][2])
dev.off()

ingrese la descripción de la imagen aquí

suncoolsu
fuente
Por favor, siéntase libre de corregirme. No creo que los eigen vecs sean perpendiculares (deben ser en teoría; ¿puedo estar tramando algo mal?).
suncoolsu
UN=(1-5 5-5 51)
Simplemente configure asp=1una relación de aspecto de 1 y flechas perpendiculares. Cambiar su código a evs <- sqrt(eigens$values)da la misma elipse que mi respuesta.
caracal
33
@ MYaseen208 Su nueva matriz no es positiva definida: tiene valores propios negativos y no es una matriz de covarianza posible. No sé qué elipse dibujar en ese caso.
caracal
@caracal gracias! ... sí, ¡me perdí la parte sqrt!
suncoolsu